LCOV - code coverage report
Current view: top level - columnphysics - icepack_algae.F90 (source / functions) Hit Total Coverage
Test: 200704-015006:b1e41d9f12:3:base,travis,quick Lines: 897 1175 76.34 %
Date: 2020-07-03 20:11:40 Functions: 10 11 90.91 %

          Line data    Source code
       1             : !=======================================================================
       2             : !
       3             : ! Compute sea ice biogeochemistry (vertical or skeletal layer)
       4             : !
       5             : ! authors: Nicole Jeffery, LANL
       6             : !          Scott Elliot,   LANL
       7             : !          Elizabeth C. Hunke, LANL
       8             : !
       9             :       module icepack_algae
      10             : 
      11             :       use icepack_kinds
      12             : 
      13             :       use icepack_parameters, only: p05, p5, c0, c1, c2, c6, c10
      14             :       use icepack_parameters, only: pi, secday, puny
      15             :       use icepack_parameters, only: hs_ssl, sk_l
      16             : 
      17             :       use icepack_parameters, only: dEdd_algae, solve_zbgc
      18             :       use icepack_parameters, only: R_dFe2dust, dustFe_sol, algal_vel
      19             :       use icepack_parameters, only: bgc_flux_type
      20             :       use icepack_parameters, only: grid_o
      21             :       use icepack_parameters, only: T_max, fsal      , fr_resp
      22             :       use icepack_parameters, only: op_dep_min       , fr_graze_s
      23             :       use icepack_parameters, only: fr_graze_e       , fr_mort2min
      24             :       use icepack_parameters, only: fr_dFe           , k_nitrif
      25             :       use icepack_parameters, only: t_iron_conv      , max_loss
      26             :       use icepack_parameters, only: max_dfe_doc1     , fr_resp_s
      27             :       use icepack_parameters, only: y_sk_DMS         , t_sk_conv
      28             :       use icepack_parameters, only: t_sk_ox
      29             : 
      30             :       use icepack_tracers, only: ntrcr, bio_index 
      31             :       use icepack_tracers, only: nt_bgc_N, nt_fbri, nt_zbgc_frac
      32             :       use icepack_tracers, only: tr_brine
      33             :       use icepack_tracers, only: tr_bgc_Nit,    tr_bgc_Am,    tr_bgc_Sil
      34             :       use icepack_tracers, only: tr_bgc_DMS,    tr_bgc_PON
      35             :       use icepack_tracers, only: tr_bgc_N,      tr_bgc_C,     tr_bgc_chl
      36             :       use icepack_tracers, only: tr_bgc_DON,    tr_bgc_Fe,    tr_zaero
      37             :       use icepack_tracers, only: nlt_bgc_Nit,   nlt_bgc_Am,   nlt_bgc_Sil
      38             :       use icepack_tracers, only: nlt_bgc_DMS,   nlt_bgc_PON
      39             :       use icepack_tracers, only: nlt_bgc_N,     nlt_bgc_C,    nlt_bgc_chl
      40             :       use icepack_tracers, only: nlt_bgc_DOC,   nlt_bgc_DON,  nlt_bgc_DIC
      41             :       use icepack_tracers, only: nlt_zaero  ,   nlt_bgc_DMSPp,nlt_bgc_DMSPd
      42             :       use icepack_tracers, only: nlt_bgc_Fed,   nlt_bgc_Fep
      43             : 
      44             :       use icepack_zbgc_shared, only: remap_zbgc, regrid_stationary
      45             :       use icepack_zbgc_shared, only: merge_bgc_fluxes
      46             :       use icepack_zbgc_shared, only: merge_bgc_fluxes_skl
      47             :       use icepack_zbgc_shared, only: phi_sk, bgc_tracer_type
      48             :       use icepack_zbgc_shared, only: zbgc_init_frac
      49             :       use icepack_zbgc_shared, only: zbgc_frac_init
      50             :       use icepack_zbgc_shared, only: tau_rel, tau_ret, thinS
      51             :       use icepack_zbgc_shared, only: r_Si2N, R_Fe2N, R_S2N, R_chl2N
      52             :       use icepack_zbgc_shared, only: chlabs, alpha2max_low, beta2max, mu_max
      53             :       use icepack_zbgc_shared, only: k_exude, K_Nit, K_Am, K_Sil, K_Fe
      54             :       use icepack_zbgc_shared, only: grow_Tdep, fr_graze, mort_pre, mort_Tdep
      55             :       use icepack_zbgc_shared, only: f_don, kn_bac, f_don_Am
      56             :       use icepack_zbgc_shared, only: f_doc, f_exude, k_bac, R_chl2N, R_C2N
      57             : 
      58             :       use icepack_warnings, only: warnstr, icepack_warnings_add
      59             :       use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
      60             : 
      61             :       use icepack_aerosol, only: update_snow_bgc
      62             : 
      63             :       implicit none
      64             : 
      65             :       private 
      66             :       public :: zbio, sklbio
      67             : 
      68             :       real (kind=dbl_kind), parameter :: & 
      69             :          exp_argmax = c10    ! maximum argument of exponential
      70             : 
      71             : !=======================================================================
      72             : 
      73             :       contains
      74             : 
      75             : !=======================================================================
      76             : 
      77      199620 :       subroutine zbio   (dt,           nblyr,       &
      78             :                          nslyr,        nilyr,       &
      79             :                          meltt,        melts,       &
      80             :                          meltb,        congel,      &
      81             :                          snoice,       nbtrcr,      &
      82             :                          fsnow,        ntrcr,       &
      83      199620 :                          trcrn,        bio_index,   &
      84             :                          aice_old,                  &
      85             :                          vice_old,     vsno_old,    &
      86             :                          vicen,        vsnon,       &
      87      199620 :                          aicen,        flux_bio_atm,& 
      88             :                          n_cat,        n_algae,     & 
      89             :                          n_doc,        n_dic,       &
      90             :                          n_don,                     &
      91             :                          n_fed,        n_fep,       &
      92             :                          n_zaero,      first_ice,   &
      93      199620 :                          hice_old,     ocean_bio,   & 
      94      199620 :                          bphin,        iphin,       &
      95      199620 :                          iDin, &
      96      199620 :                          fswthrul,                  &
      97             :                          dh_top,       dh_bot,      &
      98      199620 :                          zfswin,                    & 
      99             :                          hbri,         hbri_old,    &
     100             : !                        darcy_V,      darcy_V_chl, &
     101             :                          darcy_V, &
     102      199620 :                          bgrid,       &
     103      199620 :                          igrid,        icgrid,      &
     104             :                          bphi_min,                  &
     105      199620 :                          iTin,        &
     106      199620 :                          Zoo,                       &
     107      199620 :                          flux_bio,     dh_direct,   &
     108             :                          upNO,         upNH,        &
     109      199620 :                          fbio_snoice,  fbio_atmice, &
     110      199620 :                          PP_net,       ice_bio_net, &
     111      199620 :                          snow_bio_net, grow_net     )
     112             : 
     113             :       integer (kind=int_kind), intent(in) :: &
     114             :          nblyr,              & ! number of bio layers
     115             :          nslyr,              & ! number of snow layers
     116             :          nilyr,              & ! number of ice layers
     117             :          nbtrcr,             & ! number of distinct bio tracers
     118             :          n_cat,              & ! category number
     119             :          n_algae,            & ! number of autotrophs
     120             :          n_zaero,            & ! number of aerosols
     121             :          n_doc, n_dic,  n_don, n_fed, n_fep, &
     122             :          ntrcr                 ! number of tracers
     123             : 
     124             :       integer (kind=int_kind), dimension (nbtrcr), intent(in) :: &
     125             :          bio_index       
     126             : 
     127             :       real (kind=dbl_kind), intent(in) :: &
     128             :          dt,       &  ! time step
     129             :          hbri,     &  ! brine height  (m)
     130             :          bphi_min, &  ! surface porosity
     131             :          meltt,    &  ! thermodynamic melt/growth rates in dt (m)
     132             :          melts,    &
     133             :          meltb,    &
     134             :          congel,   &
     135             :          snoice,   &
     136             :          fsnow,    & ! snowfall rate (kg/m^2 s)
     137             :          hice_old, & ! ice height (m)
     138             :          vicen,    & ! ice volume (m)
     139             :          vsnon,    & ! snow volume (m)
     140             :          aicen,    & ! ice area fraction
     141             :          aice_old, & ! values prior to thermodynamic changes
     142             :          vice_old, &
     143             :          vsno_old, &
     144             :          darcy_V,  & ! darcy velocity
     145             : !        darcy_V_chl,& ! darcy velocity for algae
     146             :          dh_bot,     & ! change in brine bottom (m)
     147             :          dh_top,     & ! change in brine top (m)
     148             :          dh_direct     ! surface flooding or surface runoff (m)
     149             : 
     150             :       real (kind=dbl_kind), dimension (nbtrcr), intent(inout) :: &
     151             :          snow_bio_net,& ! net bio tracer in snow (mmol/m^2)
     152             :          ice_bio_net, & ! net bio tracer in ice (mmol/m^2)
     153             :          fbio_atmice, & ! bio flux from atm to ice (mmol/m^2/s)
     154             :          fbio_snoice, & ! bio flux from snow to ice  (mmol/m^2/s)
     155             :          flux_bio       ! total ocean tracer flux (mmol/m^2/s)
     156             : 
     157             :       real (kind=dbl_kind), intent(inout) :: &
     158             :          hbri_old       ! brine height  (m)
     159             : 
     160             :       real (kind=dbl_kind), dimension (nblyr+2), intent(inout) :: &
     161             :          bgrid          ! biology nondimensional vertical grid points
     162             : 
     163             :       real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: &
     164             :          igrid      , & ! biology vertical interface points
     165             :          iTin       , & ! salinity vertical interface points
     166             :          iphin      , & ! Porosity on the igrid   
     167             :          iDin           ! Diffusivity/h on the igrid (1/s)
     168             :  
     169             :       real (kind=dbl_kind), dimension (nilyr+1), intent(in) :: &
     170             :          icgrid     , & ! CICE interface coordinate   
     171             :          fswthrul       ! visible short wave radiation on icgrid (W/m^2)  
     172             : 
     173             :       real (kind=dbl_kind), dimension(nbtrcr), &
     174             :          intent(in) :: &
     175             :          flux_bio_atm   ! aerosol/bgc deposition rate (mmol/m^2 s)
     176             : 
     177             :       real (kind=dbl_kind), dimension(ntrcr), &
     178             :          intent(inout) :: &
     179             :          trcrn 
     180             : 
     181             :       real (kind=dbl_kind), dimension (nblyr+1), intent(inout) :: & 
     182             :          zfswin         ! visible Short wave flux on igrid (W/m^2)  
     183             :        
     184             :       real (kind=dbl_kind), dimension (nblyr+1), intent(inout) :: & 
     185             :          Zoo            ! N losses to the system from reaction terms
     186             :                         ! (ie. zooplankton/bacteria) (mmol/m^3)  
     187             : 
     188             :       real (kind=dbl_kind), dimension (nbtrcr), intent(in) :: &   
     189             :          !change to  inout when updating ocean fields
     190             :          ocean_bio      ! ocean concentrations (mmol/m^3) 
     191             : 
     192             :       real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
     193             :          bphin          ! Porosity on the bgrid
     194             : 
     195             :       real (kind=dbl_kind), intent(inout):: & 
     196             :          PP_net     , & ! net PP (mg C/m^2/d)  times aice
     197             :          grow_net   , & ! net specific growth (m/d) times vice
     198             :          upNO       , & ! tot nitrate uptake rate (mmol/m^2/d) times aice 
     199             :          upNH           ! tot ammonium uptake rate (mmol/m^2/d) times aice
     200             : 
     201             :       logical (kind=log_kind), intent(in) :: &
     202             :          first_ice      ! initialized values should be used
     203             : 
     204             :       ! local variables
     205             : 
     206             :       integer (kind=int_kind) :: &
     207             :          mm              ! thickness category index
     208             : 
     209             :       real (kind=dbl_kind), dimension (nblyr+1,n_algae) :: &
     210     1649268 :          upNOn      , & ! algal nitrate uptake rate  (mmol/m^3/s)
     211     1649268 :          upNHn      , & ! algal ammonium uptake rate (mmol/m^3/s)
     212     1649268 :          grow_alg       ! algal growth rate          (mmol/m^3/s)
     213             : 
     214             :       real (kind=dbl_kind), dimension (nbtrcr) :: &
     215     1279678 :          flux_bion       !tracer flux to ocean
     216             : 
     217             :       real (kind=dbl_kind),dimension(nbtrcr) :: &
     218     1279678 :          zbgc_snown, & ! aerosol contribution from snow to ice
     219     1279678 :          zbgc_atmn     ! and atm to ice concentration * volume (mmol/m^3*m)
     220             : 
     221             :       real (kind=dbl_kind), dimension(nbtrcr) :: &
     222     1279678 :          Tot_BGC_i, & ! initial column sum, ice + snow,  of BGC tracer (mmol/m^2)
     223     1279678 :          Tot_BGC_f, & ! final column sum
     224     1224292 :          flux_bio_sno !
     225             : 
     226             :       real (kind=dbl_kind) :: &
     227       48078 :          hsnow_i,  & ! initial snow thickness (m)
     228       48078 :          hsnow_f     ! final snow thickness (m)
     229             : 
     230             :       logical (kind=log_kind) :: &
     231             :          write_flux_diag
     232             : 
     233             :       real (kind=dbl_kind) :: &
     234       48078 :          a_ice
     235             : 
     236             :       character(len=*),parameter :: subname='(zbio)'
     237             : 
     238     4158976 :       zbgc_snown(:) = c0
     239     4158976 :       zbgc_atmn (:) = c0
     240     4158976 :       flux_bion (:) = c0
     241     4158976 :       flux_bio_sno(:) = c0
     242     4158976 :       Tot_BGC_i (:) = c0
     243     4158976 :       Tot_BGC_f (:) = c0
     244      199620 :       hsnow_i = c0
     245      199620 :       hsnow_f = c0
     246      199620 :       write_flux_diag = .false.
     247             :     
     248      199620 :       if (write_flux_diag) then
     249           0 :          if (aice_old > c0) then
     250           0 :             hsnow_i = vsno_old/aice_old
     251           0 :             do  mm = 1,nbtrcr
     252             :                call bgc_column_sum (nblyr, nslyr, hsnow_i, hbri_old, &
     253           0 :                               trcrn(bio_index(mm):bio_index(mm)+nblyr+2), &
     254           0 :                               Tot_BGC_i(mm))
     255           0 :                if (icepack_warnings_aborted(subname)) return
     256             :             enddo
     257             :          endif
     258             :       endif
     259             :  
     260             :       call update_snow_bgc     (dt,        nblyr,        &
     261             :                                 nslyr,                   &
     262             :                                 meltt,     melts,        &
     263             :                                 meltb,     congel,       &
     264             :                                 snoice,    nbtrcr,       &
     265             :                                 fsnow,     ntrcr,        &
     266             :                                 trcrn,     bio_index,    &
     267             :                                 aice_old,  zbgc_snown,   &
     268             :                                 vice_old,  vsno_old,     &
     269             :                                 vicen,     vsnon,        &
     270             :                                 aicen,     flux_bio_atm, &
     271      199620 :                                 zbgc_atmn, flux_bio_sno)
     272      199620 :       if (icepack_warnings_aborted(subname)) return
     273             : 
     274             :       call z_biogeochemistry   (n_cat,        dt,        &
     275             :                                 nilyr,        &
     276             :                                 nblyr,        nbtrcr,    &
     277             :                                 n_algae,      n_doc,     & 
     278             :                                 n_dic,        n_don,     &
     279             :                                 n_fed,        n_fep,     &
     280             :                                 n_zaero,      first_ice, &
     281             :                                 aicen,        vicen,     & 
     282           0 :                                 hice_old,     ocean_bio, & 
     283           0 :                                 flux_bion,    bphin,     &
     284           0 :                                 iphin,        trcrn,     &  
     285           0 :                                 iDin,  &
     286           0 :                                 fswthrul,     grow_alg,  &
     287           0 :                                 upNOn,        upNHn,     &
     288             :                                 dh_top,       dh_bot,    &
     289           0 :                                 zfswin,       hbri,      & 
     290             :                                 hbri_old,     darcy_V,   &
     291             : !                               darcy_V_chl,  bgrid,     &
     292           0 :                                 bgrid,     &
     293           0 :                                 igrid,        icgrid,    &
     294           0 :                                 bphi_min,     zbgc_snown,&
     295           0 :                                 zbgc_atmn, &
     296           0 :                                 iTin,         dh_direct, &
     297           0 :                                 Zoo,          meltb,     &
     298      199620 :                                 congel                   )
     299      199620 :       if (icepack_warnings_aborted(subname)) return
     300             :       
     301     4158976 :       do mm = 1,nbtrcr
     302     4158976 :          flux_bion(mm) = flux_bion(mm) + flux_bio_sno(mm)
     303             :       enddo
     304             : 
     305      199620 :       if (write_flux_diag) then
     306           0 :          if (aicen > c0) then
     307           0 :             hsnow_f = vsnon/aicen
     308           0 :             do mm = 1,nbtrcr
     309             :                call bgc_column_sum (nblyr, nslyr, hsnow_f, hbri, &
     310           0 :                               trcrn(bio_index(mm):bio_index(mm)+nblyr+2), &
     311           0 :                               Tot_BGC_f(mm))
     312           0 :                write(warnstr,*) subname, 'mm, Tot_BGC_i(mm), Tot_BGC_f(mm)'
     313           0 :                call icepack_warnings_add(warnstr)
     314           0 :                write(warnstr,*) subname,  mm, Tot_BGC_i(mm), Tot_BGC_f(mm)
     315           0 :                call icepack_warnings_add(warnstr)
     316           0 :                write(warnstr,*) subname, 'flux_bion(mm), flux_bio_atm(mm)'
     317           0 :                call icepack_warnings_add(warnstr)
     318           0 :                write(warnstr,*) subname,  flux_bion(mm), flux_bio_atm(mm)
     319           0 :                call icepack_warnings_add(warnstr)
     320           0 :                write(warnstr,*) subname, 'zbgc_snown(mm),zbgc_atmn(mm)'
     321           0 :                call icepack_warnings_add(warnstr)
     322           0 :                write(warnstr,*) subname,  zbgc_snown(mm),zbgc_atmn(mm)
     323           0 :                call icepack_warnings_add(warnstr)
     324           0 :                write(warnstr,*) subname, 'Tot_BGC_i(mm) + flux_bio_atm(mm)*dt - flux_bion(mm)*dt'
     325           0 :                call icepack_warnings_add(warnstr)
     326           0 :                write(warnstr,*) subname,  Tot_BGC_i(mm) + flux_bio_atm(mm)*dt - flux_bion(mm)*dt
     327           0 :                call icepack_warnings_add(warnstr)
     328             :             enddo
     329             :          endif
     330             :       endif
     331             : 
     332      199620 :       if (icepack_warnings_aborted(subname)) return
     333             : 
     334             :       call merge_bgc_fluxes   (dt,           nblyr,      &
     335           0 :                                bio_index,    n_algae,    &
     336             :                                nbtrcr,       aicen,      &    
     337             :                                vicen,        vsnon,      &
     338           0 :                                iphin,      &
     339           0 :                                trcrn,                    &
     340           0 :                                flux_bion,    flux_bio,   &
     341           0 :                                upNOn,        upNHn,      &
     342             :                                upNO,         upNH,       &
     343           0 :                                zbgc_snown,   zbgc_atmn,  &
     344           0 :                                fbio_snoice,  fbio_atmice,&
     345           0 :                                PP_net,       ice_bio_net,&
     346           0 :                                snow_bio_net, grow_alg,   &
     347      199620 :                                grow_net)
     348      199620 :       if (icepack_warnings_aborted(subname)) return
     349             :  
     350      199620 :       if (write_flux_diag) then
     351           0 :          if (aicen > c0) then
     352           0 :             if (n_cat .eq. 1) a_ice = c0
     353           0 :             a_ice = a_ice + aicen
     354           0 :             write(warnstr,*) subname, 'after merge_bgc_fluxes, n_cat:', n_cat
     355           0 :             call icepack_warnings_add(warnstr)
     356           0 :             do mm = 1,nbtrcr
     357           0 :                write(warnstr,*) subname,  'mm, flux_bio(mm):',mm,flux_bio(mm)
     358           0 :                call icepack_warnings_add(warnstr)
     359           0 :                write(warnstr,*) subname, 'fbio_snoice(mm)',fbio_snoice(mm)
     360           0 :                call icepack_warnings_add(warnstr)
     361           0 :                write(warnstr,*) subname, 'fbio_atmice(mm)',fbio_atmice(mm)
     362           0 :                call icepack_warnings_add(warnstr)
     363           0 :                write(warnstr,*) subname,  'flux_bio_atm(mm)', flux_bio_atm(mm)
     364           0 :                call icepack_warnings_add(warnstr)
     365           0 :                write(warnstr,*) subname,  'flux_bio_atm(mm)*a_ice', flux_bio_atm(mm)*a_ice
     366           0 :                call icepack_warnings_add(warnstr)
     367             :             enddo
     368             :          endif
     369             :       endif
     370             : 
     371             :       end subroutine zbio    
     372             : 
     373             : !=======================================================================
     374             : 
     375        8261 :       subroutine sklbio       (dt,       ntrcr,      &
     376             :                                nbtrcr,   n_algae,    &
     377             :                                n_doc,      &
     378             :                                n_dic,    n_don,      &
     379             :                                n_fed,    n_fep,      &
     380       16522 :                                flux_bio, ocean_bio,  &
     381             :                                aicen,      &
     382             :                                meltb,    congel,     &
     383             :                                fswthru,  first_ice,  &
     384        8261 :                                trcrn,  &
     385             :                                PP_net,   upNO,       &
     386             :                                upNH,     grow_net    )
     387             : 
     388             :       integer (kind=int_kind), intent(in) :: &
     389             :          nbtrcr,             & ! number of distinct bio tracers
     390             :          n_algae,            & ! number of autotrophs
     391             :          n_doc, n_dic,       & ! number of dissolved organic, inorganic carbon
     392             :          n_don,              & ! number of dissolved organic nitrogen
     393             :          n_fed, n_fep,       & ! number of iron
     394             :          ntrcr                 ! number of tracers
     395             : 
     396             :       logical (kind=log_kind), intent(in) :: &
     397             :          first_ice      ! initialized values should be used
     398             : 
     399             :       real (kind=dbl_kind), intent(in) :: &
     400             :          dt,       &  ! time step
     401             : !        hmix,     &  ! mixed layer depth (m)
     402             :          aicen,    &  ! ice area fraction
     403             :          meltb,    &  ! bottom melt (m)
     404             :          congel,   &  ! bottom growth (m)
     405             :          fswthru      ! visible shortwave passing to ocean(W/m^2)  
     406             : 
     407             :       real (kind=dbl_kind), dimension(ntrcr), intent(inout) :: &
     408             :          trcrn      ! bulk concentration per m^3
     409             : 
     410             :       real (kind=dbl_kind), dimension (nbtrcr), intent(inout) :: &
     411             :          flux_bio   ! ocean tracer flux (mmol/m^2/s) positive into ocean
     412             : 
     413             :       real (kind=dbl_kind), dimension (nbtrcr), intent(in) :: &
     414             :          ocean_bio  ! ocean tracer concentration (mmol/m^3)
     415             : 
     416             :       ! history output
     417             :       real (kind=dbl_kind), intent(inout):: & 
     418             :          PP_net  , & ! Bulk net PP (mg C/m^2/s)
     419             :          grow_net, & ! net specific growth (/s)
     420             :          upNO    , & ! tot nitrate uptake rate (mmol/m^2/s) 
     421             :          upNH        ! tot ammonium uptake rate (mmol/m^2/s)
     422             : 
     423             :       ! local variables
     424             : 
     425             :       real (kind=dbl_kind), dimension (n_algae) :: &
     426       33044 :          upNOn      , & ! algal nitrate uptake rate  (mmol/m^3/s)
     427       33044 :          upNHn      , & ! algal ammonium uptake rate (mmol/m^3/s)
     428       33044 :          grow_alg       ! algal growth rate          (mmol/m^3/s)
     429             : 
     430             :       real (kind=dbl_kind), dimension (nbtrcr) :: &
     431      156959 :          flux_bion       !tracer flux to ocean
     432             : 
     433             :       character(len=*),parameter :: subname='(sklbio)'
     434             : 
     435      140437 :       flux_bion (:) = c0
     436       33044 :       upNOn     (:) = c0
     437       33044 :       upNHn     (:) = c0
     438       33044 :       grow_alg  (:) = c0
     439             : 
     440             :       call skl_biogeochemistry       (dt, &
     441             :                                       n_doc,     &
     442             :                                       n_dic,     n_don,     &
     443             :                                       n_fed,     n_fep,     &
     444             :                                       nbtrcr,    n_algae,   &
     445           0 :                                       flux_bion, ocean_bio, &
     446             : !                                     hmix,      aicen,     &
     447             :                                       meltb,     congel,    &
     448             :                                       fswthru,   first_ice, &
     449           0 :                                       trcrn,     upNOn,     &
     450        8261 :                                       upNHn,     grow_alg)
     451        8261 :       if (icepack_warnings_aborted(subname)) return
     452             : 
     453             :       call merge_bgc_fluxes_skl    (nbtrcr,    n_algae,     &
     454           0 :                                     aicen,     trcrn,       &
     455           0 :                                     flux_bion, flux_bio,    &
     456           0 :                                     PP_net,    upNOn,       &
     457           0 :                                     upNHn,     upNO,        &
     458             :                                     upNH,      grow_net,    &
     459        8261 :                                     grow_alg)
     460        8261 :       if (icepack_warnings_aborted(subname)) return
     461             :  
     462             :       end subroutine sklbio    
     463             : 
     464             : !=======================================================================
     465             : !
     466             : ! skeletal layer biochemistry
     467             : ! 
     468        8261 :       subroutine skl_biogeochemistry (dt, &
     469             :                                       n_doc,        &
     470             :                                       n_dic,      n_don,        &
     471             :                                       n_fed,      n_fep,        &
     472             :                                       nbtrcr,     n_algae,      &
     473       16522 :                                       flux_bio,   ocean_bio,    &
     474             : !                                     hmix,       aicen,        &
     475             :                                       meltb,      congel,       &
     476             :                                       fswthru,    first_ice,    &
     477       16522 :                                       trcrn,      upNOn,        &
     478       16522 :                                       upNHn,      grow_alg_skl)
     479             : 
     480             :       integer (kind=int_kind), intent(in) :: &
     481             :          n_doc, n_dic,  n_don, n_fed, n_fep, &
     482             :          nbtrcr , n_algae      ! number of bgc tracers and number algae
     483             : 
     484             :       real (kind=dbl_kind), intent(in) :: &
     485             :          dt     , & ! time step 
     486             : !        hmix   , & ! mixed layer depth
     487             : !        aicen  , & ! ice area 
     488             :          meltb  , & ! bottom ice melt
     489             :          congel , & ! bottom ice growth 
     490             :          fswthru    ! shortwave passing through ice to ocean
     491             : 
     492             :       logical (kind=log_kind), intent(in) :: &
     493             :          first_ice  ! initialized values should be used
     494             : 
     495             :       real (kind=dbl_kind), dimension(:), intent(inout) :: &
     496             :          trcrn      ! bulk concentration per m^3
     497             :        
     498             :       ! history variables
     499             :   
     500             :       real (kind=dbl_kind), dimension (:), intent(out) :: &
     501             :          flux_bio   ! ocean tracer flux (mmol/m^2/s) positive into ocean
     502             : 
     503             :       real (kind=dbl_kind), dimension (:), intent(in) :: &
     504             :          ocean_bio  ! ocean tracer concentration (mmol/m^3)
     505             : 
     506             :       real (kind=dbl_kind), dimension (:), intent(out) :: &
     507             :          grow_alg_skl, & ! tot algal growth rate (mmol/m^3/s)  
     508             :          upNOn       , & !  algal NO uptake rate (mmol/m^3/s) 
     509             :          upNHn           !  algal NH uptake rate (mmol/m^3/s) 
     510             : 
     511             :       ! local variables
     512             : 
     513             :       integer (kind=int_kind) :: nn
     514             : 
     515             :       real (kind=dbl_kind), dimension(nbtrcr):: &
     516      140437 :          react        , & ! biological sources and sinks (mmol/m^3)
     517      156959 :          cinit        , & ! initial brine concentration*sk_l (mmol/m^2)
     518      140437 :          cinit_v      , & ! initial brine concentration (mmol/m^3)
     519      140437 :          congel_alg   , & ! congelation flux contribution to ice algae (mmol/m^2 s) 
     520             :                           ! (used as initialization)
     521      140437 :          f_meltn      , & ! vertical melt fraction of skeletal layer in dt
     522      140437 :          flux_bio_temp, & ! tracer flux to ocean (mmol/m^2 s)
     523      140437 :          PVflag       , & ! 1 for tracers that flow with the brine, 0 otherwise
     524      148698 :          cling            ! 1 for tracers that cling, 0 otherwise
     525             : 
     526             :       real (kind=dbl_kind) :: &
     527        8261 :          Zoo_skl      , & ! N losses from zooplankton/bacteria ... (mmol/m^3)
     528        8261 :          iTin         , &
     529        8261 :          PVt          , & ! type 'Jin2006' piston velocity (m/s) 
     530        8261 :          ice_growth   , & ! Jin2006 definition: either congel rate or bottom melt rate  (m/s)
     531        8261 :          grow_val     , & ! (m/x)
     532        8261 :          rphi_sk      , & ! 1 / skeletal layer porosity
     533        8261 :          cinit_tmp    , & ! temporary variable for concentration (mmol/m^2)
     534        8261 :          Nerror           ! change in total nitrogen from reactions
     535             : 
     536             :       real (kind=dbl_kind), parameter :: &
     537             :          PVc = 1.e-6_dbl_kind           , & ! type 'constant' piston velocity for interface (m/s) 
     538             :          PV_scale_growth = p5           , & ! scale factor in Jin code PV during ice growth
     539             :          PV_scale_melt = p05            , & ! scale factor in Jin code PV during ice melt
     540             :          growth_max = 1.85e-5_dbl_kind , & ! PVt function reaches maximum here.  (m/s)
     541             :          Tin_bot = -1.8_dbl_kind        , & ! temperature of the ice bottom (oC)
     542             :          MJ1 = 9.667e-9_dbl_kind        , & ! (m/s) coefficients in Jin2008
     543             :          MJ2 = 38.8_dbl_kind            , & ! (1) from:4.49e-4_dbl_kind*secday   
     544             :          MJ3 = 1.04e7_dbl_kind          , & ! 1/(m/s) from: 1.39e-3_dbl_kind*secday^2  
     545             :          PV_frac_max = 0.9_dbl_kind         ! Maximum Piston velocity is 90% of skeletal layer/dt
     546             : 
     547             :       logical (kind=log_kind) :: conserve_N
     548             : 
     549             :       character(len=*),parameter :: subname='(skl_biogeochemistry)'
     550             : 
     551             :       !-----------------------------------------------------------------
     552             :       ! Initialize 
     553             :       !-----------------------------------------------------------------
     554             : 
     555        8261 :       conserve_N = .true.
     556        8261 :       Zoo_skl    = c0
     557        8261 :       rphi_sk    = c1/phi_sk
     558        8261 :       PVt        = c0
     559        8261 :       iTin       = Tin_bot
     560             : 
     561      140437 :       do nn = 1, nbtrcr 
     562      132176 :          cinit     (nn) = c0
     563      132176 :          cinit_v   (nn) = c0 
     564      132176 :          congel_alg(nn) = c0
     565      132176 :          f_meltn   (nn) = c0
     566      132176 :          react     (nn) = c0
     567      132176 :          PVflag    (nn) = c1
     568      132176 :          cling     (nn) = c0
     569             : 
     570             :       !-----------------------------------------------------------------
     571             :       ! only the dominant tracer_type affects behavior
     572             :       !  < 0 is purely mobile:  > 0 stationary behavior
     573             :       ! NOTE: retention times are not used in skl model
     574             :       !-----------------------------------------------------------------
     575             : 
     576      132176 :          if (bgc_tracer_type(nn) >= c0) then  
     577       90871 :             PVflag(nn) = c0
     578       90871 :             cling (nn) = c1
     579             :          endif
     580             :  
     581      132176 :          ice_growth = (congel-meltb)/dt
     582      132176 :          if (first_ice) then 
     583         320 :             trcrn(bio_index(nn)) = ocean_bio(nn)   ! * sk_l*rphi_sk
     584             :          endif ! first_ice
     585      132176 :          cinit  (nn) = trcrn(bio_index(nn)) * sk_l * rphi_sk
     586      132176 :          cinit_v(nn) = cinit(nn)/sk_l
     587      140437 :          if (cinit(nn) < c0) then
     588           0 :             write(warnstr,*) subname,'initial sk_bgc < 0, nn,nbtrcr,cinit(nn)', &
     589           0 :                  nn,nbtrcr,cinit(nn)
     590           0 :             call icepack_warnings_add(warnstr)
     591           0 :             call icepack_warnings_add(subname//' cinit < c0')
     592           0 :             call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
     593           0 :             return
     594             :          endif  
     595             :       enddo     ! nbtrcr
     596             : 
     597        8261 :       if (icepack_warnings_aborted(subname)) return
     598             : 
     599        8261 :       if (trim(bgc_flux_type) == 'Jin2006') then  
     600             :  
     601             :       !-----------------------------------------------------------------
     602             :       ! 'Jin2006':
     603             :       ! 1. congel/melt dependent piston velocity (PV) for growth and melt
     604             :       ! 2. If congel > melt use 'congel'; if melt > congel use 'melt'
     605             :       ! 3. For algal N, PV for ice growth only provides a seeding concentration 
     606             :       ! 4. Melt affects nutrients and algae in the same manner through PV(melt)
     607             :       !-----------------------------------------------------------------
     608             : 
     609        8261 :          if (ice_growth > c0) then  ! ice_growth = congel/dt
     610        6266 :             grow_val = min(ice_growth,growth_max)  
     611             :             PVt = -min(abs(PV_scale_growth*(MJ1 + MJ2*grow_val &
     612             :                                                 - MJ3*grow_val**2)), &
     613        6266 :                            PV_frac_max*sk_l/dt)  
     614             :          else                       ! ice_growth = -meltb/dt
     615             :             PVt =  min(abs(PV_scale_melt  *(      MJ2*ice_growth &
     616             :                                                 - MJ3*ice_growth**2)), &
     617        1995 :                            PV_frac_max*sk_l/dt)
     618             :          endif
     619      140437 :          do nn = 1, nbtrcr  
     620      140437 :             if (bgc_tracer_type(nn) >= c0) then
     621       90871 :                if (ice_growth < c0) then ! flux from ice to ocean
     622             :                   ! Algae and clinging tracers melt like nutrients              
     623       21945 :                   f_meltn(nn) = PVt*cinit_v(nn) ! for algae only
     624      137852 :                elseif (ice_growth > c0 .AND. &
     625       68926 :                   cinit(nn) < ocean_bio(nn)*sk_l/phi_sk) then
     626             :                   ! Growth only contributes to seeding from ocean 
     627       36873 :                   congel_alg(nn) = (ocean_bio(nn)*sk_l/phi_sk - cinit(nn))/dt
     628             :                endif ! PVt > c0  
     629             :             endif     
     630             :          enddo
     631             : 
     632             :       else   ! bgc_flux_type = 'constant'
     633             : 
     634             :       !-----------------------------------------------------------------
     635             :       ! 'constant':
     636             :       ! 1. Constant PV for congel > melt
     637             :       ! 2. For algae, PV for ice growth only provides a seeding concentration 
     638             :       ! 3. Melt loss (f_meltn) affects algae only and is proportional to melt
     639             :       !-----------------------------------------------------------------
     640             : 
     641           0 :          if (ice_growth > c0) PVt = -PVc
     642           0 :          do nn = 1, nbtrcr
     643           0 :             if (bgc_tracer_type(nn) >= c0 ) then
     644           0 :                if (ice_growth >= c0 .AND. cinit_v(nn) < ocean_bio(nn)/phi_sk) then
     645           0 :                   congel_alg(nn) = (ocean_bio(nn)*sk_l/phi_sk - cinit(nn))/dt
     646           0 :                elseif (ice_growth < c0) then
     647           0 :                   f_meltn(nn) = min(c1, meltb/sk_l)*cinit(nn)/dt
     648             :                endif
     649             :             endif
     650             :          enddo ! nn
     651             : 
     652             :       endif  ! bgc_flux_type
     653             : 
     654             :       !-----------------------------------------------------------------
     655             :       ! begin building biogeochemistry terms
     656             :       !-----------------------------------------------------------------
     657             : 
     658      140437 :       react(:) = c0  
     659       33044 :       grow_alg_skl(:) = c0
     660             : 
     661             :       call algal_dyn (dt,              &
     662             :                       n_doc, n_dic,  n_don, n_fed, n_fep, &
     663             :                       dEdd_algae, &
     664           0 :                       fswthru,         react,     & 
     665           0 :                       cinit_v, &
     666           0 :                       grow_alg_skl(:),    n_algae,   &
     667             :                       iTin,                       &
     668           0 :                       upNOn(:),           upNHn(:),     &
     669             :                       Zoo_skl,                    &
     670        8261 :                       Nerror,          conserve_N)
     671        8261 :       if (icepack_warnings_aborted(subname)) return
     672             : 
     673             :       !-----------------------------------------------------------------
     674             :       ! compute new tracer concencentrations
     675             :       !-----------------------------------------------------------------
     676             :   
     677      140437 :       do nn = 1, nbtrcr
     678             : 
     679             :       !-----------------------------------------------------------------
     680             :       ! if PVt > 0, ie melt, then ocean_bio term drops out (MJ2006)
     681             :       ! Combine boundary fluxes
     682             :       !-----------------------------------------------------------------
     683             :            
     684      132176 :          PVflag(nn) = SIGN(PVflag(nn),PVt)
     685      132176 :          cinit_tmp = max(c0, cinit_v(nn) + react(nn))
     686      132176 :          flux_bio_temp(nn) = (PVflag(nn)*PVt*cinit_tmp &
     687      264352 :                            -  PVflag(nn)*min(c0,PVt)*ocean_bio(nn)) &
     688      264352 :                            + f_meltn(nn)*cling(nn) - congel_alg(nn)
     689             : 
     690      132176 :          if (cinit_tmp*sk_l < flux_bio_temp(nn)*dt) then
     691           0 :             flux_bio_temp(nn) = cinit_tmp*sk_l/dt*(c1-puny)
     692             :          endif
     693             : 
     694      132176 :          cinit(nn) = cinit_tmp*sk_l - flux_bio_temp(nn)*dt
     695      132176 :          flux_bio(nn) = flux_bio(nn) + flux_bio_temp(nn)*phi_sk  
     696             : 
     697             :          ! Uncomment to update ocean concentration
     698             :          ! Currently not coupled with ocean biogeochemistry
     699             : !         ocean_bio(nn) = ocean_bio(nn) + flux_bio(nn)/hmix*aicen
     700             : 
     701      132176 :          if (.not. conserve_N) then
     702           0 :               write(warnstr,*) subname, 'N not conserved in skl_bgc, Nerror:',Nerror
     703           0 :               call icepack_warnings_add(warnstr)
     704           0 :               write(warnstr,*) subname, 'sk_bgc < 0 after algal fluxes, nn,cinit,flux_bio',&
     705           0 :                                nn,cinit(nn),flux_bio(nn)
     706           0 :               call icepack_warnings_add(warnstr)
     707           0 :               write(warnstr,*) subname, 'cinit_tmp,flux_bio_temp,f_meltn,congel_alg,PVt,PVflag: '
     708           0 :               call icepack_warnings_add(warnstr)
     709           0 :               write(warnstr,*) subname, cinit_tmp,flux_bio_temp(nn),f_meltn(nn), &
     710           0 :                                congel_alg(nn),PVt,PVflag(nn)
     711           0 :               call icepack_warnings_add(warnstr)
     712           0 :               write(warnstr,*) subname, 'congel, meltb: ',congel,meltb
     713           0 :               call icepack_warnings_add(warnstr)
     714           0 :               call icepack_warnings_add(subname//' N not conserved in skl_bgc')
     715           0 :               call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
     716      132176 :          elseif (cinit(nn) < c0) then
     717           0 :               write(warnstr,*) subname, 'sk_bgc < 0 after algal fluxes, nn,cinit,flux_bio',&
     718           0 :                                nn,cinit(nn),flux_bio(nn)
     719           0 :               call icepack_warnings_add(warnstr)
     720           0 :               write(warnstr,*) subname, 'cinit_tmp,flux_bio_temp,f_meltn,congel_alg,PVt,PVflag: '
     721           0 :               call icepack_warnings_add(warnstr)
     722           0 :               write(warnstr,*) subname, cinit_tmp,flux_bio_temp(nn),f_meltn(nn), &
     723           0 :                                congel_alg(nn),PVt,PVflag(nn)
     724           0 :               call icepack_warnings_add(warnstr)
     725           0 :               write(warnstr,*) subname, 'congel, meltb: ',congel,meltb
     726           0 :               call icepack_warnings_add(warnstr)
     727           0 :               call icepack_warnings_add(subname//'sk_bgc < 0 after algal fluxes')
     728           0 :               call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
     729             :          endif
     730             : 
     731      132176 :          if (icepack_warnings_aborted(subname)) return
     732             :          
     733             :       !-----------------------------------------------------------------
     734             :       ! reload tracer array:  Bulk tracer concentration (mmol or mg per m^3)
     735             :       !-----------------------------------------------------------------
     736             : 
     737      140437 :          trcrn(bio_index(nn)) = cinit(nn) * phi_sk/sk_l
     738             :         
     739             :        enddo  !nbtrcr  
     740             : 
     741             :       end subroutine skl_biogeochemistry
     742             : 
     743             : !=======================================================================
     744             : !
     745             : ! Solve the scalar vertical diffusion equation implicitly using 
     746             : ! tridiag_solver. Calculate the diffusivity from temperature and salinity.
     747             : ! 
     748             : ! NOTE: In this subroutine, trcrn(nt_fbri) is  the volume fraction of ice with 
     749             : ! dynamic salinity or the height ratio == hinS/vicen*aicen, where hinS is the 
     750             : ! height of the brine surface relative to the bottom of the ice.  This volume fraction
     751             : ! may be > 1 in which case there is brine above the ice surface (meltponds). 
     752             : ! 
     753             : 
     754      199620 :       subroutine z_biogeochemistry (n_cat,        dt,        &
     755             :                                     nilyr,        &
     756             :                                     nblyr,        nbtrcr,    &
     757             :                                     n_algae,      n_doc,     &
     758             :                                     n_dic,        n_don,     &
     759             :                                     n_fed,        n_fep,     &
     760             :                                     n_zaero,      first_ice, &
     761             :                                     aicen,        vicen,     & 
     762      199620 :                                     hice_old,     ocean_bio, & 
     763      399240 :                                     flux_bio,     bphin,     &
     764      199620 :                                     iphin,        trcrn,     &  
     765      199620 :                                     iDin,   &
     766      399240 :                                     fswthrul,     grow_alg,  &
     767      199620 :                                     upNOn,        upNHn,     &
     768             :                                     dh_top,       dh_bot,    &
     769      199620 :                                     zfswin,       hbri,      & 
     770             :                                     hbri_old,     darcy_V,   &
     771             : !                                   darcy_V_chl,  bgrid,     &
     772      199620 :                                     bgrid,     &
     773      199620 :                                     i_grid,       ic_grid,   &
     774      199620 :                                     bphi_min,     zbgc_snow, &
     775      199620 :                                     zbgc_atm,  &
     776      199620 :                                     iTin,         dh_direct, &
     777      199620 :                                     Zoo,          meltb,     &
     778             :                                     congel                   )
     779             : 
     780             :       integer (kind=int_kind), intent(in) :: &
     781             :          n_cat,              & ! category number
     782             :          nilyr,              & ! number of ice layers
     783             :          nblyr,              & ! number of bio layers
     784             :          nbtrcr, n_algae,    & ! number of bgc tracers, number of autotrophs
     785             :          n_zaero,            & ! number of aerosols
     786             :          n_doc, n_dic,  n_don, n_fed, n_fep
     787             :                               
     788             :       logical (kind=log_kind), intent(in) :: &
     789             :          first_ice      ! initialized values should be used
     790             : 
     791             :       real (kind=dbl_kind), intent(in) :: &
     792             :          dt         , & ! time step 
     793             :          hbri       , & ! brine height  (m)
     794             :          bphi_min   , & ! surface porosity
     795             :          aicen      , & ! concentration of ice
     796             :          vicen      , & ! volume per unit area of ice  (m)
     797             :          hice_old   , & ! ice height (m)
     798             :          meltb      , & ! bottom melt in dt (m)
     799             :          congel     , & ! bottom growth in dt (m)
     800             :          darcy_V    , & ! darcy velocity
     801             : !        darcy_V_chl, & ! darcy velocity for algae
     802             :          dh_bot     , & ! change in brine bottom (m)
     803             :          dh_top     , & ! change in brine top (m)
     804             :          dh_direct      ! surface flooding or runoff (m)
     805             : 
     806             :       real (kind=dbl_kind), dimension (:), intent(inout) :: &
     807             :          bgrid      , & ! biology nondimensional vertical grid points
     808             :          flux_bio   , & ! total ocean tracer flux (mmol/m^2/s)
     809             :          zfswin     , & ! visible Short wave flux on igrid (W/m^2)  
     810             :          Zoo        , & ! N losses to the system from reaction terms
     811             :                         ! (ie. zooplankton/bacteria) (mmol/m^3)  
     812             :          trcrn          ! bulk tracer concentration (mmol/m^3)
     813             : 
     814             :       real (kind=dbl_kind), dimension (:), intent(in) :: &
     815             :          i_grid     , & ! biology vertical interface points
     816             :          iTin       , & ! salinity vertical interface points
     817             :          iphin      , & ! Porosity on the igrid   
     818             :          iDin       , & ! Diffusivity/h on the igrid (1/s)
     819             :          ic_grid    , & ! CICE interface coordinate 
     820             :          fswthrul   , & ! visible short wave radiation on icgrid (W/m^2)  
     821             :          zbgc_snow  , & ! tracer input from snow (mmol/m^3*m)
     822             :          zbgc_atm   , & ! tracer input from  atm (mmol/m^3 *m)
     823             :          ocean_bio  , & ! ocean concentrations (mmol/m^3) 
     824             :          bphin          ! Porosity on the bgrid
     825             : 
     826             :       real (kind=dbl_kind), intent(inout) :: &
     827             :          hbri_old       ! brine height  (m)
     828             : 
     829             :       real (kind=dbl_kind), dimension (:,:), intent(out) :: &
     830             :          upNOn      , & ! algal nitrate uptake rate  (mmol/m^3/s)
     831             :          upNHn      , & ! algal ammonium uptake rate (mmol/m^3/s)
     832             :          grow_alg       ! algal growth rate          (mmol/m^3/s)
     833             : 
     834             :       !-----------------------------------------------------------------------------
     835             :       ! algae absorption coefficient for 0.5 m thick layer
     836             :       ! Grenfell (1991): SA = specific absorption coefficient= 0.004 m^2/mg Chla
     837             :       ! generalizing kalg_bio(k) = SA*\sum R_chl2N(m)*trcrn(i,j,nt_bgc_N(m)+k-1)
     838             :       ! output kalg on the icgrid
     839             :       !-----------------------------------------------------------------------------
     840             : 
     841             :       ! local variables
     842             : 
     843             :       integer (kind=int_kind) :: &
     844             :          k, m, mm        ! vertical biology layer index 
     845             : 
     846             :       real (kind=dbl_kind) :: &
     847       48078 :          hin         , & ! ice thickness (m)        
     848       48078 :          hin_old     , & ! ice thickness before current melt/growth (m)
     849       48078 :          ice_conc    , & ! algal concentration in ice above hin > hinS
     850       48078 :          sum_old     , & !
     851       48078 :          sum_new     , & !
     852       48078 :          sum_tot     , & !
     853       48078 :          zspace      , & ! 1/nblyr
     854       48078 :          darcyV      , & !
     855       48078 :          dhtop       , & !
     856       48078 :          dhbot       , & !
     857       48078 :          dhflood         ! >=0 (m) surface flooding from the ocean
     858             : 
     859             :       real (kind=dbl_kind), dimension (nblyr+2) :: &
     860      783864 :          bphin_N         ! porosity for tracer model has minimum 
     861             :                          ! bphin_N >= bphimin
     862             : 
     863             :       real (kind=dbl_kind), dimension (nblyr+1) :: &
     864      735786 :          iphin_N      , & ! tracer porosity on the igrid
     865      735786 :          sbdiagz      , & ! sub-diagonal matrix elements
     866      735786 :          diagz        , & ! diagonal matrix elements
     867      735786 :          spdiagz      , & ! super-diagonal matrix elements
     868      735786 :          rhsz         , & ! rhs of tri-diagonal matrix equation
     869      735786 :          ML_diag      , & ! lumped mass matrix
     870      735786 :          D_spdiag     , & ! artificial diffusion matrix
     871      735786 :          D_sbdiag     , & ! artificial diffusion matrix
     872      735786 :          biomat_low   , & ! Low order solution
     873      735786 :          Nerror           ! Change in N after reactions
     874             : 
     875             :       real (kind=dbl_kind), dimension(nblyr+1,nbtrcr):: &
     876     8707806 :          react            ! biological sources and sinks for equation matrix
     877             : 
     878             :       real (kind=dbl_kind), dimension(nblyr+1,nbtrcr):: &
     879     8707806 :          in_init_cons , & ! Initial bulk concentration*h (mmol/m^2)
     880     8707806 :          biomat_cons  , & ! Matrix output of (mmol/m^2)
     881     8707806 :          biomat_brine     ! brine concentration (mmol/m^3)
     882             : 
     883             :       real (kind=dbl_kind), dimension(nbtrcr):: &
     884     1279678 :          C_top,         & ! bulk top tracer source: h phi C(meltwater) (mmol/m^2)
     885     1279678 :          C_bot,         & ! bulk bottom tracer source: h phi C(ocean) (mmol/m^2)
     886     1279678 :          Source_top,    & ! For cons: (+) top tracer source into ice (mmol/m^2/s)
     887     1279678 :          Source_bot,    & ! For cons: (+) bottom tracer source into ice (mmol/m^2/s)
     888     1279678 :          Sink_bot,      & ! For cons: (+ or -) remaining bottom flux into ice(mmol/m^2/s)
     889     1279678 :          Sink_top,      & ! For cons: (+ or -) remaining bottom flux into ice(mmol/m^2/s)
     890     1279678 :          exp_ret,       & ! exp dt/retention frequency
     891     1279678 :          exp_rel,       & ! exp dt/release frequency
     892     1375834 :          atm_add_cons , & ! zbgc_snow+zbgc_atm (mmol/m^3*m)
     893     1279678 :          dust_Fe      , & ! contribution of dust surface flux to dFe (umol/m*3*m)
     894     1279678 :          source           ! mmol/m^2 surface input from snow/atmosphere
     895             : 
     896             :       real (kind=dbl_kind), dimension (ntrcr+2) :: &
     897    11718710 :          trtmp0       , & ! temporary, remapped tracers
     898    11718710 :          trtmp            ! temporary, remapped tracers
     899             : 
     900             :       logical (kind=log_kind), dimension(nblyr+1) :: &
     901      399240 :          conserve_N
     902             : 
     903             :       real (kind=dbl_kind), dimension(nblyr+1):: &  ! temporary variables for
     904      735786 :          Diff         , & ! diffusivity 
     905      735786 :          initcons     , & ! initial concentration
     906      735786 :          biocons      , & !  new concentration
     907      735786 :          dmobile      , & !
     908      735786 :          initcons_mobile,&!
     909      735786 :          initcons_stationary
     910             :  
     911             :       real (kind=dbl_kind):: &
     912       48078 :          top_conc         ! 1% (min_bgc) of surface concentration 
     913             :                           ! when hin > hbri:  just used in sw calculation
     914             : 
     915             :       real (kind=dbl_kind):: &
     916       48078 :          bio_tmp, &       ! temporary variable
     917       48078 :          exp_min          ! temporary exp var
     918             : 
     919             :       real (kind=dbl_kind):: &
     920       48078 :          Sat_conc   , & ! adsorbing saturation concentration  (mmols/m^3)
     921       48078 :          phi_max    , & ! maximum porosity
     922       48078 :          S_col      , & ! surface area of collector (um^2)
     923       48078 :          P_b        , & ! projected area of diatoms (um^2)
     924       48078 :          V_c        , & ! volume of collector  (um^3)
     925       48078 :          V_alg          ! volume of algae (um^3)
     926             : 
     927             :       real (kind=dbl_kind), dimension(nbtrcr) :: & 
     928     1224292 :          mobile           ! c1 if mobile, c0 otherwise
     929             : 
     930             :       ! local parameters
     931             :          
     932             :       real (kind=dbl_kind), parameter :: &
     933             :          accuracy = 1.0e-14_dbl_kind, &
     934             :          r_c  = 3.0e3_dbl_kind     , & ! ice crystal radius (um)
     935             :          r_bac= 15.0_dbl_kind    , & ! diatom large radius (um)
     936             :          r_alg= 10.0_dbl_kind    , & ! diatom small radius (um)
     937             :          N_vol = 0.04e-12_dbl_kind  , & ! (g) Nitrogen per um^3
     938             :          Ng_to_mmol =0.0140067_dbl_kind , & ! (g/mmol) Nitrogen
     939             :          f_s = c1 , &  ! fracton of sites available for saturation
     940             :          f_a = c1 , &  ! fraction of collector available for attachment
     941             :          f_v = 0.7854  ! fraction of algal coverage on area availabel for attachment
     942             :                        ! 4(pi r^2)/(4r)^2  [Johnson et al, 1995, water res. research]
     943             :           
     944             :       integer, parameter :: &
     945             :          nt_zfswin = 1    ! for interpolation of short wave to bgrid
     946             : 
     947             :       character(len=*),parameter :: subname='(z_biogeochemistry)'
     948             : 
     949             :   !-------------------------------------
     950             :   ! Initialize 
     951             :   !----------------------------------- 
     952             : 
     953      199620 :       zspace = c1/real(nblyr,kind=dbl_kind)
     954    35833824 :       in_init_cons(:,:) = c0
     955     4158976 :       atm_add_cons(:) = c0
     956      199620 :       dhtop = c0
     957      199620 :       dhbot = c0
     958      199620 :       darcyV = c0
     959     4158976 :       C_top(:) = c0
     960     4158976 :       mobile(:) = c0
     961     1796580 :       conserve_N(:) = .true.
     962             : 
     963     4158976 :       do m = 1, nbtrcr
     964    35833824 :          do k  = 1, nblyr+1
     965             : 
     966    31674848 :             bphin_N(nblyr+2) =c1
     967    31674848 :             bphin_N(k) = bphin(k)
     968    31674848 :             iphin_N(k) = iphin(k)
     969    31674848 :             bphin_N(1) = bphi_min
     970             : 
     971    31674848 :             if (first_ice) then
     972      599200 :                trcrn(bio_index(m) + k-1) = ocean_bio(m)*zbgc_init_frac(m)
     973      599200 :                in_init_cons(k,m) = trcrn(bio_index(m) + k-1)*hbri_old
     974    31075648 :             elseif (abs(trcrn(bio_index(m) + k-1)) < puny) then               
     975     5053763 :                trcrn(bio_index(m) + k-1) = c0
     976     5053763 :                in_init_cons(k,m) = c0
     977             :             else
     978    26021885 :                in_init_cons(k,m) = trcrn(bio_index(m) + k-1)* hbri_old
     979             :             endif         ! first_ice
     980             : 
     981    31674848 :             if (trcrn(bio_index(m) + k-1) < c0  ) then
     982           0 :                write(warnstr,*) subname,'zbgc initialization error, first ice = ', first_ice
     983           0 :                call icepack_warnings_add(warnstr)
     984           0 :                write(warnstr,*) subname,'Category,m:',n_cat,m
     985           0 :                call icepack_warnings_add(warnstr)
     986           0 :                write(warnstr,*) subname,'hbri,hbri_old' 
     987           0 :                call icepack_warnings_add(warnstr)
     988           0 :                write(warnstr,*) subname, hbri,hbri_old
     989           0 :                call icepack_warnings_add(warnstr)
     990           0 :                write(warnstr,*) subname,'trcrn(bio_index(m) + k-1)'
     991           0 :                call icepack_warnings_add(warnstr)
     992           0 :                write(warnstr,*) subname, trcrn(bio_index(m) + k-1)
     993           0 :                call icepack_warnings_add(warnstr)
     994           0 :                call icepack_warnings_add(subname//' zbgc initialization error')
     995           0 :                call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
     996             :             endif 
     997    35634204 :             if (icepack_warnings_aborted(subname)) return
     998             :         enddo         !k
     999             :       enddo           !m
    1000             : 
    1001             :       !-----------------------------------------------------------------
    1002             :       !     boundary conditions
    1003             :       !-----------------------------------------------------------------
    1004             : 
    1005      199620 :       ice_conc = c0
    1006      199620 :       hin = vicen/aicen
    1007      199620 :       hin_old = hice_old       
    1008             : 
    1009             :       !-----------------------------------------------------------------
    1010             :       !    calculate the saturation concentration for attachment: Sat_conc
    1011             :       !-----------------------------------------------------------------
    1012             : 
    1013     1796580 :       phi_max = maxval(bphin_N(2:nblyr+1))
    1014      199620 :       S_col   = 4.0_dbl_kind*pi*r_c**2
    1015      199620 :       P_b     = pi*r_bac**2    !*10-6 for colloids
    1016      199620 :       V_c     = 4.0_dbl_kind*pi*r_c**3/3.0_dbl_kind*(1.0e-6_dbl_kind)**3  ! (m^3) sphere
    1017      199620 :       V_alg   = pi/6.0_dbl_kind*r_bac*r_alg**2       ! prolate spheroid (*10-9 for colloids)
    1018      199620 :       Sat_conc= f_s*f_a*f_v*(c1-phi_max)/V_c*S_col/P_b*N_vol*V_alg/Ng_to_mmol
    1019             :                             !mmol/m^3 (algae, don, hum...) and umols/m^3 for colloids 
    1020             : 
    1021             :       !-----------------------------------------------------------------
    1022             :       !    convert surface dust flux (n_zaero > 2) to dFe(1) flux    
    1023             :       !-----------------------------------------------------------------
    1024             : 
    1025     4158976 :       dust_Fe(:) = c0
    1026             : 
    1027      199620 :       if (tr_zaero .and. n_zaero > 2 .and. tr_bgc_Fe) then
    1028           0 :        do m = 3,n_zaero
    1029           0 :          dust_Fe(nlt_bgc_Fed(1)) = dust_Fe(nlt_bgc_Fed(1)) + &
    1030           0 :               (zbgc_snow(nlt_zaero(m)) + zbgc_atm(nlt_zaero(m))) * &
    1031           0 :                R_dFe2dust * dustFe_sol
    1032             :         ! dust_Fe(nlt_zaero(m)) = -(zbgc_snow(nlt_zaero(m)) + zbgc_atm(nlt_zaero(m))) * &
    1033             :         !       dustFe_sol
    1034             :        enddo  
    1035             :       endif
    1036             : 
    1037     4158976 :       do m = 1,nbtrcr 
    1038             :       !-----------------------------------------------------------------
    1039             :       !   time constants for mobile/stationary phase changes
    1040             :       !-----------------------------------------------------------------
    1041             : 
    1042     3959356 :          exp_rel(m) = c0
    1043     3959356 :          exp_ret(m) = c0
    1044     3959356 :          if (tau_ret(m) > c0) then
    1045     3959356 :             exp_min = min(dt/tau_ret(m),exp_argmax)
    1046     3959356 :             exp_ret(m) = exp(-exp_min)
    1047             :          endif
    1048     3959356 :          if (tau_rel(m) > c0) then
    1049     3959356 :             exp_min = min(dt/tau_rel(m),exp_argmax)
    1050     3959356 :             exp_rel(m) = exp(-exp_min)
    1051             :          endif
    1052     3959356 :          if (m .ne. nlt_bgc_N(1)) then  
    1053     3759736 :             if (hin_old  > hin) then  !melting
    1054     1874644 :                exp_ret(m) = c1
    1055             :             else                              !not melting
    1056     1885092 :                exp_rel(m) = c1
    1057             :             endif  
    1058      199620 :          elseif (tr_bgc_N .and. hin_old > hin + algal_vel*dt) then
    1059       38598 :                exp_ret(m) = c1
    1060      161022 :          elseif (tr_bgc_N) then
    1061      161022 :                exp_rel(m) = c1
    1062             :          endif
    1063             : 
    1064     3959356 :          dhtop      = dh_top
    1065     3959356 :          dhbot      = dh_bot
    1066     3959356 :          darcyV     = darcy_V
    1067     3959356 :          C_top(m)   = in_init_cons(1,m)*trcrn(nt_zbgc_frac+m-1)!mobile fraction
    1068     3959356 :          source(m)  = abs(zbgc_snow(m) + zbgc_atm(m) + dust_Fe(m))
    1069     3959356 :          dhflood  = max(c0,-dh_direct)                              ! ocean water flooding surface
    1070             : 
    1071     3959356 :          if (dhtop+darcyV/bphin_N(1)*dt < -puny) then !snow/top ice melt
    1072      635340 :              C_top(m) = (zbgc_snow(m)+zbgc_atm(m) + dust_Fe(m))/abs(dhtop &
    1073     3083060 :                         + darcyV/bphin_N(1)*dt + puny)*hbri_old    
    1074     1169472 :          elseif (dhtop+darcyV/bphin_N(1)*dt >= -puny .and. &
    1075      293176 :                         abs((zbgc_snow(m)+zbgc_atm(m) + dust_Fe(m)) + &
    1076      293176 :                         ocean_bio(m)*bphin_N(1)*dhflood) >  puny) then
    1077        1674 :               atm_add_cons(m) =  abs(zbgc_snow(m) + zbgc_atm(m)+ dust_Fe(m)) + &
    1078        5694 :                                       ocean_bio(m)*bphin_N(1)*dhflood      
    1079             :          else   ! only positive fluxes 
    1080      870602 :               atm_add_cons(m) =  abs(zbgc_snow(m) + zbgc_atm(m)+ dust_Fe(m))
    1081             :          endif
    1082             : 
    1083     4158976 :          C_bot(m) = ocean_bio(m)*hbri_old*iphin_N(nblyr+1)            
    1084             : 
    1085             :       enddo             ! m
    1086             : 
    1087             :       !-----------------------------------------------------------------
    1088             :       ! Interpolate shortwave flux, fswthrul (defined at top to bottom with nilyr+1 
    1089             :       !  evenly spaced  with spacing = (1/nilyr) to grid variable zfswin:
    1090             :       !-----------------------------------------------------------------
    1091             : 
    1092    48543416 :       trtmp(:) = c0 
    1093    48543416 :       trtmp0(:)= c0
    1094     1796580 :       zfswin(:) = c0
    1095             : 
    1096     1796580 :       do k = 1, nilyr+1
    1097             :          ! contains cice values (fswthrul(1) is surface value)
    1098             :          ! and fwsthrul(nilyr+1) is output
    1099     1796580 :          trtmp0(nt_zfswin+k-1) = fswthrul(k) 
    1100             :       enddo   !k
    1101             : 
    1102             :       call remap_zbgc(nilyr+1,  &
    1103             :                       nt_zfswin,                  &
    1104           0 :                       trtmp0(1:ntrcr),  trtmp(1:ntrcr+2), &
    1105             :                       0,                nblyr+1,  &
    1106             :                       hin,              hbri,     &
    1107           0 :                       ic_grid(1:nilyr+1),         &
    1108      199620 :                       i_grid(1:nblyr+1),ice_conc  )
    1109             : 
    1110      199620 :       if (icepack_warnings_aborted(subname)) return
    1111             : 
    1112     1796580 :       do k = 1,nblyr+1
    1113     1796580 :          zfswin(k) = trtmp(nt_zfswin+k-1)
    1114             :       enddo
    1115             : 
    1116             :       !-----------------------------------------------------------------
    1117             :       ! Initialze Biology  
    1118             :       !----------------------------------------------------------------- 
    1119             : 
    1120     4158976 :       do mm = 1, nbtrcr
    1121     3959356 :          mobile(mm) = c0
    1122     3959356 :          if (bgc_tracer_type(mm) .GE. c0) mobile(mm) = c1
    1123             : 
    1124    35833824 :          do k = 1, nblyr+1
    1125    35634204 :             biomat_cons(k,mm) = in_init_cons(k,mm)
    1126             :          enddo  !k
    1127             :       enddo  !mm
    1128             : 
    1129             :       !-----------------------------------------------------------------
    1130             :       ! Compute FCT
    1131             :       !----------------------------------------------------------------- 
    1132             : 
    1133     4158976 :       do mm = 1, nbtrcr 
    1134             : 
    1135     3959356 :          if (hbri_old > thinS .and. hbri > thinS) then 
    1136    35634060 :             do k = 1,nblyr+1
    1137    31674720 :                initcons_mobile(k) = in_init_cons(k,mm)*trcrn(nt_zbgc_frac+mm-1)
    1138    31674720 :                initcons_stationary(k) = mobile(mm)*(in_init_cons(k,mm)-initcons_mobile(k))
    1139     7428000 :                dmobile(k) = mobile(mm)*(initcons_mobile(k)*(exp_ret(mm)-c1) + &
    1140    31674720 :                                     initcons_stationary(k)*(c1-exp_rel(mm)))
    1141    31674720 :                initcons_mobile(k) = max(c0,initcons_mobile(k) + dmobile(k))
    1142    31674720 :                initcons_stationary(k) = max(c0,initcons_stationary(k) - dmobile(k))
    1143    31674720 :                if (initcons_stationary(k)/hbri_old > Sat_conc) then
    1144           0 :                   initcons_mobile(k) = initcons_mobile(k) + initcons_stationary(k) - Sat_conc*hbri_old
    1145           0 :                    initcons_stationary(k) = Sat_conc*hbri_old
    1146             :                endif
    1147             : 
    1148    31674720 :                Diff(k) = iDin(k) 
    1149    31674720 :                initcons(k) = initcons_mobile(k)                          
    1150    35634060 :                biocons(k) =  initcons_mobile(k)
    1151             :             enddo
    1152             : 
    1153             :             call compute_FCT_matrix &
    1154             :                                 (initcons,sbdiagz, dt, nblyr,  &
    1155           0 :                                 diagz, spdiagz, rhsz, bgrid,   & 
    1156             :                                 darcyV,    dhtop,      &
    1157             :                                 dhbot,   iphin_N,              &
    1158             :                                 Diff, hbri_old,                &
    1159      928500 :                                 atm_add_cons(mm), bphin_N,     &
    1160      928500 :                                 C_top(mm), C_bot(mm),          &
    1161      928500 :                                 Source_bot(mm), Source_top(mm),&
    1162      928500 :                                 Sink_bot(mm),Sink_top(mm),     &
    1163     4887840 :                                 D_sbdiag, D_spdiag, ML_diag)
    1164     3959340 :             if (icepack_warnings_aborted(subname)) return
    1165             : 
    1166             :             call tridiag_solverz &
    1167             :                                (nblyr+1, sbdiagz,               &
    1168             :                                 diagz,   spdiagz,               &
    1169     3959340 :                                 rhsz,    biocons)
    1170     3959340 :             if (icepack_warnings_aborted(subname)) return
    1171             : 
    1172             :             call check_conservation_FCT &
    1173             :                                (initcons,    &
    1174             :                                 biocons,     &
    1175             :                                 biomat_low,               &
    1176      928500 :                                 Source_top(mm),        &
    1177      928500 :                                 Source_bot(mm),        &
    1178      928500 :                                 Sink_bot(mm),          &
    1179      928500 :                                 Sink_top(mm),          &
    1180      928500 :                                 dt, flux_bio(mm),     &
    1181             :                                 nblyr, &
    1182     3959340 :                                 source(mm))
    1183     3959340 :             if (icepack_warnings_aborted(subname)) return
    1184             : 
    1185             :             call compute_FCT_corr & 
    1186             :                                 (initcons,   &
    1187             :                                  biocons, dt, nblyr, &
    1188     3959340 :                                  D_sbdiag, D_spdiag, ML_diag)  
    1189     3959340 :             if (icepack_warnings_aborted(subname)) return
    1190             : 
    1191     3959340 :             top_conc = c0        ! or frazil ice concentration
    1192             :  
    1193             :             ! assume diatoms actively maintain there relative position in the ice
    1194             : 
    1195     3959340 :             if (mm .ne. nlt_bgc_N(1)) then    
    1196             :                        
    1197             :                call regrid_stationary & 
    1198             :                                 (initcons_stationary,    hbri_old,    &
    1199             :                                  hbri,                   dt,          &
    1200             :                                  ntrcr,                               &
    1201             :                                  nblyr,                  top_conc,    &
    1202      880423 :                                  i_grid,                 flux_bio(mm),&
    1203     4640144 :                                  meltb,                  congel)
    1204     3759721 :                if (icepack_warnings_aborted(subname)) return
    1205             : 
    1206      199619 :             elseif (tr_bgc_N .and. mm .eq. nlt_bgc_N(1)) then  
    1207      199619 :                if (meltb > algal_vel*dt .or. aicen < 0.001_dbl_kind) then
    1208             : 
    1209             :                   call regrid_stationary &
    1210             :                                 (initcons_stationary,    hbri_old,    &
    1211             :                                  hbri,                   dt,          &
    1212             :                                  ntrcr,                               &
    1213             :                                  nblyr,                  top_conc,    &
    1214       11023 :                                  i_grid,                 flux_bio(mm),&
    1215       57942 :                                  meltb,                  congel)       
    1216       46919 :                   if (icepack_warnings_aborted(subname)) return
    1217             : 
    1218             :                endif
    1219             :             endif
    1220             : 
    1221    35634060 :             biomat_cons(:,mm) =  biocons(:) +  initcons_stationary(:)
    1222             : 
    1223     3959340 :             sum_old = (biomat_low(1) + biomat_low(nblyr+1))*zspace/c2
    1224     3959340 :             sum_new = (biocons(1)+ biocons(nblyr+1))*zspace/c2
    1225     3959340 :             sum_tot = (biomat_cons(1,mm) + biomat_cons(nblyr+1,mm))*zspace/c2
    1226    27715380 :             do k = 2,nblyr
    1227    23756040 :                sum_old = sum_old + biomat_low(k)*zspace
    1228    23756040 :                sum_new = sum_new + biocons(k)*zspace
    1229    27715380 :                sum_tot = sum_tot + biomat_cons(k,mm)*zspace
    1230             :             enddo
    1231     3959340 :             trcrn(nt_zbgc_frac+mm-1) = zbgc_frac_init(mm)
    1232     3959340 :             if (sum_tot > c0 .and. mobile(mm) > c0) trcrn(nt_zbgc_frac+mm-1) = sum_new/sum_tot
    1233             : 
    1234             :             if (abs(sum_new-sum_old) > accuracy*sum_old .or. &
    1235           0 :                 minval(biocons(:)) < c0  .or. minval(initcons_stationary(:)) < c0 &
    1236    75227460 :                 .or. icepack_warnings_aborted()) then
    1237           0 :                 write(warnstr,*) subname,'zbgc FCT tracer solution failed'
    1238           0 :                 call icepack_warnings_add(warnstr)
    1239           0 :                 write(warnstr,*) subname,'sum_new,sum_old:',sum_new,sum_old
    1240           0 :                 call icepack_warnings_add(warnstr)
    1241           0 :                 write(warnstr,*) subname,'mm,biocons(:):',mm,biocons(:)
    1242           0 :                 call icepack_warnings_add(warnstr)
    1243           0 :                 write(warnstr,*) subname,'biomat_low:',biomat_low
    1244           0 :                 call icepack_warnings_add(warnstr)
    1245           0 :                 write(warnstr,*) subname,'Diff(:):',Diff(:)
    1246           0 :                 call icepack_warnings_add(warnstr)
    1247           0 :                 write(warnstr,*) subname,'dmobile(:):',dmobile(:)
    1248           0 :                 call icepack_warnings_add(warnstr)
    1249           0 :                 write(warnstr,*) subname,'mobile(mm):',mobile(mm)
    1250           0 :                 call icepack_warnings_add(warnstr)
    1251           0 :                 write(warnstr,*) subname,'initcons_stationary(:):',initcons_stationary(:)
    1252           0 :                 call icepack_warnings_add(warnstr)
    1253           0 :                 write(warnstr,*) subname, 'trcrn(nt_zbgc_frac+mm-1):',trcrn(nt_zbgc_frac+mm-1)
    1254           0 :                 call icepack_warnings_add(warnstr)
    1255           0 :                 write(warnstr,*) subname, 'in_init_cons(:,mm):',in_init_cons(:,mm)
    1256           0 :                 call icepack_warnings_add(warnstr)
    1257           0 :                 write(warnstr,*) subname, 'exp_ret( mm),exp_rel( mm)',exp_ret( mm),exp_rel( mm)
    1258           0 :                 call icepack_warnings_add(warnstr)
    1259           0 :                 write(warnstr,*) subname,'darcyV,dhtop,dhbot'
    1260           0 :                 call icepack_warnings_add(warnstr)
    1261           0 :                 write(warnstr,*) subname,darcyV,dhtop,dhbot
    1262           0 :                 call icepack_warnings_add(warnstr)
    1263           0 :                 write(warnstr,*) subname,'Category,mm:',n_cat,mm
    1264           0 :                 call icepack_warnings_add(warnstr)
    1265           0 :                 call icepack_warnings_add(subname//'zbgc FCT tracer solution failed')
    1266           0 :                 call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
    1267             :             endif
    1268     3959340 :             if (icepack_warnings_aborted(subname)) return
    1269             : 
    1270             :          else              
    1271             :   
    1272           0 :             call thin_ice_flux(hbri,hbri_old,biomat_cons(:,mm), &
    1273          16 :                                flux_bio(mm),source(mm), &
    1274          32 :                                dt, nblyr,ocean_bio(mm))
    1275          16 :             if (icepack_warnings_aborted(subname)) return
    1276             : 
    1277             :          endif ! thin or not
    1278             : 
    1279    35833824 :          do k = 1,nblyr+1 
    1280    35634204 :             biomat_brine(k,mm) =  biomat_cons(k,mm)/hbri/iphin_N(k) 
    1281             :          enddo ! k
    1282             :       enddo ! mm 
    1283             : 
    1284    35833824 :       react(:,:) = c0  
    1285     5589360 :       grow_alg(:,:) = c0
    1286             : 
    1287      199620 :       if (solve_zbgc) then
    1288     1796580 :          do k = 1, nblyr+1   
    1289             :             call algal_dyn (dt,              &
    1290             :                          n_doc, n_dic,  n_don, n_fed, n_fep, &
    1291             :                          dEdd_algae, &
    1292      384624 :                          zfswin(k),        react(k,:),     & 
    1293           0 :                          biomat_brine(k,:), &
    1294           0 :                          grow_alg(k,:),    n_algae,        &
    1295      384624 :                          iTin(k),                          &
    1296           0 :                          upNOn(k,:),       upNHn(k,:),     &
    1297      384624 :                          Zoo(k),                           &
    1298     2366208 :                          Nerror(k),        conserve_N(k))
    1299     1796580 :             if (icepack_warnings_aborted(subname)) return
    1300             :                          
    1301             :          enddo ! k
    1302             :       endif    ! solve_zbgc
    1303             : 
    1304             :       !-----------------------------------------------------------------
    1305             :       ! Update the tracer variable
    1306             :       !-----------------------------------------------------------------
    1307             :     
    1308     4158976 :       do m = 1,nbtrcr
    1309    35833824 :          do k = 1,nblyr+1                  ! back to bulk quantity
    1310    31674848 :             bio_tmp = (biomat_brine(k,m) + react(k,m))*iphin_N(k) 
    1311             :                      
    1312    31674848 :             if (.not. conserve_N(k)) then  
    1313           0 :                 write(warnstr,*) subname, 'N in algal_dyn not conserved'
    1314           0 :                 call icepack_warnings_add(warnstr)
    1315           0 :                 write(warnstr,*) subname, 'Nerror(k):', Nerror(k)
    1316           0 :                 call icepack_warnings_add(warnstr)
    1317           0 :                 write(warnstr,*) subname, 'k,m,hbri,hbri_old,bio_tmp,biomat_cons(k,m),ocean_bio(m)'
    1318           0 :                 call icepack_warnings_add(warnstr)
    1319           0 :                 write(warnstr,*) subname,  k,m,hbri,hbri_old,bio_tmp,biomat_cons(k,m),ocean_bio(m)
    1320           0 :                 call icepack_warnings_add(warnstr)
    1321           0 :                 write(warnstr,*) subname, 'react(k,m),iphin_N(k),biomat_brine(k,m)'
    1322           0 :                 call icepack_warnings_add(warnstr)
    1323           0 :                 write(warnstr,*) subname,  react(k,m),iphin_N(k),biomat_brine(k,m)
    1324           0 :                 call icepack_warnings_add(warnstr)
    1325           0 :                 call icepack_warnings_add(subname//' N in algal_dyn not conserved')
    1326           0 :                 call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
    1327    31674848 :             elseif (abs(bio_tmp) < puny) then  
    1328     5211208 :                bio_tmp = c0
    1329    26463640 :             elseif (bio_tmp > 1.0e6_dbl_kind) then
    1330           0 :                 write(warnstr,*) subname, 'very large bgc value'
    1331           0 :                 call icepack_warnings_add(warnstr)
    1332           0 :                 write(warnstr,*) subname, 'k,m,hbri,hbri_old,bio_tmp,biomat_cons(k,m),ocean_bio(m)'
    1333           0 :                 call icepack_warnings_add(warnstr)
    1334           0 :                 write(warnstr,*) subname,  k,m,hbri,hbri_old,bio_tmp,biomat_cons(k,m),ocean_bio(m)
    1335           0 :                 call icepack_warnings_add(warnstr)
    1336           0 :                 write(warnstr,*) subname, 'react(k,m),iphin_N(k),biomat_brine(k,m)'
    1337           0 :                 call icepack_warnings_add(warnstr)
    1338           0 :                 write(warnstr,*) subname,  react(k,m),iphin_N(k),biomat_brine(k,m)
    1339           0 :                 call icepack_warnings_add(warnstr)
    1340           0 :                 call icepack_warnings_add(subname//' very large bgc value')
    1341           0 :                 call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
    1342    26463640 :             elseif (bio_tmp < c0) then
    1343           0 :                 write(warnstr,*) subname, 'negative bgc'
    1344           0 :                 call icepack_warnings_add(warnstr)
    1345           0 :                 write(warnstr,*) subname, 'k,m,nlt_bgc_Nit,hbri,hbri_old'
    1346           0 :                 call icepack_warnings_add(warnstr)
    1347           0 :                 write(warnstr,*) subname,  k,m,nlt_bgc_Nit,hbri,hbri_old
    1348           0 :                 call icepack_warnings_add(warnstr)
    1349           0 :                 write(warnstr,*) subname, 'bio_tmp,biomat_cons(k,m),ocean_bio(m)'
    1350           0 :                 call icepack_warnings_add(warnstr)
    1351           0 :                 write(warnstr,*) subname,  bio_tmp,biomat_cons(k,m),ocean_bio(m)
    1352           0 :                 call icepack_warnings_add(warnstr)
    1353           0 :                 write(warnstr,*) subname, 'react(k,m),iphin_N(k),biomat_brine(k,m)'
    1354           0 :                 call icepack_warnings_add(warnstr)
    1355           0 :                 write(warnstr,*) subname,  react(k,m),iphin_N(k),biomat_brine(k,m)
    1356           0 :                 call icepack_warnings_add(warnstr)
    1357           0 :                 write(warnstr,*) subname, 'exp_ret( m),exp_ret( m)',exp_ret( m),exp_ret( m)
    1358           0 :                 call icepack_warnings_add(warnstr)
    1359           0 :                 call icepack_warnings_add(subname//'negative bgc')
    1360           0 :                 call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
    1361             :             endif
    1362    31674848 :             if (icepack_warnings_aborted()) then
    1363           0 :                 write(warnstr,*) subname, 'trcrn(nt_zbgc_frac+m-1):',trcrn(nt_zbgc_frac+m-1)
    1364           0 :                 call icepack_warnings_add(warnstr)
    1365           0 :                 write(warnstr,*) subname, 'in_init_cons(k,m):',in_init_cons(k,m)
    1366           0 :                 call icepack_warnings_add(warnstr)
    1367           0 :                 write(warnstr,*) subname, 'trcrn(bio_index(m) + k-1)'
    1368           0 :                 call icepack_warnings_add(warnstr)
    1369           0 :                 write(warnstr,*) subname,  trcrn(bio_index(m) + k-1)
    1370           0 :                 call icepack_warnings_add(warnstr)
    1371           0 :                 write(warnstr,*) subname, 'Category,m:',n_cat,m
    1372           0 :                 call icepack_warnings_add(warnstr)
    1373           0 :                 return
    1374             :             endif
    1375    31674848 :             trcrn(bio_index(m)+k-1) = max(c0, bio_tmp)
    1376    35634204 :             if (ocean_bio(m) .le. c0 .and. flux_bio(m) < c0) then
    1377             :            !     if (flux_bio(m) < -1.0e-12_dbl_kind) then
    1378             :            !       write(warnstr,*) subname, 'no ocean_bio but flux_bio < c0'
    1379             :            !       call icepack_warnings_add(warnstr)
    1380             :            !       write(warnstr,*) subname, 'm,ocean_bio(m),flux_bio(m)'
    1381             :            !       call icepack_warnings_add(warnstr)
    1382             :            !       write(warnstr,*) subname, m,ocean_bio(m),flux_bio(m)
    1383             :            !       call icepack_warnings_add(warnstr)
    1384             :            !       write(warnstr,*) subname, 'setting flux_bio(m) = c0'
    1385             :            !       call icepack_warnings_add(warnstr)
    1386             :            !       call icepack_warnings_add(subname//' flux_bio < 0 when ocean_bio = 0')
    1387             :            !       call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
    1388             :            !     endif
    1389      644706 :                 flux_bio(m) = max(c0,flux_bio(m))
    1390             :             endif
    1391             :          enddo        ! k
    1392             :       enddo        ! m   
    1393             :    
    1394             : 770 format (I6,D16.6)        
    1395             : 781 format (I6,I6,I6)
    1396             : 790 format (I6,I6)
    1397             : 791 format (f24.17)
    1398             : 792 format (2D16.6)
    1399             : 793 format (3D16.6)
    1400             : 794 format (4D15.5)
    1401             : 800 format (F10.4)
    1402             : 
    1403             :       end subroutine z_biogeochemistry
    1404             : 
    1405             : !=======================================================================
    1406             : !
    1407             : ! Do biogeochemistry from subroutine algal_dynamics
    1408             : ! authors: Scott Elliott, LANL
    1409             : !          Nicole Jeffery, LANL
    1410             : 
    1411     1605221 :       subroutine algal_dyn (dt,           &
    1412             :                             n_doc, n_dic,  n_don, n_fed, n_fep, &
    1413             :                             dEdd_algae,   &
    1414     1605221 :                             fswthru,      reactb,       & 
    1415     1605221 :                             ltrcrn,       &
    1416     1605221 :                             grow_alg,     n_algae,      &
    1417             :                             T_bot,                      &
    1418     1605221 :                             upNOn,        upNHn,        &
    1419             :                             Zoo,                        &
    1420             :                             Nerror,       conserve_N)      
    1421             : 
    1422             :       integer (kind=int_kind), intent(in) :: &
    1423             :          n_doc, n_dic,  n_don, n_fed, n_fep, &
    1424             :          n_algae    ! number of autotrophic types
    1425             : 
    1426             :       real (kind=dbl_kind), intent(in) :: &
    1427             :          dt      , & ! time step 
    1428             :          T_bot   , & ! ice temperature (oC)
    1429             :          fswthru     ! average shortwave passing through current ice layer (W/m^2)
    1430             : 
    1431             :       real (kind=dbl_kind), intent(inout) :: &
    1432             :          Zoo,     & ! N losses from zooplankton/bacteria... (mmol/m^3)
    1433             :          Nerror     ! Change in N after reactions (mmol/m^3)
    1434             : 
    1435             :       real (kind=dbl_kind), dimension (:), intent(out) :: &
    1436             :          grow_alg,& !  algal growth rate   (mmol/m^3/s)
    1437             :          upNOn,   & !  algal NO uptake rate   (mmol/m^3/s)
    1438             :          upNHn      !  algal NH uptake rate   (mmol/m^3/s)
    1439             : 
    1440             :       real (kind=dbl_kind), dimension(:), intent(inout) :: &
    1441             :          reactb     ! biological reaction terms (mmol/m3)
    1442             : 
    1443             :       real (kind=dbl_kind), dimension(:), intent(in) :: &
    1444             :          ltrcrn     ! brine concentrations  in layer (mmol/m^3) 
    1445             : 
    1446             :       logical (kind=log_kind), intent(inout) :: & 
    1447             :          conserve_N
    1448             : 
    1449             :       logical (kind=log_kind), intent(in) :: & 
    1450             :          dEdd_algae  ! .true.  chla impact on shortwave computed in dEdd
    1451             : 
    1452             :       !  local variables
    1453             :       !------------------------------------------------------------------------------------
    1454             :       !            3 possible autotrophs nt_bgc_N(1:3):  diatoms, flagellates, phaeocystis
    1455             :       !                2 types of dissolved organic carbon nt_bgc_DOC(1:2): 
    1456             :       !                        polysaccharids, lipids
    1457             :       !                1 DON (proteins)
    1458             :       !                1 particulate iron (nt_bgc_Fe) n_fep
    1459             :       !                1 dossp;ved orpm m+fed 
    1460             :       ! Limiting macro/micro nutrients: nt_bgc_Nit -> nitrate, nt_bgc_NH -> ammonium, 
    1461             :       !                        nt_bgc_Sil -> silicate, nt_bgc_Fe -> dissolved iron   
    1462             :       ! --------------------------------------------------------------------------------------
    1463             : 
    1464             : !     real (kind=dbl_kind),  parameter, dimension(max_algae) :: &
    1465             : !        alpha2max_high  = (/ 0.25_dbl_kind, 0.25_dbl_kind, 0.25_dbl_kind/) ! light limitation (1/(W/m^2))
    1466             : 
    1467             :       integer (kind=int_kind) :: k, n
    1468             : 
    1469             :       real (kind=dbl_kind), dimension(n_algae) :: &
    1470     3996212 :          Nin        , &     ! algal nitrogen concentration on volume (mmol/m^3) 
    1471             : !        Cin        , &     ! algal carbon concentration on volume (mmol/m^3)
    1472     3996212 :          chlin              ! algal chlorophyll concentration on volume (mg/m^3)
    1473             : 
    1474             :       real (kind=dbl_kind), dimension(n_doc) :: &
    1475     3603327 :          DOCin              ! dissolved organic carbon concentration on volume (mmolC/m^3) 
    1476             : 
    1477             : !     real (kind=dbl_kind), dimension(n_dic) :: &
    1478             : !        DICin              ! dissolved inorganic carbon concentration on volume (mmolC/m^3) 
    1479             : 
    1480             :       real (kind=dbl_kind), dimension(n_don) :: &  !proteins
    1481     3210442 :          DONin              ! dissolved organic nitrogen concentration on volume (mmolN/m^3) 
    1482             : 
    1483             :       real (kind=dbl_kind), dimension(n_fed) :: &  !iron
    1484     3210442 :          Fedin              ! dissolved iron concentration on volume (umol/m^3) 
    1485             : 
    1486             :       real (kind=dbl_kind), dimension(n_fep) :: &  !iron
    1487     3210442 :          Fepin              ! algal nitrogen concentration on volume (umol/m^3) 
    1488             : 
    1489             :       real (kind=dbl_kind) :: &
    1490      392885 :          Nitin      , &     ! nitrate concentration on volume (mmol/m^3) 
    1491      392885 :          Amin       , &     ! ammonia/um concentration on volume (mmol/m^3) 
    1492      392885 :          Silin      , &     ! silicon concentration on volume (mmol/m^3) 
    1493             : !        DMSPpin    , &     ! DMSPp concentration on volume (mmol/m^3)
    1494      392885 :          DMSPdin    , &     ! DMSPd concentration on volume (mmol/m^3)
    1495      392885 :          DMSin      , &     ! DMS concentration on volume (mmol/m^3)
    1496             : !        PONin      , &     ! PON concentration on volume (mmol/m^3)
    1497      392885 :          op_dep     , &     ! bottom layer attenuation exponent (optical depth)
    1498      392885 :          Iavg_loc           ! bottom layer attenuated Fswthru (W/m^2)
    1499             : 
    1500             :       real (kind=dbl_kind), dimension(n_algae) :: &
    1501     3996212 :          L_lim    , &  ! overall light limitation
    1502     3996212 :          Nit_lim  , &  ! overall nitrate limitation
    1503     3996212 :          Am_lim   , &  ! overall ammonium limitation
    1504     3996212 :          N_lim    , &  ! overall nitrogen species limitation
    1505     3996212 :          Sil_lim  , &  ! overall silicon limitation
    1506     3996212 :          Fe_lim   , &  ! overall iron limitation
    1507     3996212 :          fr_Nit   , &  ! fraction of local ecological growth as nitrate
    1508     3996212 :          fr_Am    , &  ! fraction of local ecological growth as ammonia
    1509     3996212 :          growmax_N, &  ! maximum growth rate in N currency (mmol/m^3/s)
    1510     3996212 :          grow_N   , &  ! true growth rate in N currency (mmol/m^3/s)
    1511             : !        potU_Nit , &  ! potential nitrate uptake (mmol/m^3/s)
    1512     3996212 :          potU_Am  , &  ! potential ammonium uptake (mmol/m^3/s)
    1513     3996212 :          U_Nit    , &  ! actual nitrate uptake (mmol/m^3/s)
    1514     3996212 :          U_Am     , &  ! actual ammonium uptake (mmol/m^3/s)
    1515     3996212 :          U_Sil    , &  ! actual silicon uptake (mmol/m^3/s)
    1516     3996212 :          U_Fe     , &  ! actual iron uptake   (umol/m^3/s)
    1517     3996212 :          U_Nit_f  , &  ! fraction of Nit uptake due to each algal species
    1518     3996212 :          U_Am_f   , &  ! fraction of Am uptake due to each algal species
    1519     3996212 :          U_Sil_f  , &  ! fraction of Sil uptake due to each algal species
    1520     3996212 :          U_Fe_f        ! fraction of Fe uptake due to each algal species
    1521             : 
    1522             :       real (kind=dbl_kind) :: &
    1523      392885 :          dTemp        , &  ! sea ice temperature minus sst (oC) < 0
    1524      392885 :          U_Nit_tot    , &  ! actual nitrate uptake (mmol/m^3/s)
    1525      392885 :          U_Am_tot     , &  ! actual ammonium uptake (mmol/m^3/s)
    1526      392885 :          U_Sil_tot    , &  ! actual silicon uptake (mmol/m^3/s)
    1527      392885 :          U_Fe_tot     , &  ! actual iron uptake   (umol/m^3/s)
    1528      392885 :          nitrif       , &  ! nitrification (mmol/m^3/s)
    1529      392885 :          mort_N       , &  ! total algal mortality (mmol N/m^3)
    1530      392885 :          mort_C       , &  ! total algal mortality (mmol C/m^3)
    1531      392885 :          graze_N      , &  ! total algae grazed (mmol N/m^3)
    1532      392885 :          graze_C      , &  ! total algae grazed (mmol C/m^3)
    1533      392885 :          exude_C      , &  ! total carbon exuded by algae (mmol C/m^3)
    1534      392885 :          resp_N       , &  ! total N in respiration (mmol N/m^3)
    1535      392885 :          growth_N          ! total algal growth (mmol N/m^3)
    1536             : !        fr_graze_p   , &  ! fraction of N grazed that becomes protein
    1537             : !                          !  (rest is assimilated) < (1-fr_graze_a)
    1538             : !                          !  and fr_graze_a*fr_graze_e becomes ammonia
    1539             : !        fr_mort_p         ! fraction of N mortality that becomes protein 
    1540             : !                          ! < (1-fr_mort2min)
    1541             : 
    1542             :       real (kind=dbl_kind), dimension(n_algae) :: &
    1543     3996212 :          resp     , &  ! respiration (mmol/m^3/s)
    1544     3996212 :          graze    , &  ! grazing (mmol/m^3/s)
    1545     3996212 :          mort          ! sum of mortality and excretion (mmol/m^3/s)
    1546             : 
    1547             : !  source terms underscore s, removal underscore r
    1548             : 
    1549             :       real (kind=dbl_kind), dimension(n_algae) :: &
    1550     3996212 :          N_s       , &  ! net algal nitrogen sources (mmol/m^3)
    1551     3996212 :          N_r            ! net algal nitrogen removal (mmol/m^3)
    1552             : 
    1553             :       real (kind=dbl_kind), dimension(n_doc) :: &
    1554     3603327 :          DOC_r      , &  ! net DOC removal (mmol/m^3)
    1555     3603327 :          DOC_s           ! net DOC sources (mmol/m^3)
    1556             : 
    1557             :       real (kind=dbl_kind), dimension(n_don) :: &
    1558     3210442 :          DON_r      , &  ! net DON removal (mmol/m^3)
    1559     3210442 :          DON_s           ! net DON sources (mmol/m^3)
    1560             : 
    1561             :       real (kind=dbl_kind), dimension(n_fed) :: &
    1562     3210442 :          Fed_r_l     , &  ! removal due to loss of binding saccharids (nM)
    1563     3210442 :          Fed_r       , &  ! net Fed removal (nM)
    1564     3210442 :          Fed_s       , &  ! net Fed sources (nM)
    1565     3210442 :          rFed             ! ratio of dissolved Fe to tot Fed
    1566             : 
    1567             :       real (kind=dbl_kind), dimension(n_fep) :: &
    1568     3210442 :          Fep_r       , &  ! net Fep removal (nM)
    1569     3210442 :          Fep_s       , &  ! net Fep sources (nM)
    1570     3210442 :          rFep             ! ratio of particulate Fe to tot Fep
    1571             : 
    1572             :       real (kind=dbl_kind) :: &
    1573      392885 :          dN        , &  ! change in N (mmol/m^3)
    1574             : !        N_s_p     , &  ! algal nitrogen photosynthesis (mmol/m^3)
    1575             : !        N_r_g     , &  ! algal nitrogen losses to grazing (mmol/m^3)
    1576             : !        N_r_r     , &  ! algal nitrogen losses to respiration (mmol/m^3)
    1577             : !        N_r_mo    , &  ! algal nitrogen losses to mortality (mmol/m^3)
    1578      392885 :          Nit_s_n   , &  ! nitrate from nitrification (mmol/m^3)
    1579             : !        Nit_s_r   , &  ! nitrate from respiration (mmol/m^3)
    1580      392885 :          Nit_r_p   , &  ! nitrate uptake by algae (mmol/m^3)
    1581      392885 :          Nit_s     , &  ! net nitrate sources (mmol/m^3)
    1582      392885 :          Nit_r     , &  ! net nitrate removal (mmol/m^3)
    1583      392885 :          Am_s_e    , &  ! ammonium source from excretion (mmol/m^3)
    1584      392885 :          Am_s_r    , &  ! ammonium source from respiration (mmol/m^3)
    1585      392885 :          Am_s_mo   , &  ! ammonium source from mort/remin (mmol/m^3) 
    1586      392885 :          Am_r_p    , &  ! ammonium uptake by algae (mmol/m^3)
    1587      392885 :          Am_s      , &  ! net ammonium sources (mmol/m^3)
    1588      392885 :          Am_r      , &  ! net ammonium removal (mmol/m^3)
    1589      392885 :          Sil_r_p   , &  ! silicon uptake by algae (mmol/m^3)
    1590      392885 :          Sil_r     , &  ! net silicon removal (mmol/m^3)
    1591      392885 :          Fe_r_p         ! iron uptake by algae  (nM)
    1592             : !        DOC_r_c   , &  ! net doc removal from bacterial consumption (mmol/m^3)
    1593             : !        doc_s_m   , &  ! protein source due to algal mortality (mmol/m^3)
    1594             : !        doc_s_g        ! protein source due to grazing (mmol/m^3)         
    1595             : 
    1596             :       real (kind=dbl_kind) :: &
    1597      392885 :          DMSPd_s_r , &  ! skl dissolved DMSP from respiration (mmol/m^3)
    1598      392885 :          DMSPd_s_mo, &  ! skl dissolved DMSP from MBJ algal mortexc (mmol/m^3)
    1599      392885 :          DMSPd_r   , &  ! skl dissolved DMSP conversion (mmol/m^3) DMSPD_sk_r
    1600      392885 :          DMSPd_s   , &  ! net skl dissolved DMSP sources (mmol/m^3)
    1601      392885 :          DMS_s_c   , &  ! skl DMS source via conversion (mmol/m^3)
    1602      392885 :          DMS_r_o   , &  ! skl DMS losses due to oxidation (mmol/m^3)
    1603      392885 :          DMS_s     , &  ! net skl DMS sources (mmol/m^3)
    1604      392885 :          DMS_r     , &  ! net skl DMS removal (mmol/m^3)
    1605      392885 :          Fed_tot   , &  ! total dissolved iron from all sources (nM)
    1606      392885 :          Fed_tot_r , &  ! total dissolved iron losses (nM)
    1607      392885 :          Fed_tot_s , &  ! total dissolved iron sources (nM)
    1608      392885 :          Fep_tot   , &  ! total particulate iron from all sources (nM)
    1609             : !        Fep_tot_r , &  ! total particulate iron losses (nM)
    1610      392885 :          Fep_tot_s , &  ! total particulate iron sources (nM)
    1611      392885 :          Zoo_s_a   , &  ! N Losses due to zooplankton assimilation (mmol/m^3)
    1612      392885 :          Zoo_s_s   , &  ! N Losses due to grazing spillage (mmol/m^3)
    1613      392885 :          Zoo_s_m   , &  ! N Losses due to algal mortality (mmol/m^3)
    1614      392885 :          Zoo_s_b        ! N losses due to bacterial recycling of DON (mmol/m^3)
    1615             : 
    1616             :       character(len=*),parameter :: subname='(algal_dyn)'
    1617             : 
    1618             :       !-----------------------------------------------------------------------
    1619             :       ! Initialize
    1620             :       !-----------------------------------------------------------------------
    1621             : 
    1622     1605221 :        conserve_N = .true.
    1623     6420884 :        Nin(:)     = c0
    1624             : !      Cin(:)     = c0
    1625     6420884 :        chlin(:)   = c0
    1626     4815663 :        DOCin(:)   = c0
    1627             : !      DICin(:)   = c0
    1628     3210442 :        DONin(:)   = c0
    1629     3210442 :        Fedin(:)   = c0
    1630     3210442 :        Fepin(:)   = c0
    1631     1605221 :        Nitin      = c0
    1632     1605221 :        Amin       = c0
    1633     1605221 :        Silin      = c0
    1634             : !      DMSPpin    = c0
    1635     1605221 :        DMSPdin    = c0
    1636     1605221 :        DMSin      = c0
    1637     1605221 :        U_Am_tot   = c0
    1638     1605221 :        U_Nit_tot  = c0
    1639     1605221 :        U_Sil_tot  = c0
    1640     1605221 :        U_Fe_tot   = c0
    1641     6420884 :        U_Am_f(:)  = c0
    1642     6420884 :        U_Nit_f(:) = c0
    1643     6420884 :        U_Sil_f(:) = c0
    1644     6420884 :        U_Fe_f(:)  = c0
    1645     4815663 :        DOC_s(:)   = c0
    1646     4815663 :        DOC_r(:)   = c0
    1647             : !      DOC_r_c    = c0
    1648     1605221 :        nitrif     = c0 
    1649     1605221 :        mort_N     = c0
    1650     1605221 :        mort_C     = c0
    1651     1605221 :        graze_N    = c0
    1652     1605221 :        graze_C    = c0
    1653     1605221 :        exude_C    = c0
    1654     1605221 :        resp_N     = c0
    1655     1605221 :        growth_N   = c0
    1656     1605221 :        Nit_r      = c0 
    1657     1605221 :        Am_s       = c0
    1658     1605221 :        Am_r       = c0 
    1659     1605221 :        Sil_r      = c0
    1660     3210442 :        Fed_r(:)   = c0
    1661     3210442 :        Fed_s(:)   = c0
    1662     3210442 :        Fep_r(:)   = c0
    1663     3210442 :        Fep_s(:)   = c0
    1664     1605221 :        DMSPd_s    = c0 
    1665     1605221 :        dTemp      = min(T_bot-T_max,c0)
    1666     1605221 :        Fed_tot    = c0
    1667     1605221 :        Fed_tot_r  = c0
    1668     1605221 :        Fed_tot_s  = c0
    1669     3210442 :        rFed(:)    = c0
    1670     1605221 :        Fep_tot    = c0
    1671             : !      Fep_tot_r  = c0
    1672     1605221 :        Fep_tot_s  = c0
    1673     3210442 :        rFep(:)    = c0
    1674             :      
    1675     1605221 :        Nitin     = ltrcrn(nlt_bgc_Nit)
    1676     1605221 :        op_dep = c0
    1677     6420884 :        do k = 1, n_algae
    1678     4815663 :           Nin(k)   = ltrcrn(nlt_bgc_N(k))
    1679     4815663 :           chlin(k) = R_chl2N(k)* Nin(k)  
    1680     6420884 :           op_dep = op_dep + chlabs(k)*chlin(k)
    1681             :        enddo
    1682     1605221 :        if (tr_bgc_C)   then
    1683             :         ! do k = 1, n_algae
    1684             :         !     Cin(k)=  ltrcrn(nlt_bgc_C(k))
    1685             :         ! enddo
    1686     4815663 :          do k = 1, n_doc
    1687     4815663 :              DOCin(k)= ltrcrn(nlt_bgc_DOC(k))
    1688             :          enddo
    1689             : !        do k = 1, n_dic
    1690             : !            DICin(k)= ltrcrn(nlt_bgc_DIC(k))
    1691             : !        enddo
    1692             :        endif
    1693     1605221 :        if (tr_bgc_Am)  Amin     = ltrcrn(nlt_bgc_Am)
    1694     1605221 :        if (tr_bgc_Sil) Silin    = ltrcrn(nlt_bgc_Sil)
    1695     1605221 :        if (tr_bgc_DMS) then
    1696             :        !      DMSPpin  = ltrcrn(nlt_bgc_DMSPp)
    1697     1605221 :              DMSPdin  = ltrcrn(nlt_bgc_DMSPd)
    1698     1605221 :              DMSin    = ltrcrn(nlt_bgc_DMS) 
    1699             :        endif
    1700             : !      if (tr_bgc_PON) then
    1701             : !         PONin    = c0 
    1702             : !         PONin    = ltrcrn(nlt_bgc_PON) 
    1703             : !      endif
    1704     1605221 :        if (tr_bgc_DON) then
    1705     3210442 :          do k = 1, n_don
    1706     3210442 :              DONin(k) = ltrcrn(nlt_bgc_DON(k))
    1707             :          enddo
    1708             :        endif
    1709     1605221 :        if (tr_bgc_Fe ) then
    1710      148698 :          do k = 1, n_fed 
    1711      148698 :              Fedin(k) = ltrcrn(nlt_bgc_Fed(k))
    1712             :          enddo
    1713      148698 :          do k = 1, n_fep 
    1714      148698 :              Fepin(k) = ltrcrn(nlt_bgc_Fep(k))
    1715             :          enddo
    1716             :        endif
    1717             : 
    1718             :       !-----------------------------------------------------------------------
    1719             :       ! Total iron from all pools
    1720             :       !-----------------------------------------------------------------------
    1721             : 
    1722     3210442 :        do k = 1,n_fed
    1723     3210442 :          Fed_tot = Fed_tot + Fedin(k)
    1724             :        enddo
    1725     3210442 :        do k = 1,n_fep
    1726     3210442 :          Fep_tot = Fep_tot + Fepin(k)
    1727             :        enddo
    1728     1605221 :        if (Fed_tot > puny) then
    1729      148698 :        do k = 1,n_fed
    1730      148698 :          rFed(k) = Fedin(k)/Fed_tot
    1731             :        enddo
    1732             :        endif
    1733     1605221 :        if (Fep_tot > puny) then
    1734      148698 :        do k = 1,n_fep
    1735      148698 :          rFep(k) = Fepin(k)/Fep_tot
    1736             :        enddo
    1737             :        endif
    1738             : 
    1739             :       !-----------------------------------------------------------------------
    1740             :       ! Light limitation  (op_dep) defined above
    1741             :       !-----------------------------------------------------------------------
    1742             : 
    1743     1605221 :        if (op_dep > op_dep_min .and. .not. dEdd_algae) then
    1744     1584270 :          Iavg_loc = fswthru * (c1 - exp(-op_dep)) / op_dep
    1745             :        else
    1746       20951 :          Iavg_loc = fswthru
    1747             :        endif
    1748             : 
    1749     6420884 :        do k = 1, n_algae
    1750             :           ! With light inhibition ! Maybe include light inhibition for diatoms but phaeocystis
    1751    19262652 :           L_lim = (c1 - exp(-alpha2max_low(k)*Iavg_loc)) * exp(-beta2max(k)*Iavg_loc)
    1752             : 
    1753             :           ! Without light inhibition
    1754             :           ! L_lim(k) = (c1 - exp(-alpha2max_low(k)*Iavg_loc))
    1755             : 
    1756             :       !-----------------------------------------------------------------------
    1757             :       ! Nutrient limitation
    1758             :       !-----------------------------------------------------------------------
    1759             : 
    1760     4815663 :           Nit_lim(k) = Nitin/(Nitin + K_Nit(k))
    1761     4815663 :           Am_lim(k)  = c0
    1762     4815663 :           N_lim(k) = Nit_lim(k)
    1763     4815663 :           if (tr_bgc_Am) then
    1764     4815663 :              Am_lim(k) = Amin/(Amin + K_Am(k))
    1765     4815663 :              N_lim(k)  = min(c1, Nit_lim(k) + Am_lim(k))  
    1766             :           endif
    1767     4815663 :           Sil_lim(k) = c1
    1768     4815663 :           if (tr_bgc_Sil .and. K_Sil(k) > c0) Sil_lim(k) = Silin/(Silin + K_Sil(k))
    1769             : 
    1770             :       !-----------------------------------------------------------------------
    1771             :       ! Iron limitation
    1772             :       !-----------------------------------------------------------------------
    1773             : 
    1774     4815663 :           Fe_lim(k) = c1         
    1775     4815663 :           if (tr_bgc_Fe  .and. K_Fe (k) > c0) Fe_lim (k) = Fed_tot/(Fed_tot + K_Fe(k))
    1776             :   
    1777             :       !----------------------------------------------------------------------------
    1778             :       ! Growth and uptake computed within the bottom layer 
    1779             :       ! Note here per A93 discussions and MBJ model, salinity is a universal 
    1780             :       ! restriction.  Comparison with available column nutrients inserted 
    1781             :       ! but in tests had no effect.
    1782             :       ! Primary production reverts to SE form, see MBJ below and be careful
    1783             :       !----------------------------------------------------------------------------
    1784             : 
    1785     4815663 :           growmax_N(k) = mu_max(k) / secday * exp(grow_Tdep(k) * dTemp)* Nin(k) *fsal
    1786     4815663 :           grow_N(k)    = min(L_lim(k), N_lim(k), Sil_lim(k), Fe_lim(k)) * growmax_N(k)
    1787             : !         potU_Nit(k)  = Nit_lim(k)* growmax_N(k)
    1788     4815663 :           potU_Am(k)   = Am_lim(k)* growmax_N(k) 
    1789     4815663 :           U_Am(k)      = min(grow_N(k), potU_Am(k))
    1790     4815663 :           U_Nit(k)     = grow_N(k) - U_Am(k)
    1791     4815663 :           U_Sil(k)     = R_Si2N(k) * grow_N(k)
    1792     4815663 :           U_Fe (k)     = R_Fe2N(k) * grow_N(k)
    1793             : 
    1794     4815663 :           U_Am_tot     = U_Am_tot  + U_Am(k)
    1795     4815663 :           U_Nit_tot    = U_Nit_tot + U_Nit(k)
    1796     4815663 :           U_Sil_tot    = U_Sil_tot + U_Sil(k)
    1797     6420884 :           U_Fe_tot     = U_Fe_tot  + U_Fe(k)
    1798             :        enddo
    1799     6420884 :        do k = 1, n_algae
    1800     4815663 :           if (U_Am_tot > c0) U_Am_f(k) = U_Am(k)/U_Am_tot
    1801     4815663 :           if (U_Nit_tot > c0) U_Nit_f(k) = U_Nit(k)/U_Nit_tot
    1802     4815663 :           if (U_Sil_tot > c0) U_Sil_f(k) = U_Sil(k)/U_Sil_tot
    1803     6420884 :           if (U_Fe_tot > c0) U_Fe_f(k) = U_Fe(k)/U_Fe_tot
    1804             :        enddo
    1805             : 
    1806     1605221 :        if (tr_bgc_Sil) U_Sil_tot = min(U_Sil_tot, max_loss * Silin/dt)
    1807     1605221 :        if (tr_bgc_Fe)  U_Fe_tot  = min(U_Fe_tot, max_loss * Fed_tot/dt)
    1808     1605221 :        U_Nit_tot = min(U_Nit_tot, max_loss * Nitin/dt)  
    1809     1605221 :        U_Am_tot  = min(U_Am_tot,  max_loss * Amin/dt)    
    1810             : 
    1811     6420884 :        do k = 1, n_algae
    1812     4815663 :           U_Am(k)  = U_Am_f(k)*U_Am_tot
    1813     4815663 :           U_Nit(k) = U_Nit_f(k)*U_Nit_tot
    1814     4815663 :           U_Sil(k) = U_Sil_f(k)*U_Sil_tot
    1815     4815663 :           U_Fe(k)  = U_Fe_f(k)*U_Fe_tot
    1816             : 
    1817     4815663 :           if (R_Si2N(k) > c0) then
    1818     1605221 :              grow_N(k) = min(U_Sil(k)/R_Si2N(k),U_Nit(k) + U_Am(k), U_Fe(k)/R_Fe2N(k))
    1819             :           else
    1820     3210442 :              grow_N(k) = min(U_Nit(k) + U_Am(k),U_Fe(k)/R_Fe2N(k))
    1821             :           endif
    1822             : 
    1823     4815663 :           fr_Am(k) = c0
    1824     4815663 :           if (tr_bgc_Am) then
    1825     4815663 :              fr_Am(k) = p5
    1826     4815663 :              if (grow_N(k) > c0) fr_Am(k) = min(U_Am(k)/grow_N(k), c1)
    1827             :           endif
    1828     4815663 :           fr_Nit(k) = c1 - fr_Am(k)
    1829     4815663 :           U_Nit(k)  = fr_Nit(k) * grow_N(k)
    1830     4815663 :           U_Am(k)   = fr_Am(k)  * grow_N(k)
    1831     4815663 :           U_Sil(k)  = R_Si2N(k) * grow_N(k)
    1832     4815663 :           U_Fe (k)  = R_Fe2N(k) * grow_N(k)
    1833             :     
    1834             :       !-----------------------------------------------------------------------
    1835             :       ! Define reaction terms
    1836             :       !-----------------------------------------------------------------------
    1837             : 
    1838             :       ! Since the framework remains incomplete at this point,
    1839             :       ! it is assumed as a starting expedient that 
    1840             :       ! DMSP loss to melting results in 10% conversion to DMS
    1841             :       ! which is then given a ten day removal constant.
    1842             :       ! Grazing losses are channeled into rough spillage and assimilation
    1843             :       ! then following ammonia there is some recycling.
    1844             : 
    1845             :       !--------------------------------------------------------------------
    1846             :       ! Algal reaction term
    1847             :       ! N_react = (grow_N*(c1 - fr_graze-fr_resp) - mort)*dt  
    1848             :       !--------------------------------------------------------------------
    1849             : 
    1850     4815663 :           resp(k)   = fr_resp  * grow_N(k)  
    1851     4815663 :           graze(k)  = fr_graze(k) * grow_N(k)
    1852     1178655 :           mort(k)   = min(max_loss * Nin(k)/dt, &
    1853     4815663 :                           mort_pre(k)*exp(mort_Tdep(k)*dTemp) * Nin(k)/secday)
    1854             :  
    1855             :         ! history variables
    1856     4815663 :           grow_alg(k) = grow_N(k)
    1857     4815663 :           upNOn(k) = U_Nit(k)
    1858     4815663 :           upNHn(k) = U_Am(k)
    1859             : 
    1860             : !         N_s_p  = grow_N(k) * dt  
    1861             : !         N_r_g  = graze(k)  * dt 
    1862             : !         N_r_r  = resp(k)   * dt
    1863             : !         N_r_mo = mort(k)   * dt
    1864     4815663 :           N_s(k)    = (c1- fr_resp - fr_graze(k)) * grow_N(k) *dt   !N_s_p
    1865     4815663 :           N_r(k)    = mort(k) * dt                                  !N_r_g  + N_r_mo + N_r_r 
    1866             : 
    1867     4815663 :           graze_N   = graze_N + graze(k)
    1868     4815663 :           graze_C   = graze_C + R_C2N(k)*graze(k)
    1869     4815663 :           mort_N    = mort_N + mort(k)      
    1870     4815663 :           mort_C    = mort_C + R_C2N(k)*mort(k)
    1871     4815663 :           resp_N    = resp_N + resp(k)
    1872     6420884 :           growth_N  = growth_N + grow_N(k)
    1873             :  
    1874             :       enddo ! n_algae
    1875             :       !--------------------------------------------------------------------
    1876             :       ! Ammonium source: algal grazing, respiration, and mortality
    1877             :       !--------------------------------------------------------------------
    1878             : 
    1879     1605221 :       Am_s_e  = graze_N*(c1-fr_graze_s)*fr_graze_e*dt
    1880     1605221 :       Am_s_mo = mort_N*fr_mort2min*dt
    1881     1605221 :       Am_s_r  = resp_N*dt
    1882     1605221 :       Am_s    = Am_s_r + Am_s_e + Am_s_mo
    1883             : 
    1884             :       !--------------------------------------------------------------------
    1885             :       ! Nutrient net loss terms: algal uptake
    1886             :       !--------------------------------------------------------------------
    1887             :         
    1888     6420884 :       do k = 1, n_algae
    1889     4815663 :          Am_r_p  = U_Am(k)   * dt
    1890     4815663 :          Am_r    = Am_r + Am_r_p
    1891     4815663 :          Nit_r_p = U_Nit(k)  * dt
    1892     4815663 :          Nit_r   = Nit_r + Nit_r_p
    1893     4815663 :          Sil_r_p = U_Sil(k) * dt
    1894     4815663 :          Sil_r   = Sil_r + Sil_r_p
    1895     4815663 :          Fe_r_p  = U_Fe (k) * dt
    1896     4815663 :          Fed_tot_r = Fed_tot_r + Fe_r_p
    1897     6420884 :          exude_C = exude_C + k_exude(k)* R_C2N(k)*Nin(k) / secday
    1898             :       enddo
    1899             : 
    1900             :       !--------------------------------------------------------------------
    1901             :       ! nitrification
    1902             :       !--------------------------------------------------------------------
    1903             : 
    1904     1605221 :       nitrif  = k_nitrif /secday * Amin
    1905     1605221 :       Am_r    = Am_r +  nitrif*dt
    1906     1605221 :       Nit_s_n = nitrif * dt  ! source from NH4
    1907     1605221 :       Nit_s   = Nit_s_n
    1908             : 
    1909             :       !--------------------------------------------------------------------
    1910             :       ! PON:  currently using PON to shadow nitrate
    1911             :       !
    1912             :       ! N Losses are counted in Zoo.  These arise from mortality not 
    1913             :       ! remineralized (Zoo_s_m), assimilated grazing not excreted (Zoo_s_a), 
    1914             :       !spilled N not going to DON (Zoo_s_s) and  bacterial recycling 
    1915             :       ! of DON (Zoo_s_b). 
    1916             :       !--------------------------------------------------------------------
    1917             : 
    1918     1605221 :       if (tr_bgc_Am) then
    1919     1605221 :          Zoo_s_a = graze_N*(c1-fr_graze_e)*(c1-fr_graze_s) *dt
    1920     1605221 :          Zoo_s_s = graze_N*fr_graze_s*dt
    1921     1605221 :          Zoo_s_m = mort_N*dt  -  Am_s_mo
    1922             :       else
    1923           0 :          Zoo_s_a = graze_N*dt*(c1-fr_graze_s)
    1924           0 :          Zoo_s_s = graze_N*fr_graze_s*dt
    1925           0 :          Zoo_s_m = mort_N*dt 
    1926             :       endif
    1927             : 
    1928     1605221 :       Zoo_s_b = c0
    1929             : 
    1930             :       !--------------------------------------------------------------------
    1931             :       ! DON (n_don = 1)
    1932             :       ! Proteins   
    1933             :       !--------------------------------------------------------------------
    1934             : 
    1935     3210442 :       DON_r(:) = c0
    1936     3210442 :       DON_s(:) = c0
    1937             : 
    1938     1605221 :       if (tr_bgc_DON) then
    1939     3210442 :       do n = 1, n_don
    1940     1605221 :          DON_r(n) = kn_bac(n)/secday * DONin(n) * dt
    1941     1605221 :          DON_s(n) = graze_N*f_don(n)*fr_graze_s * dt
    1942     1605221 :          Zoo_s_s  = Zoo_s_s - DON_s(n)
    1943     3210442 :          Zoo_s_b  = Zoo_s_b + DON_r(n)*(c1-f_don_Am(n))
    1944             :          !Am_s     = Am_s + DON_r(n)*f_don_Am(n)
    1945             :       enddo
    1946             :       endif
    1947             :      
    1948     1605221 :       Zoo = Zoo_s_a + Zoo_s_s + Zoo_s_m + Zoo_s_b
    1949             : 
    1950             :       !--------------------------------------------------------------------
    1951             :       ! DOC
    1952             :       ! polysaccharids, lipids
    1953             :       !--------------------------------------------------------------------
    1954             : 
    1955     4815663 :       do n = 1, n_doc
    1956             :           
    1957     3210442 :          DOC_r(n) = k_bac(n)/secday * DOCin(n) * dt
    1958      785770 :          DOC_s(n) = f_doc(n)*(fr_graze_s *graze_C + mort_C)*dt &
    1959     4815663 :                   + f_exude(n)*exude_C
    1960             :       enddo
    1961             : 
    1962             :       !--------------------------------------------------------------------
    1963             :       ! Iron sources from remineralization  (follows ammonium but reduced)
    1964             :       ! only Fed_s(1)  has remineralized sources
    1965             :       !--------------------------------------------------------------------
    1966             :       
    1967     1605221 :       Fed_s(1) = Fed_s(1) + Am_s * R_Fe2N(1) * fr_dFe   ! remineralization source
    1968             : 
    1969             :       !--------------------------------------------------------------------
    1970             :       !  Conversion to dissolved Fe from Particulate requires DOC(1)
    1971             :       !  Otherwise the only source of dFe is from remineralization
    1972             :       !--------------------------------------------------------------------
    1973             : 
    1974     1605221 :       if (tr_bgc_C .and. tr_bgc_Fe) then
    1975       74349 :         if (DOCin(1) > c0) then 
    1976      148698 :         if (Fed_tot/DOCin(1) > max_dfe_doc1) then             
    1977           0 :           do n = 1,n_fed                                  ! low saccharid:dFe ratio leads to
    1978           0 :              Fed_r_l(n)  = Fedin(n)/t_iron_conv*dt/secday ! loss of bioavailable Fe to particulate fraction
    1979           0 :              Fep_tot_s   = Fep_tot_s + Fed_r_l(n)
    1980           0 :              Fed_r(n)    = Fed_r_l(n)                     ! removal due to particulate scavenging
    1981             :           enddo
    1982           0 :           do n = 1,n_fep
    1983           0 :              Fep_s(n) = rFep(n)* Fep_tot_s                ! source from dissolved Fe
    1984             :           enddo
    1985       74349 :         elseif (Fed_tot/DOCin(1) < max_dfe_doc1) then  
    1986      148698 :           do n = 1,n_fep                                  ! high saccharid:dFe ratio leads to
    1987       74349 :              Fep_r(n)  = Fepin(n)/t_iron_conv*dt/secday   ! gain of bioavailable Fe from particulate fraction
    1988      148698 :              Fed_tot_s = Fed_tot_s + Fep_r(n)
    1989             :           enddo  
    1990      148698 :           do n = 1,n_fed
    1991      148698 :              Fed_s(n) = Fed_s(n) + rFed(n)* Fed_tot_s     ! source from particulate Fe
    1992             :           enddo    
    1993             :         endif
    1994             :         endif !Docin(1) > c0
    1995     1530872 :       elseif (tr_bgc_Fe) then
    1996           0 :         do n = 1,n_fed
    1997           0 :            Fed_r(n) = Fed_r(n) + rFed(n)*Fed_tot_r        ! scavenging + uptake
    1998             :         enddo 
    1999             : 
    2000             :       ! source from algal mortality/grazing and fraction of remineralized nitrogen that does 
    2001             :       ! not become immediately bioavailable
    2002             : 
    2003           0 :          do n = 1,n_fep
    2004           0 :             Fep_s(n) = Fep_s(n) + rFep(n)* (Am_s * R_Fe2N(1) * (c1-fr_dFe))   
    2005             :          enddo ! losses not direct to Fed 
    2006             :       endif
    2007             : 
    2008             :       !--------------------------------------------------------------------
    2009             :       ! Sulfur cycle begins here
    2010             :       !--------------------------------------------------------------------
    2011             :       ! Grazing losses are channeled into rough spillage and assimilation
    2012             :       ! then onward and the MBJ mortality channel is included
    2013             :       ! It is assumed as a starting expedient that 
    2014             :       ! DMSP loss to melting gives partial conversion to DMS in product layer
    2015             :       ! which then undergoes Stefels removal.
    2016             : 
    2017             :       !--------------------------------------------------------------------
    2018             :       ! DMSPd  reaction term  (DMSPd conversion is outside the algal loop)
    2019             :       ! DMSPd_react = R_S2N*((fr_graze_s+fr_excrt_2S*fr_graze_e*fr_graze_a)
    2020             :       !                      *fr_graze*grow_N + fr_mort2min*mort)*dt
    2021             :       !             - [\DMSPd]/t_sk_conv*dt
    2022             :       !--------------------------------------------------------------------
    2023     6420884 :       do k = 1,n_algae
    2024     4815663 :          DMSPd_s_r = fr_resp_s  * R_S2N(k) * resp(k)   * dt  !respiration fraction to DMSPd
    2025     4815663 :          DMSPd_s_mo= fr_mort2min * R_S2N(k)* mort(k)   * dt  !mortality and extracellular excretion
    2026     6420884 :          DMSPd_s = DMSPd_s + DMSPd_s_r + DMSPd_s_mo 
    2027             :       enddo
    2028     1605221 :       DMSPd_r = (c1/t_sk_conv) * (c1/secday)  * (DMSPdin) * dt
    2029             : 
    2030             :       !--------------------------------------------------------------------
    2031             :       ! DMS reaction term + DMSPd loss term 
    2032             :       ! DMS_react = ([\DMSPd]*y_sk_DMS/t_sk_conv - c1/t_sk_ox *[\DMS])*dt
    2033             :       !--------------------------------------------------------------------
    2034             : 
    2035     1605221 :       DMS_s_c = y_sk_DMS * DMSPd_r
    2036     1605221 :       DMS_r_o = DMSin * dt / (t_sk_ox * secday)
    2037     1605221 :       DMS_s   = DMS_s_c
    2038     1605221 :       DMS_r   = DMS_r_o
    2039             : 
    2040             :       !-----------------------------------------------------------------------
    2041             :       ! Load reaction array
    2042             :       !-----------------------------------------------------------------------
    2043             : 
    2044     1605221 :       dN = c0
    2045     6420884 :       do k = 1,n_algae
    2046     4815663 :          reactb(nlt_bgc_N(k))  = N_s(k) - N_r(k)
    2047     6420884 :          dN = dN + reactb(nlt_bgc_N(k))
    2048             :       enddo
    2049     1605221 :       if (tr_bgc_C) then
    2050             :        ! do k = 1,n_algae
    2051             :        !      reactb(nlt_bgc_C(k))  = R_C2N(k)*reactb(nlt_bgc_N(k))
    2052             :        ! enddo
    2053     4815663 :          do k = 1,n_doc
    2054     4815663 :             reactb(nlt_bgc_DOC(k))= DOC_s(k) - DOC_r(k)
    2055             :          enddo
    2056             :       endif
    2057     1605221 :             reactb(nlt_bgc_Nit)   = Nit_s   - Nit_r
    2058     1605221 :             dN = dN + reactb(nlt_bgc_Nit)
    2059     1605221 :       if (tr_bgc_Am)  then
    2060     1605221 :             reactb(nlt_bgc_Am)    = Am_s    - Am_r
    2061     1605221 :             dN = dN + reactb(nlt_bgc_Am)
    2062             :       endif
    2063     1605221 :       if (tr_bgc_Sil) then
    2064     1605221 :             reactb(nlt_bgc_Sil)   =  - Sil_r
    2065             :       endif
    2066     1605221 :       if (tr_bgc_DON) then
    2067     3210442 :          do k = 1,n_don
    2068     1605221 :             reactb(nlt_bgc_DON(k))= DON_s(k) - DON_r(k)  
    2069     3210442 :             dN = dN + reactb(nlt_bgc_DON(k))
    2070             :          enddo
    2071             :       endif
    2072     1605221 :       if (tr_bgc_Fe ) then
    2073      148698 :          do k = 1,n_fed
    2074      148698 :             reactb(nlt_bgc_Fed(k))= Fed_s (k) - Fed_r (k) 
    2075             :          enddo
    2076      148698 :          do k = 1,n_fep
    2077      148698 :             reactb(nlt_bgc_Fep(k))= Fep_s (k) - Fep_r (k) 
    2078             :          enddo
    2079             :       endif
    2080     1605221 :       if (tr_bgc_DMS) then
    2081     1605221 :          reactb(nlt_bgc_DMSPd) = DMSPd_s - DMSPd_r
    2082     1605221 :          reactb(nlt_bgc_DMS)   = DMS_s   - DMS_r
    2083             :       endif
    2084     1605221 :       Nerror = dN + Zoo
    2085             :       ! if (abs(Nerror) > max(reactb(:))*1.0e-5) then
    2086             :       !      conserve_N = .false.
    2087             :       !      write(warnstr,*) subname, 'Conservation error!'
    2088             :       !      call icepack_warnings_add(warnstr)
    2089             :       !      write(warnstr,*) subname, 'Nerror,dN, DONin(1),kn_bac(1),secday,dt,n_doc'
    2090             :       !      call icepack_warnings_add(warnstr)
    2091             :       !      write(warnstr,*) subname, Nerror,dN, DONin(1),kn_bac(1),secday,dt,n_doc
    2092             :       !      call icepack_warnings_add(warnstr)
    2093             :       !      write(warnstr,*) subname, 'reactb(nlt_bgc_Nit),reactb(nlt_bgc_N(1)),reactb(nlt_bgc_N(2)'
    2094             :       !      call icepack_warnings_add(warnstr)
    2095             :       !      write(warnstr,*) subname, reactb(nlt_bgc_Nit),reactb(nlt_bgc_N(1)),reactb(nlt_bgc_N(2))
    2096             :       !      call icepack_warnings_add(warnstr)
    2097             :       !      write(warnstr,*) subname, 'reactb(nlt_bgc_Am),reactb(nlt_bgc_DON(1)), DON_r(1),DON_s(1)'
    2098             :       !      call icepack_warnings_add(warnstr)
    2099             :       !      write(warnstr,*) subname, reactb(nlt_bgc_Am),reactb(nlt_bgc_DON(1)),DON_r(1),DON_s(1)
    2100             :       !      call icepack_warnings_add(warnstr)
    2101             :       !      write(warnstr,*) subname, 'Zoo:',Zoo
    2102             :       ! endif
    2103             :           
    2104     1605221 :       end subroutine algal_dyn
    2105             : 
    2106             : !=======================================================================
    2107             : !
    2108             : ! Find ice-ocean flux when ice is thin and internal dynamics/reactions are
    2109             : !             assumed to be zero
    2110             : !
    2111             : ! authors     Nicole Jeffery, LANL
    2112             : 
    2113          16 :       subroutine thin_ice_flux (hin, hin_old, Cin, flux_o_tot, &
    2114             :                                 source, dt, nblyr, &
    2115             :                                 ocean_bio) 
    2116             : 
    2117             :       integer (kind=int_kind), intent(in) :: &
    2118             :          nblyr    ! number of bio layers
    2119             : 
    2120             :       real (kind=dbl_kind), dimension(nblyr+1), intent(inout) :: &
    2121             :          Cin               ! initial concentration*hin_old*phin
    2122             : 
    2123             :       real (kind=dbl_kind), intent(in) :: &
    2124             :          hin_old   , &     ! brine thickness (m) 
    2125             :          hin       , &     ! new brine thickness (m)
    2126             :          dt        , &     ! time step
    2127             :          source    , &     ! atm, ocean, dust flux (mmol/m^2)
    2128             :          ocean_bio         ! ocean tracer concentration (mmol/m^3)
    2129             : 
    2130             :       real (kind=dbl_kind), intent(inout) :: &
    2131             :          flux_o_tot        ! tracer flux, gravity+molecular drainage flux ,
    2132             :                            ! and boundary flux to ocean (mmol/m^2/s)  
    2133             :                            ! positive into the ocean  
    2134             : 
    2135             :      ! local variables
    2136             : 
    2137             :      integer (kind=int_kind) :: &
    2138             :          k                 ! vertical biology layer index
    2139             :    
    2140             :      real (kind=dbl_kind) :: &
    2141          16 :          sum_bio,   & ! initial bio mass (mmol/m^2)
    2142          16 :          zspace,    & ! 1/nblyr
    2143          16 :          dC,        & ! added ocean bio mass (mmol/m^2)
    2144          16 :          dh           ! change in thickness  (m)  
    2145             :    
    2146             :      character(len=*),parameter :: subname='(thin_ice_flux)'
    2147             : 
    2148          16 :      zspace = c1/real(nblyr,kind=dbl_kind)
    2149             : 
    2150          16 :      dC = c0
    2151          16 :      sum_bio = c0
    2152          16 :      dh = hin-hin_old
    2153             :  
    2154          16 :      if (dh .le. c0) then  ! keep the brine concentration fixed
    2155          16 :        sum_bio = (Cin(1)+Cin(nblyr+1))/hin_old*zspace*p5
    2156          16 :        Cin(1) = Cin(1)/hin_old*hin 
    2157          16 :        Cin(nblyr+1) = Cin(nblyr+1)/hin_old*hin
    2158         112 :        do k = 2, nblyr
    2159          96 :         sum_bio = sum_bio + Cin(k)/hin_old*zspace
    2160         112 :         Cin(k) = Cin(k)/hin_old*hin + dC 
    2161             :        enddo
    2162             :      else
    2163           0 :        dC = dh*ocean_bio
    2164           0 :        do k = 1, nblyr+1
    2165           0 :          Cin(k) = Cin(k) + dC
    2166             :        enddo
    2167             :      endif
    2168             :      
    2169          16 :      flux_o_tot = - dh*sum_bio/dt - dC/dt + source/dt 
    2170             : 
    2171          16 :      end subroutine thin_ice_flux
    2172             : 
    2173             : !=======================================================================
    2174             : !
    2175             : ! Compute matrix elements for the low order solution of FEM-FCT scheme
    2176             : ! Predictor step
    2177             : !
    2178             : ! July, 2014 by N. Jeffery
    2179             : !
    2180    11878020 :       subroutine compute_FCT_matrix  (C_in, sbdiag, dt,  nblyr,     &
    2181    11878020 :                                       diag, spdiag, rhs, bgrid,     &
    2182             :                                       darcyV, dhtop, dhbot, &
    2183     3959340 :                                       iphin_N, iDin, hbri_old,      &
    2184     3959340 :                                       atm_add, bphin_N,             &
    2185             :                                       C_top, C_bot, Qbot, Qtop,     &
    2186             :                                       Sink_bot, Sink_top,           &
    2187     3959340 :                                       D_sbdiag, D_spdiag, ML)
    2188             : 
    2189             :       integer (kind=int_kind), intent(in) :: &
    2190             :          nblyr           ! number of bio layers
    2191             : 
    2192             :       real (kind=dbl_kind), dimension(nblyr+1), intent(in) :: &
    2193             :          C_in            ! Initial (bulk) concentration*hbri_old (mmol/m^2)
    2194             :                          ! conserved quantity on i_grid
    2195             : 
    2196             :       real (kind=dbl_kind), intent(in) :: &
    2197             :          dt              ! time step
    2198             :      
    2199             :       real (kind=dbl_kind), dimension (nblyr+1), intent(inout) :: &
    2200             :          iDin            ! Diffusivity on the igrid (1/s)
    2201             : 
    2202             :       real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: &
    2203             :          iphin_N         ! Porosity with min condition on igrid
    2204             : 
    2205             :       real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
    2206             :          bphin_N, &      ! Porosity with min condition on igrid
    2207             :          bgrid 
    2208             : 
    2209             :       real (kind=dbl_kind), dimension (nblyr+1), &
    2210             :          intent(out) :: &
    2211             :          sbdiag      , & ! sub-diagonal matrix elements
    2212             :          diag        , & ! diagonal matrix elements
    2213             :          spdiag      , & ! super-diagonal matrix elements
    2214             :          rhs         , & ! rhs of tri-diagonal matrix eqn.
    2215             :          ML,           & ! lumped mass matrix
    2216             :          D_spdiag, D_sbdiag     ! artificial diffusion matrix
    2217             : 
    2218             :       real (kind=dbl_kind), intent(in) :: &
    2219             :          dhtop         , & ! Change in brine top (m)
    2220             :          dhbot         , & ! Change in brine bottom (m)
    2221             :          hbri_old      , & ! brine height (m)
    2222             :          atm_add       , & ! atm-ice flux
    2223             :          C_top         , & ! bulk surface source (mmol/m^2)
    2224             :          C_bot         , & ! bulk bottom source (mmol/m^2)
    2225             :          darcyV            ! Darcy velocity (m/s
    2226             : 
    2227             :       real (kind=dbl_kind), intent(inout) :: &   ! positive into ice
    2228             :          Qtop         , & ! top  flux source (mmol/m^2/s)
    2229             :          Qbot         , & ! bottom flux  source (mmol/m^2/s)
    2230             :          Sink_bot     , & ! rest of bottom flux (sink or source) (mmol/m^2/s)
    2231             :          Sink_top         ! rest oftop flux (sink or source) (mmol/m^2/s)
    2232             : 
    2233             :       ! local variables
    2234             : 
    2235             :       real (kind=dbl_kind) :: &
    2236     1857000 :          vel, vel2, dphi_dx, zspace
    2237             : 
    2238             :       integer (kind=int_kind) :: &
    2239             :          k                ! vertical index
    2240             : 
    2241             :       real (kind=dbl_kind), dimension (nblyr+1) :: &
    2242    21846180 :          Q_top,  Q_bot,              & ! surface and bottom source
    2243    29274180 :          K_diag, K_spdiag, K_sbdiag, & ! advection matrix
    2244    29274180 :          S_diag, S_spdiag, S_sbdiag, & ! diffusion matrix
    2245    21846180 :          D_diag, iDin_phi
    2246             : 
    2247             :       real (kind=dbl_kind), dimension (nblyr) :: &
    2248    20917680 :          kvector1, kvectorn1
    2249             : 
    2250             :       character(len=*),parameter :: subname='(compute_FCT_matrix)'
    2251             : 
    2252             : !---------------------------------------------------------------------
    2253             : !  Diag (jj) solve for j = 1:nblyr+1
    2254             : !  spdiag(j) == (j,j+1) solve for j = 1:nblyr otherwise 0
    2255             : !  sbdiag(j) == (j,j-1) solve for j = 2:nblyr+1 otherwise 0
    2256             : !---------------------------------------------------------------------
    2257    31674720 :       kvector1(:) = c0
    2258     3959340 :       kvector1(1) = c1
    2259    31674720 :       kvectorn1(:) = c1
    2260     3959340 :       kvectorn1(1) = c0
    2261             : 
    2262     3959340 :       zspace = c1/real(nblyr,kind=dbl_kind)
    2263     3959340 :       Qbot = c0
    2264     3959340 :       Qtop = c0
    2265     3959340 :       Sink_bot = c0
    2266     3959340 :       Sink_top = c0
    2267             : 
    2268             : ! compute the lumped mass matrix
    2269             : 
    2270    35634060 :       ML(:) = zspace
    2271     3959340 :       ML(1) = zspace/c2
    2272     3959340 :       ML(nblyr+1) = zspace/c2
    2273             : 
    2274             : ! compute matrix K: K_diag , K_sbdiag, K_spdiag
    2275             : ! compute matrix S: S_diag, S_sbdiag, S_spdiag
    2276             : 
    2277    35634060 :       K_diag(:) = c0
    2278    35634060 :       D_diag(:) = c0
    2279    35634060 :       D_spdiag(:) = c0
    2280    35634060 :       D_sbdiag(:) = c0
    2281    35634060 :       K_spdiag(:) = c0
    2282    35634060 :       K_sbdiag(:) = c0
    2283    35634060 :       S_diag(:) = c0
    2284    35634060 :       S_spdiag(:) = c0
    2285    35634060 :       S_sbdiag(:) = c0
    2286    35634060 :       iDin_phi(:) = c0
    2287             : 
    2288     3959340 :       iDin_phi(1) = c0  !element 1
    2289     3959340 :       iDin_phi(nblyr+1) = iDin(nblyr+1)/iphin_N(nblyr+1)  !outside ice at bottom
    2290     3959340 :       iDin_phi(nblyr) = p5*(iDin(nblyr+1)/iphin_N(nblyr+1)+iDin(nblyr)/iphin_N(nblyr))
    2291             : 
    2292     3959340 :       vel = (bgrid(2)*dhbot - (bgrid(2)-c1)*dhtop)/dt+darcyV/bphin_N(2)
    2293     3959340 :       K_diag(1) = p5*vel/hbri_old   
    2294     3959340 :       dphi_dx = (iphin_N(nblyr+1) - iphin_N(nblyr))/(zspace)
    2295     3959340 :       vel = (bgrid(nblyr+1)*dhbot - (bgrid(nblyr+1)-c1)*dhtop)/dt  +darcyV/bphin_N(nblyr+1)  
    2296     3959340 :       vel = vel/hbri_old   
    2297     3959340 :       vel2 = (dhbot/hbri_old/dt +darcyV/hbri_old) 
    2298     1857000 :       K_diag(nblyr+1) =   min(c0, vel2) - iDin_phi(nblyr+1)/(zspace+ grid_o/hbri_old) &
    2299     5816340 :                    + p5*(-vel + iDin_phi(nblyr)/bphin_N(nblyr+1)*dphi_dx) 
    2300             : 
    2301    27715380 :       do k = 1, nblyr-1
    2302    23756040 :          vel = (bgrid(k+1)*dhbot - (bgrid(k+1)-c1)*dhtop)/dt+darcyV/bphin_N(k+1)
    2303    23756040 :          iDin_phi(k+1) = p5*(iDin(k+2)/iphin_N(k+2) + iDin(k+1)/iphin_N(k+1))
    2304    23756040 :          dphi_dx =  (iphin_N(k+1) - iphin_N(k))/(zspace)
    2305     5571000 :          K_spdiag(k)= p5*(vel/hbri_old - &
    2306    23756040 :                          iDin_phi(k)/(bphin_N(k+1))*dphi_dx)
    2307             : 
    2308    23756040 :          vel = (bgrid(k+1)*dhbot - (bgrid(k+1)-c1)*dhtop)/dt  +darcyV/bphin_N(k+1)   
    2309    23756040 :          dphi_dx = c0
    2310    23756040 :          dphi_dx = kvectorn1(k)*(iphin_N(k+1) - iphin_N(k))/(zspace) 
    2311     5571000 :          K_sbdiag(k+1)= -p5*(vel/hbri_old- &
    2312    29327040 :                          iDin_phi(k)/bphin_N(k+1)*dphi_dx)
    2313    23756040 :          K_diag(k) = K_diag(1)*kvector1(k) + (K_spdiag(k) + K_sbdiag(k))*kvectorn1(k)
    2314             : 
    2315    23756040 :          S_diag(k+1) =   -(iDin_phi(k)+ iDin_phi(k+1))/zspace
    2316    23756040 :          S_spdiag(k)   = iDin_phi(k)/zspace 
    2317    27715380 :          S_sbdiag(k+1) = iDin_phi(k)/zspace
    2318             :       enddo
    2319             : 
    2320             :       ! k = nblyr
    2321             : 
    2322     3959340 :       vel = (bgrid(nblyr+1)*dhbot - (bgrid(nblyr+1)-c1)*dhtop)/dt+darcyV/bphin_N(nblyr+1)
    2323     3959340 :       dphi_dx =  (iphin_N(nblyr+1) - iphin_N(nblyr))/(zspace)
    2324      928500 :       K_spdiag(nblyr)= p5*(vel/hbri_old - &
    2325     4887840 :                          iDin_phi(nblyr)/(bphin_N(nblyr+1))*dphi_dx)
    2326     3959340 :       vel = (bgrid(nblyr+1)*dhbot - (bgrid(nblyr+1)-c1)*dhtop)/dt  +darcyV/bphin_N(nblyr+1)   
    2327     3959340 :       dphi_dx = kvectorn1(nblyr)*(iphin_N(nblyr+1) - iphin_N(nblyr))/(zspace) 
    2328      928500 :       K_sbdiag(nblyr+1)= -p5*(vel/hbri_old- &
    2329     4887840 :                          iDin_phi(nblyr)/bphin_N(nblyr+1)*dphi_dx)
    2330     3959340 :       K_diag(nblyr) = K_spdiag(nblyr) + K_sbdiag(nblyr)
    2331     3959340 :       S_diag(nblyr+1) = -iDin_phi(nblyr)/zspace 
    2332     3959340 :       S_spdiag(nblyr)   = iDin_phi(nblyr)/zspace 
    2333     3959340 :       S_sbdiag(nblyr+1) = iDin_phi(nblyr)/zspace
    2334             :        
    2335             : ! compute matrix artificial D: D_spdiag, D_diag  (D_spdiag(k) = D_sbdiag(k+1))
    2336             : 
    2337    31674720 :       do k = 1,nblyr
    2338    27715380 :          D_spdiag(k)    = max(-K_spdiag(k), c0, -K_sbdiag(k+1))
    2339    31674720 :          D_sbdiag(k+1)  = D_spdiag(k)
    2340             :       enddo
    2341    35634060 :       do k = 1,nblyr+1
    2342    35634060 :          D_diag(k) = D_diag(k) - D_spdiag(k) - D_sbdiag(k)
    2343             :       enddo
    2344             : 
    2345             : ! compute Q_top and Q_bot: top and bottom sources 
    2346             : 
    2347     3959340 :       vel2 = -(dhtop/hbri_old/dt +darcyV/bphin_N(1)/hbri_old)
    2348             : 
    2349    35634060 :       Q_top(:) = c0
    2350     3959340 :       Q_top(1) = max(c0,vel2*C_top + atm_add/dt)
    2351     3959340 :       Qtop = Q_top(1)
    2352             : 
    2353     3959340 :       vel = (dhbot/hbri_old/dt +darcyV/hbri_old)  ! going from iphin_N(nblyr+1) to c1 makes a difference
    2354             : 
    2355    35634060 :       Q_bot(:) = c0
    2356     1857000 :       Q_bot(nblyr+1) = max(c0,vel*C_bot) + iDin_phi(nblyr+1)*C_bot&  
    2357     4887840 :                       /(zspace + grid_o/hbri_old)
    2358             :  
    2359     3959340 :       Qbot = Q_bot(nblyr+1)
    2360             :  
    2361     3959340 :       Sink_bot = K_diag(nblyr+1) +  K_spdiag(nblyr)
    2362     3959340 :       Sink_top = K_diag(1) + K_sbdiag(2)
    2363             : 
    2364             : !compute matrix elements  (1 to nblyr+1)
    2365             : 
    2366    35634060 :      spdiag = -dt * (D_spdiag + K_spdiag + S_spdiag)
    2367    35634060 :      sbdiag = -dt * (D_sbdiag + K_sbdiag + S_sbdiag)
    2368    35634060 :      diag = ML - dt *  (D_diag + K_diag + S_diag)
    2369    35634060 :      rhs = ML * C_in + dt * Q_top + dt* Q_bot 
    2370             :       
    2371     3959340 :      end subroutine compute_FCT_matrix
    2372             : 
    2373             : !=======================================================================
    2374             : !
    2375             : ! Compute matrices for final solution FCT for passive tracers
    2376             : ! Corrector step
    2377             : !
    2378             : ! July, 2014 by N. Jeffery
    2379             : !
    2380     3959340 :       subroutine compute_FCT_corr    (C_in, C_low, dt,  nblyr, &
    2381     3959340 :                                       D_sbdiag, D_spdiag, ML)
    2382             : 
    2383             :       integer (kind=int_kind), intent(in) :: &
    2384             :          nblyr           ! number of bio layers
    2385             : 
    2386             :       real (kind=dbl_kind), dimension(nblyr+1), intent(in) :: &
    2387             :          C_in            ! Initial (bulk) concentration*hbri_old (mmol/m^2)
    2388             :                          ! conserved quantity on igrid
    2389             : 
    2390             :       real (kind=dbl_kind), dimension(nblyr+1), intent(inout) :: &
    2391             :          C_low           ! Low order solution (mmol/m^2) corrected
    2392             : 
    2393             :       real (kind=dbl_kind), intent(in) :: &
    2394             :          dt              ! time step
    2395             : 
    2396             :       real (kind=dbl_kind), dimension (nblyr+1), &
    2397             :          intent(in) :: &
    2398             :          D_sbdiag      , & ! sub-diagonal artificial diffusion matrix elements
    2399             :          ML            , & ! Lumped mass diagonal matrix elements
    2400             :          D_spdiag          ! super-diagonal artificial diffusion matrix elements
    2401             : 
    2402             :       ! local variables
    2403             : 
    2404             :       real (kind=dbl_kind) :: &
    2405      928500 :         zspace
    2406             : 
    2407             :       integer (kind=int_kind) :: &
    2408             :          k                ! vertical index
    2409             : 
    2410             :       real (kind=dbl_kind), dimension (nblyr+1) :: &
    2411    21846180 :          M_spdiag, M_sbdiag,         & ! mass matrix
    2412    29274180 :          F_diag, F_spdiag, F_sbdiag, & ! anti-diffusive matrix
    2413    21846180 :          Pplus, Pminus,              & !
    2414    21846180 :          Qplus, Qminus,              & !
    2415    21846180 :          Rplus, Rminus,              & !
    2416    22774680 :          a_spdiag, a_sbdiag            ! weightings of F
    2417             : 
    2418             :       character(len=*),parameter :: subname='(compute_FCT_corr)'
    2419             : 
    2420             : !---------------------------------------------------------------------
    2421             : !  Diag (jj) solve for j = 1:nblyr+1
    2422             : !  spdiag(j) == (j,j+1) solve for j = 1:nblyr otherwise 0
    2423             : !  sbdiag(j) == (j,j-1) solve for j = 2:nblyr+1 otherwise 0
    2424             : !---------------------------------------------------------------------
    2425             : 
    2426     3959340 :       zspace = c1/real(nblyr,kind=dbl_kind) 
    2427             : 
    2428             : ! compute the mass matrix
    2429             : 
    2430    35634060 :       M_spdiag(:) = zspace/c6
    2431     3959340 :       M_spdiag(nblyr+1) = c0
    2432    35634060 :       M_sbdiag(:) = zspace/c6
    2433     3959340 :       M_sbdiag(1) = c0
    2434             : 
    2435             : ! compute off matrix F:
    2436             : 
    2437    35634060 :       F_diag(:) = c0
    2438    35634060 :       F_spdiag(:) = c0
    2439    35634060 :       F_sbdiag(:) = c0
    2440             : 
    2441    31674720 :       do k = 1, nblyr 
    2442    19498500 :          F_spdiag(k) = M_spdiag(k)*(C_low(k)-C_in(k) - C_low(k+1)+ C_in(k+1))/dt &
    2443    40714380 :                      + D_spdiag(k)*(C_low(k)-C_low(k+1))
    2444    25998000 :          F_sbdiag(k+1) =  M_sbdiag(k+1)*(C_low(k+1)-C_in(k+1) - C_low(k)+ C_in(k))/dt &
    2445    53713380 :                        + D_sbdiag(k+1)*(C_low(k+1)-C_low(k))
    2446             : 
    2447    27715380 :          if (F_spdiag(k)*(C_low(k) - C_low(k+1)) > c0) F_spdiag(k) = c0
    2448    31674720 :          if (F_sbdiag(k+1)*(C_low(k+1) - C_low(k)) > c0) F_sbdiag(k+1) = c0
    2449             :       enddo
    2450             : 
    2451    35634060 :       if (maxval(abs(F_spdiag)) > c0) then
    2452             : 
    2453             : ! compute the weighting factors: a_spdiag, a_sbdiag
    2454             : 
    2455    22248675 :       a_spdiag(:) = c0
    2456    22248675 :       a_sbdiag(:) = c0
    2457             : 
    2458     2472075 :       Pplus(1)  = max(c0, F_spdiag(1))  
    2459     2472075 :       Pminus(1) = min(c0, F_spdiag(1))
    2460     2472075 :       Pplus(nblyr+1)  = max(c0, F_sbdiag(nblyr+1))
    2461     2472075 :       Pminus(nblyr+1) = min(c0, F_sbdiag(nblyr+1))
    2462     2472075 :       Qplus(1) = max(c0,C_low(2)-C_low(1))
    2463     2472075 :       Qminus(1)= min(c0,C_low(2)-C_low(1))
    2464     2472075 :       Qplus(nblyr+1) = max(c0,C_low(nblyr)-C_low(nblyr+1))
    2465     2472075 :       Qminus(nblyr+1)= min(c0,C_low(nblyr)-C_low(nblyr+1))
    2466     2472075 :       Rplus(1)  = min(c1, ML(1)*Qplus(1)/dt/(Pplus(1)+puny))
    2467     2472075 :       Rminus(1) = min(c1, ML(1)*Qminus(1)/dt/(Pminus(1)-puny))
    2468     2472075 :       Rplus(nblyr+1)  = min(c1, ML(nblyr+1)*Qplus(nblyr+1)/dt/(Pplus(nblyr+1)+puny))
    2469     2472075 :       Rminus(nblyr+1) = min(c1, ML(nblyr+1)*Qminus(nblyr+1)/dt/(Pminus(nblyr+1)-puny))
    2470    17304525 :       do k = 2,nblyr
    2471    14832450 :          Pplus(k)  = max(c0,F_spdiag(k)) + max(c0,F_sbdiag(k))
    2472    14832450 :          Pminus(k) = min(c0,F_spdiag(k)) + min(c0,F_sbdiag(k))
    2473    14832450 :          Qplus(k)  = max(c0, max(C_low(k+1)-C_low(k),C_low(k-1)-C_low(k)))
    2474    14832450 :          Qminus(k) = min(c0, min(C_low(k+1)-C_low(k),C_low(k-1)-C_low(k)))
    2475    14832450 :          Rplus(k)  = min(c1, ML(k)*Qplus(k)/dt/(Pplus(k)+puny))
    2476    17304525 :          Rminus(k) = min(c1, ML(k)*Qminus(k)/dt/(Pminus(k)-puny))
    2477             :       enddo
    2478             :      
    2479    19776600 :       do k = 1, nblyr 
    2480    17304525 :          a_spdiag(k) = min(Rminus(k),Rplus(k+1))
    2481    17304525 :          if (F_spdiag(k) > c0) a_spdiag(k) = min(Rplus(k),Rminus(k+1))
    2482    17304525 :          a_sbdiag(k+1) = min(Rminus(k+1),Rplus(k))
    2483    19776600 :          if (F_sbdiag(k+1) > c0) a_sbdiag(k+1) = min(Rplus(k+1),Rminus(k))
    2484             :       enddo
    2485             : 
    2486             : !compute F_diag:
    2487             : 
    2488     2472075 :       F_diag(1) = a_spdiag(1)*F_spdiag(1)
    2489     2472075 :       F_diag(nblyr+1) = a_sbdiag(nblyr+1)* F_sbdiag(nblyr+1)
    2490     2472075 :       C_low(1) = C_low(1) + dt*F_diag(1)/ML(1)
    2491     2472075 :       C_low(nblyr+1) = C_low(nblyr+1) + dt*F_diag(nblyr+1)/ML(nblyr+1)
    2492             : 
    2493    17304525 :       do k = 2,nblyr
    2494    14832450 :          F_diag(k) = a_spdiag(k)*F_spdiag(k) + a_sbdiag(k)*F_sbdiag(k)
    2495    17304525 :          C_low(k) = C_low(k) + dt*F_diag(k)/ML(k)
    2496             :       enddo
    2497             :       
    2498             :       endif  !F_spdiag is nonzero
    2499             : 
    2500     3959340 :       end subroutine compute_FCT_corr
    2501             :   
    2502             : !=======================================================================
    2503             : !
    2504             : ! Tridiagonal matrix solver-- for salinity
    2505             : !
    2506             : ! authors William H. Lipscomb, LANL
    2507             : !         C. M. Bitz, UW
    2508             : !
    2509     7918680 :       subroutine tridiag_solverz (nmat,     sbdiag,   &
    2510     7918680 :                                  diag,      spdiag,   &
    2511     3959340 :                                  rhs,       xout)
    2512             : 
    2513             :       integer (kind=int_kind), intent(in) :: &
    2514             :          nmat            ! matrix dimension
    2515             : 
    2516             :       real (kind=dbl_kind), dimension (nmat), intent(in) :: &
    2517             :          sbdiag      , & ! sub-diagonal matrix elements
    2518             :          diag        , & ! diagonal matrix elements
    2519             :          spdiag      , & ! super-diagonal matrix elements
    2520             :          rhs             ! rhs of tri-diagonal matrix eqn.
    2521             : 
    2522             :       real (kind=dbl_kind), dimension (nmat), intent(inout) :: &
    2523             :          xout            ! solution vector
    2524             : 
    2525             :       ! local variables     
    2526             : 
    2527             :       integer (kind=int_kind) :: &
    2528             :          k               ! row counter
    2529             : 
    2530             :       real (kind=dbl_kind) :: &
    2531      928500 :          wbeta           ! temporary matrix variable
    2532             : 
    2533             :       real (kind=dbl_kind), dimension(nmat):: &
    2534    14418180 :          wgamma          ! temporary matrix variable
    2535             : 
    2536             :       character(len=*),parameter :: subname='(tridiag_solverz)'
    2537             : 
    2538     3959340 :       wbeta = diag(1)
    2539     3959340 :       xout(1) = rhs(1) / wbeta
    2540             : 
    2541    31674720 :       do k = 2, nmat
    2542    27715380 :          wgamma(k) = spdiag(k-1) / wbeta
    2543    27715380 :          wbeta = diag(k) - sbdiag(k)*wgamma(k)
    2544    31674720 :          xout(k) = (rhs(k) - sbdiag(k)*xout(k-1)) / wbeta
    2545             :       enddo                     ! k
    2546             : 
    2547    31674720 :       do k = nmat-1, 1, -1
    2548    31674720 :          xout(k) = xout(k) - wgamma(k+1)*xout(k+1)
    2549             :       enddo                     ! k
    2550             : 
    2551     3959340 :       end subroutine tridiag_solverz
    2552             : 
    2553             : !=======================================================================
    2554             : !
    2555             : ! authors     Nicole Jeffery, LANL
    2556             : 
    2557     3959340 :       subroutine check_conservation_FCT (C_init, C_new, C_low, S_top, &
    2558             :                                          S_bot, L_bot, L_top, dt,     &
    2559             :                                          fluxbio, nblyr, &
    2560             :                                          source) 
    2561             : 
    2562             :       integer (kind=int_kind), intent(in) :: &
    2563             :          nblyr             ! number of bio layers
    2564             : 
    2565             :       real (kind=dbl_kind), dimension(nblyr+1), intent(in) :: &
    2566             :          C_init        , & ! initial bulk concentration * h_old (mmol/m^2)
    2567             :          C_new             ! new bulk concentration * h_new (mmol/m^2)
    2568             : 
    2569             :       real (kind=dbl_kind), dimension(nblyr+1), intent(out) :: &
    2570             :          C_low             ! define low order solution = C_new
    2571             : 
    2572             :       real (kind=dbl_kind),  intent(in) :: &
    2573             :          S_top         , & ! surface flux into ice (mmol/m^2/s)
    2574             :          S_bot         , & ! bottom flux into ice (mmol/m^2/s)
    2575             :          L_bot         , & ! remaining  bottom flux into ice (mmol/m^2/s)
    2576             :          L_top         , & ! remaining  top  flux into ice (mmol/m^2/s)
    2577             :          dt            , & 
    2578             :          source            ! nutrient source from snow and atmosphere (mmol/m^2)
    2579             : 
    2580             :       real (kind=dbl_kind), intent(inout) :: &
    2581             :          fluxbio            ! (mmol/m^2/s)  positive down (into the ocean)
    2582             : 
    2583             :       ! local variables
    2584             : 
    2585             :       integer (kind=int_kind) :: &
    2586             :          k
    2587             : 
    2588             :       real (kind=dbl_kind) :: &
    2589      928500 :          diff_dt     , &
    2590      928500 :          C_init_tot  , &
    2591      928500 :          C_new_tot   , &
    2592      928500 :          zspace      , &  !1/nblyr
    2593      928500 :          accuracy          ! centered difference is Order(zspace^2)
    2594             : 
    2595             :       character(len=*),parameter :: subname='(check_conservation_FCT)'
    2596             : 
    2597     3959340 :       zspace = c1/real(nblyr,kind=dbl_kind)
    2598             : 
    2599             :       !-------------------------------------
    2600             :       !  Ocean flux: positive into the ocean
    2601             :       !-------------------------------------
    2602     3959340 :       C_init_tot = (C_init(1) + C_init(nblyr+1))*zspace*p5
    2603     3959340 :       C_new_tot = (C_new(1) + C_new(nblyr+1))*zspace*p5
    2604     3959340 :       C_low(1) = C_new(1)
    2605     3959340 :       C_low(nblyr+1) = C_new(nblyr+1)
    2606             : 
    2607    27715380 :       do k = 2, nblyr
    2608    23756040 :          C_init_tot = C_init_tot + C_init(k)*zspace
    2609    23756040 :          C_new_tot = C_new_tot + C_new(k)*zspace
    2610    27715380 :          C_low(k) = C_new(k)
    2611             :       enddo
    2612             : 
    2613     3959340 :       accuracy = 1.0e-14_dbl_kind*max(c1, C_init_tot, C_new_tot)
    2614     3959340 :       fluxbio = (C_init_tot - C_new_tot + source)/dt
    2615     3959340 :       diff_dt =C_new_tot - C_init_tot - (S_top+S_bot+L_bot*C_new(nblyr+1)+L_top*C_new(1))*dt
    2616             : 
    2617    35634060 :       if (minval(C_low) < c0) then
    2618           0 :          write(warnstr,*) subname, 'Positivity of zbgc low order solution failed: C_low:',C_low
    2619           0 :          call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
    2620             :       endif
    2621             :            
    2622     3959340 :       if (abs(diff_dt) > accuracy ) then
    2623           0 :          call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
    2624           0 :          write(warnstr,*) subname, 'Conservation of zbgc low order solution failed: diff_dt:',&
    2625           0 :                         diff_dt
    2626           0 :          write(warnstr,*) subname, 'Total initial tracer', C_init_tot
    2627           0 :          write(warnstr,*) subname, 'Total final1  tracer', C_new_tot
    2628           0 :          write(warnstr,*) subname, 'bottom final tracer', C_new(nblyr+1)
    2629           0 :          write(warnstr,*) subname, 'top final tracer', C_new(1)
    2630           0 :          write(warnstr,*) subname, 'Near bottom final tracer', C_new(nblyr)
    2631           0 :          write(warnstr,*) subname, 'Near top final tracer', C_new(2)
    2632           0 :          write(warnstr,*) subname, 'Top flux*dt into ice:', S_top*dt
    2633           0 :          write(warnstr,*) subname, 'Bottom flux*dt into ice:', S_bot*dt
    2634           0 :          write(warnstr,*) subname, 'Remaining bot flux*dt into ice:', L_bot*C_new(nblyr+1)*dt
    2635           0 :          write(warnstr,*) subname, 'S_bot*dt + L_bot*C_new(nblyr+1)*dt'
    2636           0 :          write(warnstr,*) subname,  S_bot*dt + L_bot*C_new(nblyr+1)*dt
    2637           0 :          write(warnstr,*) subname, 'fluxbio*dt:', fluxbio*dt
    2638           0 :          write(warnstr,*) subname, 'fluxbio:', fluxbio
    2639           0 :          write(warnstr,*) subname, 'Remaining top flux*dt into ice:', L_top*C_new(1)*dt
    2640             :       endif
    2641             :          
    2642     3959340 :       end subroutine check_conservation_FCT
    2643             : 
    2644             : !=======================================================================
    2645             : 
    2646             : ! For each grid cell, sum field over all ice and snow layers
    2647             : !
    2648             : ! author: Nicole Jeffery, LANL
    2649             : 
    2650           0 :       subroutine bgc_column_sum (nblyr, nslyr, hsnow, hbrine, xin, xout)
    2651             : 
    2652             :       integer (kind=int_kind), intent(in) :: &
    2653             :          nblyr, &         ! number of ice layers
    2654             :          nslyr            ! number of snow layers
    2655             : 
    2656             :       real (kind=dbl_kind), dimension(nblyr+3), intent(in) :: &
    2657             :          xin              ! input field
    2658             : 
    2659             :       real (kind=dbl_kind), intent(in) :: &
    2660             :          hsnow, &         ! snow thickness
    2661             :          hbrine           ! brine height 
    2662             : 
    2663             :       real (kind=dbl_kind), intent(out) :: &
    2664             :          xout             ! output field
    2665             : 
    2666             :       ! local variables
    2667             : 
    2668             :       real (kind=dbl_kind) :: &
    2669           0 :          dzssl, &        ! snow surface layer (m)
    2670           0 :          dzint, &        ! snow interior depth (m)
    2671           0 :          hslyr, &        ! snow layer thickness (m)
    2672           0 :          zspace          ! brine layer thickness/hbrine
    2673             : 
    2674             :       integer (kind=int_kind) :: &
    2675             :          n                ! category/layer index
    2676             : 
    2677             :       character(len=*),parameter :: subname='(bgc_column_sum)'
    2678             : 
    2679           0 :       hslyr      = hsnow/real(nslyr,kind=dbl_kind)
    2680           0 :       dzssl      = min(hslyr*p5, hs_ssl)
    2681           0 :       dzint      = max(c0,hsnow - dzssl)
    2682           0 :       zspace     = c1/real(nblyr,kind=dbl_kind)  
    2683             : 
    2684           0 :       xout = c0
    2685           0 :       xout = (xin(1) + xin(nblyr+1))*hbrine*p5*zspace
    2686           0 :       do n = 2, nblyr
    2687           0 :          xout = xout + xin(n)*zspace*hbrine
    2688             :       enddo                 ! n
    2689           0 :       xout = xout + dzssl*xin(nblyr+2) + dzint*xin(nblyr+3)
    2690             : 
    2691           0 :       end subroutine bgc_column_sum
    2692             : 
    2693             : !=======================================================================
    2694             : 
    2695             :       end module icepack_algae
    2696             : 
    2697             : !=======================================================================

Generated by: LCOV version 1.14-6-g40580cd