! *************************************************************************************************
! biogem_box.f90
! C-GOLDSTEIn/BioGeM
! MISCELLANEOUS MECHANICS OF THE SYSTEM
! *************************************************************************************************


MODULE biogem_box

  use gem_carbchem
  
  USE biogem_lib
  IMPLICIT NONE
  SAVE


CONTAINS


  ! ******************************************
  ! *** OCEAN-ATMOSPHERE EXCHANGE ROUTINES ***
  ! ******************************************


  ! *** calculate solubility coefficients ***
  ! <<< NON-GENERIC ALGORITHM >>>
  subroutine sub_calc_solconst(dum_i,dum_j)
    ! dummy arguments
    integer,INTENT(in)::dum_i,dum_j
    ! local variables
    integer::l,ia
    REAL::loc_T,loc_rT,loc_Tr100,loc_S
    real::loc_rho,junkst
    ! calculate local constants
    ! NOTE: pressure in units of (bar) (1 m depth approx = 1 dbar pressure)
    ! NOTE: temperature in K
    loc_T     = ocn(io_T,dum_i,dum_j,n_kmax)

!km#ifdef tneg5sol
!km       junkst = loc_T - 5.0
!km       if (junkst > tfsw(n_kmax) ) then
!km          loc_T = junkst
!km       else
!km          loc_T = tfsw(n_kmax)
!km       endif
!km#endif

    !km 4/2012 - solubility perturbation temp is now a biogem_config.par parameter; default=0 of course
    !km print*,'in biogem_box; what is solubility_dtemp, solubility_dsal? ', solubility_dtemp, solubility_dsal
    !km print*,'  loc_T ', loc_T

    junkst = loc_T + solubility_dtemp
    if (junkst > tfsw(n_kmax) ) then
       loc_T = junkst
    else
       loc_T = tfsw(n_kmax)
    endif
    !km print*,'  junkst, loc_T ',junkst, loc_T 

#ifdef solTconst
!seasonal constant 5/12/09
       if (biostep_test == 0 ) then
          print*,'error on seasonal previous state read'
          print*,'biostep_stop = ',biostep_test
          stop
       endif

       loc_T = ocn_T_season1(dum_i,dum_j,n_kmax,biostep)
#endif

    if (loc_T <  const_zeroC)         loc_T = const_zeroC
    if (loc_T > (const_zeroC + 50.0)) loc_T = const_zeroC + 50.0

    loc_S = ocn(io_S,dum_i,dum_j,n_kmax)
    !km 4/2012 - solubility salinity perturbation as solubility_dtemp; biogem_config.par parameter; defaul=1.0
    loc_S = loc_S * solubility_dsal

    if (loc_S < 20.0) loc_S = 20.0
    if (loc_S > 40.0) loc_S = 40.0
    loc_rT    = 1.0/loc_T
    loc_Tr100 = loc_T/100.0
    loc_rho   = phys_ocn(ipo_rho,dum_i,dum_j,n_kmax)
    ! calculate Solubility Coefficients (mol/(kg atm))
    ! NOTE: for CO2 and N2O, the soluability coefficient is in units of mol/(kg atm)
    !       rather than as a Bunsen Solubility Coefficient (see Wanninkohf [1992])
    !       => convert units for others
    ! NOTE: for CFC-11 and CFC-12, the soluability coefficient is in units of mol/(kg atm)
    !       rather than as a Bunsen Solubility Coefficient (see Wanninkohf [1992])
    !       (actaully, it is not really this simple and K should be corrected for water vapour pressure and lame things like that)
    DO l=3,n_iamax
       ia = conv_iselected_ia(l)
       IF (atm_type(ia) == 1) then
          SELECT CASE (ia)
          CASE (ia_pCO2,ia_pN2O)
             ocnatm_airsea_solconst(ia,dum_i,dum_j) = EXP( &
                  & par_bunsen_coef(1,ia) + par_bunsen_coef(2,ia)*(100*loc_rT) + par_bunsen_coef(3,ia)*LOG(loc_Tr100) + &
                  & loc_S*(par_bunsen_coef(4,ia) + par_bunsen_coef(5,ia)*(loc_Tr100) + par_bunsen_coef(6,ia)*(loc_Tr100)**2) &
                  &  )
          CASE (ia_pCFC11,ia_pCFC12)
             ocnatm_airsea_solconst(ia,dum_i,dum_j) = EXP( &
                  & par_bunsen_coef(1,ia) + par_bunsen_coef(2,ia)*(100*loc_rT) + par_bunsen_coef(3,ia)*LOG(loc_Tr100) + &
                  & loc_S*(par_bunsen_coef(4,ia) + par_bunsen_coef(5,ia)*(loc_Tr100) + par_bunsen_coef(6,ia)*(loc_Tr100)**2) &
                  &  )
          CASE default
             ocnatm_airsea_solconst(ia,dum_i,dum_j) = EXP( &
                  & par_bunsen_coef(1,ia) + par_bunsen_coef(2,ia)*(100*loc_rT) + par_bunsen_coef(3,ia)*LOG(loc_Tr100) + &
                  & loc_S*(par_bunsen_coef(4,ia) + par_bunsen_coef(5,ia)*(loc_Tr100) + par_bunsen_coef(6,ia)*(loc_Tr100)**2) &
                  &  )/ &
                  & (loc_rho*const_V)
          END SELECT
       end if
    end do
  end subroutine sub_calc_solconst


  ! *** calculate piston velocity ***
  ! <<< GENERIC ALGORITHM >>>
  SUBROUTINE sub_calc_pv(dum_i,dum_j)
    ! dummy arguments
    INTEGER::dum_i,dum_j
    ! local variables
    integer::l,ia
    REAL::loc_Sc
    REAL::loc_TC,loc_TC2,loc_TC3
    real::loc_u2, junkst
    ! set local variables
    ! temperature powers
    ! NOTE: temeprature must be converted to the correct units (degrees C)
    ! NOTE: valid temperature range is 0 - 30 C for the Schmidt number empirical fit - see; Wanninkhof et al. [1992]
    loc_TC  = ocn(io_T,dum_i,dum_j,n_kmax) - const_zeroC
#ifdef tneg5sol
       junkst = loc_TC - 5.0
       if (junkst > tfsw(n_kmax) ) then
          loc_TC = junkst
       else
          loc_TC = tfsw(n_kmax)
       endif
#endif
#ifdef solTconst
!seasonal constant 5/12/09
       if (biostep_test == 0 ) then
          print*,'error on seasonal previous state read'
          print*,'biostep_stop = ',biostep_test
          stop
       endif
  
       loc_TC = ocn_T_season1(dum_i,dum_j,n_kmax,biostep) - const_zeroC
#endif
    if (loc_TC <  0.0) loc_TC =  0.0 
    if (loc_TC > 30.0) loc_TC = 30.0 
    loc_TC2 = loc_TC*loc_TC
    loc_TC3 = loc_TC2*loc_TC
    ! wind speed^2
    loc_u2 = phys_ocnatm(ipoa_u,dum_i,dum_j)**2
    !  calculate piston velocity
    DO l=3,n_iamax
       ia = conv_iselected_ia(l)
       IF (atm_type(ia) == 1) then
          if (.NOT. ocnatm_airsea_B(ia)) then
             ! calculate gas transfer Schmidt number
             loc_Sc =                           &
                  & par_Sc_coef(1,ia)         - &
                  & par_Sc_coef(2,ia)*loc_TC  + &
                  & par_Sc_coef(3,ia)*loc_TC2 - &
                  & par_Sc_coef(4,ia)*loc_TC3
             ! calculate CO2 gas transfer velocity (piston velocity)
             ! NOTE: from Wanninkhof [1992] equation 1/3
             ! NOTE: convert from units of (cm hr-1) to (m yr-1)
             ! NOTE: pre-calculate 1.0/660 (= 1.515E-3)
             ocnatm_airsea_pv(ia,dum_i,dum_j) = conv_cm_m*conv_yr_hr*par_gastransfer_a*loc_u2*(loc_Sc*1.515E-3)**(-0.5)
          end if
       end if
    end do
  END SUBROUTINE sub_calc_pv


  ! *** calculate air-sea exchange ***
  ! <<< NON GENERIC ALGORITHM >>>
  FUNCTION fun_calc_ocnatm_flux(dum_i,dum_j,dum_atm,dum_dt)
    ! result variable
    REAL,dimension(0:n_atm)::fun_calc_ocnatm_flux ! units of (mol yr-1)
    ! dummy arguments
    INTEGER,INTENT(in)::dum_i,dum_j
    REAL,dimension(0:n_atm),INTENT(in)::dum_atm
    REAL,INTENT(in)::dum_dt
    ! local variables
    integer::l,ia,io
    REAL,dimension(0:n_atm)::loc_focnatm,loc_fatmocn
    real::loc_alpha_k,loc_alpha_alpha,loc_alpha_1,loc_alpha_2,loc_alpha_sum
    real::loc_alpha_sa,loc_alpha_as
    real::loc_rho,loc_TC,junkst
    real::loc_r13C_ocn,loc_r14C_ocn,loc_r15N_ocn
    real::loc_r13C_atm,loc_r14C_atm,loc_r15N_atm
    real::loc_ocn,loc_atm
    real::loc_A
    real::loc_r_dflux_deqm
    REAL,dimension(0:n_atm)::loc_dflux
    REAL,dimension(0:n_ocn)::loc_deqm

    ! *** INITIALIZE VARIABLES ***
    loc_focnatm(:) = 0.0
    loc_fatmocn(:) = 0.0
    loc_rho = phys_ocn(ipo_rho,dum_i,dum_j,n_kmax)
    loc_TC = ocn(io_T,dum_i,dum_j,n_kmax) - const_zeroC
#ifdef tneg5sol
       junkst = loc_TC - 5.0
       if (junkst > tfsw(n_kmax) ) then
          loc_TC = junkst
       else
          loc_TC = tfsw(n_kmax)
       endif
#endif
    ! area available for air-sea gas transfer
    loc_A = (1.0 - phys_ocnatm(ipoa_seaice,dum_i,dum_j))*phys_ocnatm(ipoa_A,dum_i,dum_j)
    ! *** calculate air-sea gas exchange fluxes ***
    DO l=3,n_iamax
       ia = conv_iselected_ia(l)
       if (.NOT. ocnatm_airsea_A(ia)) then   !in general, do this
          ! set corresponding ocean tracer
          ! NOTE: assume that there is a one-to-one mapping from atm tracers to ocn tracers,
          !       with the corresponding ocean tracer index given in the i=1 index position of the conv_atm_ocn_i array
          io = conv_atm_ocn_i(1,ia)
          SELECT CASE (atm_type(ia))
          CASE (1)
             ! calculate bulk gas exchange
             ! NOTE: check for special case of CO2; only CO2(aq) is relevant to air-sea gas exchange
             ! NOTE: local atmospheric tracer value has Bunsen Solubility Coefficient implicit in its value 
             if (io == io_DIC) then
                loc_ocn = carb(ic_conc_CO2,dum_i,dum_j,n_kmax)
             else
                loc_ocn = ocn(io,dum_i,dum_j,n_kmax)
             endif
             loc_atm = ocnatm_airsea_solconst(ia,dum_i,dum_j)*dum_atm(ia)
             ! make sure nothing 'nasty' can happen if a tracer has a -ve concentration
             if (loc_ocn < const_real_nullsmall) loc_ocn = 0.0
             if (loc_atm < const_real_nullsmall) loc_atm = 0.0
             ! calculate gas exchange fluxes ocn->atm and atm->ocn
             ! NOTE: units of (mol yr-1)
             ! NOTE: the soluability coefficient must be converted to units of mol/(kg atm) from a Bunsen Solubility Coefficient
             loc_focnatm(ia) = ocnatm_airsea_pv(ia,dum_i,dum_j)*loc_A*loc_rho*loc_ocn
             loc_fatmocn(ia) = ocnatm_airsea_pv(ia,dum_i,dum_j)*loc_A*loc_rho*loc_atm

             ! check for 'excessive' gas transfer (i.e., with the potential to lead to numerical instability)
             ! => rescale the fluxes so that the ocean surface is brought exactly into equilibrium
             ! NOTE: in the case of DIC, only CO2(aq) is considered
             ! NOTE: no account is taken of molar changes due to ocean circulation and biological activity in making this check
             ! calculate the molar magnitude of ocean deficit or surfit w.r.t. the atmosphere
             loc_deqm(io) = abs(phys_ocn(ipo_dD,dum_i,dum_j,n_kmax)*phys_ocnatm(ipoa_A,dum_i,dum_j)*loc_rho*(loc_atm - loc_ocn))
             ! calculate the molar transfer that would normally then be applied
             loc_dflux(ia) = abs(dum_dt*(loc_focnatm(ia) - loc_fatmocn(ia)))
             ! ensure that molar transfer does not exceed the current disequilibrium
             ! (i.e., ensure that a +ve disequilibrium is not turned into a larger -ve disequilibrium at the next time-step)

          case default
             ! calculate derived isotopic exchange
             ! NOTE: assume that associated bulk tracer flux has already been calculated (i.e., earlier during this routine)
             SELECT CASE (ia)
             CASE (ia_pCO2_13C)
                ! isotopic fluxes - 13C
                loc_r13C_ocn = ocn(io_DIC_13C,dum_i,dum_j,n_kmax)/ocn(io_DIC,dum_i,dum_j,n_kmax)
                loc_r13C_atm = dum_atm(ia_pCO2_13C)/dum_atm(ia_pCO2)
                ! overall 13C fractionationscheme following Marchal et al. [1998]
                ! NOTE: fractionation factors all take from Zhang et al. [1995]
                ! NOTE: some notation borrowed from Yamahaka and Tajika [1996]
                ! kinetic fractionation
                loc_alpha_k = 0.99912
                ! equilibrium fractionation air/sea
                ! fractionation between aqueous CO2 and gaseous CO2
                loc_alpha_alpha = 0.99869 + 4.9E-6 * loc_TC
                ! equilibrium fractionation sea/air
                loc_alpha_1 = (0.99869 + loc_TC*4.9E-6)/(1.01078 - loc_TC*114.0E-6)
                loc_alpha_2 = (0.99869 + loc_TC*4.9E-6)/(1.00722 - loc_TC*52.0E-6)
                loc_alpha_sum = &
                     & ( &
                     &   carb(ic_conc_CO2,dum_i,dum_j,n_kmax) + &
                     &   loc_alpha_1*carb(ic_conc_HCO3,dum_i,dum_j,n_kmax) + &
                     &   loc_alpha_2*carb(ic_conc_CO3,dum_i,dum_j,n_kmax) &
                     & ) &
                     & /ocn(io_DIC,dum_i,dum_j,n_kmax)
                loc_alpha_as = loc_alpha_alpha*loc_alpha_k
                loc_alpha_sa = loc_alpha_sum*loc_alpha_k
                ! calculate fluxesHFliBOmG
                loc_fatmocn(ia_pCO2_13C) = loc_alpha_as*loc_r13C_atm*loc_fatmocn(ia_pCO2)
                loc_focnatm(ia_pCO2_13C) = loc_alpha_sa*loc_r13C_ocn*loc_focnatm(ia_pCO2)
             CASE (ia_pCO2_14C)
                ! isotopic fluxes - 14C
                ! NOTE: assume that 13C fractionations have already been calculated (i.e., earlier during this routine)
                ! NOTE: const_13C_14C_fracratio defines the ratio of fractionation of 14C:12C relative to 13C:12C
                loc_r14C_ocn = ocn(io_DIC_14C,dum_i,dum_j,n_kmax)/ocn(io_DIC,dum_i,dum_j,n_kmax)
                loc_r14C_atm = dum_atm(ia_pCO2_14C)/dum_atm(ia_pCO2)
                ! calculate fluxes
                loc_fatmocn(ia_pCO2_14C) = (loc_alpha_as**2)*loc_r14C_atm*loc_fatmocn(ia_pCO2)
                loc_focnatm(ia_pCO2_14C) = (loc_alpha_sa**2)*loc_r14C_ocn*loc_focnatm(ia_pCO2)
             CASE (ia_pN2_15N)
                ! km isotopic fluxes - 15N of N2; assume no fractionation - 1Mar07
                loc_r15N_ocn = ocn(io_N2_15N,dum_i,dum_j,n_kmax)/ocn(io_N2,dum_i,dum_j,n_kmax)
                loc_r15N_atm = dum_atm(ia_pN2_15N)/dum_atm(ia_pN2)
                ! calculate fluxes
                loc_fatmocn(ia_pN2_15N) = loc_r15N_atm*loc_fatmocn(ia_pN2)
                loc_focnatm(ia_pN2_15N) = loc_r15N_ocn*loc_focnatm(ia_pN2)
             case default
                ! \/\/\/ INSERT CODE TO DEAL WITH ADDITIONAL ISOTOPES \/\/\/
                !
                ! /\/\/\ INSERT CODE TO DEAL WITH ADDITIONAL ISOTOPES /\/\/\
             end SELECT
          end SELECT
          ! calculate net gas transfer and set results variable
          fun_calc_ocnatm_flux(ia) = loc_focnatm(ia) - loc_fatmocn(ia)
       end if
    end do

  END FUNCTION fun_calc_ocnatm_flux


  ! *** calculate (surface ocean) biological tracer uptake**
  ! <<< CONTAINS NON-GENERIC ALGORITHMS >>>
!  SUBROUTINE sub_calc_bio_uptake(dum_i,dum_j,dum_dt,dum_docn_restore,dum_ocn_flux,dum_biostep)
  SUBROUTINE sub_calc_bio_uptake(dum_i,dum_j,dum_dt,dum_docn_restore,dum_ocn_flux)
    ! dummy arguments
    INTEGER,INTENT(in)::dum_i,dum_j
    real,intent(in)::dum_dt
    real,dimension(0:n_ocn,n_maxk),intent(inout)::dum_docn_restore
    real,dimension(0:n_ocn,n_maxk),intent(inout)::dum_ocn_flux
    ! local variables
    INTEGER::k,l,io,is
    integer::loc_i,loc_tot_i,dum_biostep
    real,dimension(n_ocn)::loc_bio_uptake
    real::loc_dPO4,loc_ficefree,loc_kI
    real::loc_dPO4_lg,loc_dPO4_sm,loc_dPO4_dia_used,loc_dPO4_lg_dia,loc_dPO4_lg_nd      
    real::loc_alpha_1,loc_alpha_2,loc_alpha_sum
    real::loc_Kq,loc_delta,loc_alpha,loc_TC
    real::loc_r13C,loc_r14C,loc_r15N
    real::loc_junk1,loc_junk2,junkst
    real::loc_Si_conc                                    !Hannah 2 July 2010 -- Concentration of Si in water
    real::loc_fet_max,loc_fet_min,loc_a1,loc_a2,loc_a3,loc_a4,loc_max_si2n

    integer::loc_mm_index_lg, loc_mm_index_sm
    character*15 limnut

    real::loc_Cd,loc_FeT         !COPY from Andy SUN 2009/02/25

!    ! Michealis-Menten kinetics [PO4] & [NO3] by Chikamoto 05/12/06
    real::po4_MM, no3_MM, co2_MM, fe_MM
    ! Hannah's large and small classes
    real::po4_MM_lg, no3_MM_lg, co2_MM_lg, fe_MM_lg, sio2_MM
    real::po4_MM_sm, no3_MM_sm, co2_MM_sm, fe_MM_sm 
    real::po4_lim, no3_lim, co2_lim, fe_lim, fe_lim_lg, fe_lim_sm, sio2_lim
    real::loc_limiter_lg, loc_limiter_sm
    real::loc_rate_lim_lg, loc_rate_lim_sm 
    real::loc_conc_lim, loc_conc_lim_lg, loc_conc_lim_sm
    real::loc_bio_part_red_POC_POFe_lg                   !Hannah 6 July 2010
    real::loc_bio_part_red_POC_POFe_sm                   !Hannah 6 July 2010
!    real::po4_MM, po4_MM_lg, po4_MM_sm, no3_MM, no3_MM_lg, no3_MM_sm, co2_MM, co2_MM_lg, co2_MM_sm, fe_MM, fe_MM_lg, fe_MM_sm, po4_lim, no3_lim, co2_lim, fe_lim
   ! 30Si fractionation by Chikamoto 05/12/06 
    real::loc_r30Si 

    ! local variables for productivity modifiers - km 07/21/06
    real,parameter::decayI=20.0             ! exponential decay : 20 m
    real,parameter::halfsatI=20.0           ! light half saturation : 20 Wm-2

    real,parameter::c1=1.d0
    real,parameter::c0=0.d0
    real::meanI                             ! average light within 1st layer (Wm-2)
    real::zlayer, z1, z2                    ! dz, top z of grid box, bottom z of grid box
    real::loc_temp                          ! temp dependence of production
    real::loc_biomass,loc_biomass_lg, loc_biomass_sm                       ! biomass dependence (Maier-Reimer, 1993)
    real::loc_biomass_lg_dia,loc_biomass_lg_nd                             !Hannah 2 July 2010
    real::loc_mld                           ! mixed layer depth dependence
    real::dtime                             ! time step

#ifdef wor16_2in75
    real,parameter::zcrit=75.0              ! critical depth as per OCMIP-2 : 75 m
#elif wor16_2in75d
    real,parameter::zcrit=75.0
#elif wor16_2in100
    real,parameter::zcrit=100.0
#elif wor16_2in100ec
    real,parameter::zcrit=100.0
#elif wor16_2in100ecmm
    real,parameter::zcrit=100.0
#elif wor16_2in100d
    real,parameter::zcrit=100.0
#endif

    levels: DO k=n_kmax+1-nlayer_prod,n_kmax            ! k=15,16 = top 2 layers for 2in100

    ! *** INITIALIZE VARIABLES ***
    ! NOTE: kI is AJR's, which I think is not right, unless it is redefined in LMTBtau

    loc_bio_uptake(:) = c0
    loc_TC = ocn(io_T,dum_i,dum_j,k) - const_zeroC

#ifdef tneg5prod
       junkst = loc_TC - 5.0
       if (junkst > tfsw(k) ) then
          loc_TC = junkst
       else
          loc_TC = tfsw(k)
       endif
#endif

    loc_ficefree = (1.0 - phys_ocnatm(ipoa_seaice,dum_i,dum_j))
!kstcrack 070905
!no cracks    if (loc_ficefree < par_crack ) loc_ficefree = par_crack
!kstcrack
!    loc_kI = phys_ocnatm(ipoa_solfor,dum_i,dum_j)/goldstein_solconst   ! AJR's
    po4_MM = c0
    po4_MM_lg = c0
    po4_MM_sm = c0
    no3_MM = c0
    no3_MM_lg = c0
    no3_MM_sm = c0
    co2_MM = c0
    co2_MM_lg = c0
    co2_MM_sm = c0
    fe_MM = c0
    fe_MM_lg = c0
    fe_MM_sm = c0
    sio2_MM = c0
    loc_biomass = c0
    loc_biomass_lg = c0
    loc_biomass_sm = c0
    loc_junk1 = c0
    loc_junk2 = c0
    loc_bio_part_red_POC_POFe_lg = c0  !Hannah 6 July 2010: Do I need this?
    loc_bio_part_red_POC_POFe_sm = c0  !Hannah 6 July 2010: Do I need this?
!    loc_biomass_lg_dia = c0        !Hannah 2 July 2010
!    loc_biomass_lg_nd = c0         !Hannah 2 July 2010
    dtime =dum_dt

    if (LMTBtau) then !moved this up here 9/1/10kst
       ! KM modifiers
       zlayer = goldstein_dz(k)*goldstein_dsc
       if (k==n_kmax) then
          z1 = 0.0
          z2 = zlayer
       else
          z1 = goldstein_dz(n_kmax)*goldstein_dsc
          z2 = zlayer + goldstein_dz(n_kmax)*goldstein_dsc
       endif       
       meanI = phys_ocnatm(ipoa_solfor,dum_i,dum_j)*decayI/zlayer*(exp(-z1/decayI)-exp(-z2/decayI))
       loc_kI = meanI/(meanI+halfsatI)
 
       loc_temp = (loc_TC+2)/(loc_TC+10)
#ifdef prodTconst      !kst      keep temperature constant for only production calculation
       if (biostep_test == 0 ) then
          print*,'error on seasonal previous state read'
          print*,'biostep_stop = ',biostep_test
          stop
       endif
       loc_temp = (ocn_T_season1(dum_i,dum_j,k,biostep)-const_zeroC+2.0 )/( ocn_T_season1(dum_i,dum_j,k,biostep)-const_zeroC+10.0)
#endif

!       loc_temp = min(1.0, en_tempwt*loc_temp)     !kst: use this to alter temperature dependence for ensemble runs
       loc_mld = min(1.0,zcrit/mldz(dum_i,dum_j))
    endif

    ! *** ADJUST PARTICULATE COMPOSITION 'REDFIELD' RATIOS ***
    ! CaCO3:POC 'rain ratio'
    ! NOTE: a correction is made for the fact that a proportion of the POM is transformed into DOM,
    !      whereas the initially calculated CaCO3 and opal fluxes do not change
    !      => re-scale CaCO3 and opal ratios so that the prescribed export ratio value better reflects final export composition
    if (opt_force(iopt_force_CaCO3toPOCrainratio)) then              !usually false
       bio_part_red(is_POC,is_CaCO3,dum_i,dum_j) = (1.0 - par_bio_red_DOMfrac)*par_bio_CaCO3toPOCrainratio(dum_i,dum_j)
      ! print*,'bio_part_red(is_POC,is_CaCO3,10,10)',bio_part_red(is_POC,is_CaCO3,10,10)
    else
#ifdef CO3const
!seasonal constant 5/12/09
       if (biostep_test == 0 ) then
          print*,'error on seasonal previous state read'
          print*,'biostep_stop = ',biostep_test
          stop
       endif
       if (CO3_carb_ohm_season1(dum_i,dum_j,k,biostep) > 1.0) then
          bio_part_red(is_POC,is_CaCO3,dum_i,dum_j) = (1.0 - par_bio_red_DOMfrac)*par_bio_red_POC_CaCO3* &
               & (CO3_carb_ohm_season1(dum_i,dum_j,k,biostep) - 1.0)**par_bio_red_POC_CaCO3_pP
       else
          bio_part_red(is_POC,is_CaCO3,dum_i,dum_j) = 0.0
       endif

#else
       if (carb(ic_ohm_cal,dum_i,dum_j,k) > 1.0) then
          bio_part_red(is_POC,is_CaCO3,dum_i,dum_j) = (1.0 - par_bio_red_DOMfrac)*par_bio_red_POC_CaCO3* &
               & (carb(ic_ohm_cal,dum_i,dum_j,k) - 1.0)**par_bio_red_POC_CaCO3_pP
         ! print*,'bio_part_red(is_POC,is_CaCO3,dum_i,dum_j)',bio_part_red(is_POC,is_CaCO3,dum_i,dum_j)
       else                                                          !if not saturated wrt cal, no CaCO3 is formed
          bio_part_red(is_POC,is_CaCO3,dum_i,dum_j) = 0.0
       end if
#endif          
    end if


       
!!! 7/9/2010 - comment out these lines (opal production) with the introduction of 5N2T w/ Si depletion linked to organic C
!  this opal:POC will get modified in 5N2T (which requires Fe, too)

    if (sed_select(is_opal)) then
       ! modify Si:C uptake according to H4SiO4 depletion
       ! NOTE: this is just to ensure that SiO2 does not drop below zero,
       !       and does not consistute a true 'nutrient limitation' of productivity
       if (ocn(io_SiO2,dum_i,dum_j,k) > const_real_nullsmall) then
          bio_part_red(is_POC,is_opal,dum_i,dum_j) = (1.0 - par_bio_red_DOMfrac)*par_bio_red_POC_opal*  &
               & ocn(io_SiO2,dum_i,dum_j,k)/(par_bio_c0_SiO2 + ocn(io_SiO2,dum_i,dum_j,k))
       else
          bio_part_red(is_POC,is_opal,dum_i,dum_j) = 0.0
       end if
    end if

!COPY from Andy Sun 2009/02/25
   ! modify C:Fe cellular quotient according to Fe limitation
    ! NOTE: following Ridgwell [2001] (mean parameter values from diatom and coccolithophorid parameterizations)
    ! NOTE: default uniform Fe:C ratio has already been set during model initialization
    ! NOTE: for options 'bio_PFeSi*' this has no effect, as below two different ratios will be used for sicliceous
    ! and non-siliceous phytoplankton. Finally an appropriate mixture of these two ratios (according to the ratio
    ! of siliceous and non-siliceous phytoplankton produced) will be used to update the value computed here                needs to be checked....
    ! KEY: par_bio_FetoC_pP == power in [FeT] dependent Fe:C ratio equation [Ridgwell, 2001] (-0.4225)
    !      par_bio_FetoC_K  == scaling in [FeT] dependent Fe:C ratio equation [Ridgwell, 2001] (103684.0)
    !      par_bio_FetoC_C  == constant in [FeT] dependent Fe:C ratio equation [Ridgwell, 2001] (0.0)

!   below a threshold (par_part_red_FeTmin, as [FeT] decr. Fe:C decr (more efficent, more C per Fe)
!NTOE:  par_part_red_FetoCmax is actually CtoFe, C:Fe

    if (sed_select(is_POFe)) then
       loc_FeT = ocn(io_Fe,dum_i,dum_j,k) + ocn(io_FeL,dum_i,dum_j,k)           
       !  for now:  hardwire in, then make these adjustable and read in biobio_5:
!       par_bio_FetoC_C  = 16500
!       par_bio_FetoC_K = 90478.0
!       par_bio_FetoC_pP = -0.65
!NOTE:  for the following parameters, the POC/POFE ratio will not exceed ~340 000 regardless of the max in bio_bio_5, because
!       of the 'usual' surface [FeT] - it doesn't reach less than about 
       par_bio_FetoC_C  = 0.
!!       par_bio_FetoC_K = 94500.
       par_bio_FetoC_K = 103684.
       par_bio_FetoC_pP = -0.4225
       if (.NOT. opt_bio(iopt_bio_Fe_fixedFetoC) ) then   !usually not fixed. ie: this if is done
          if (loc_FeT > par_part_red_FeTmin) then       
             bio_part_red(is_POC,is_POFe,dum_i,dum_j) = 1.0/ &
                  & MIN(par_part_red_FetoCmax, (par_bio_FetoC_C + par_bio_FetoC_K* (1.0E9*loc_FeT)**par_bio_FetoC_pP))    !Hannah
             bio_part_red(is_POFe,is_POC,dum_i,dum_j) = 1.0/bio_part_red(is_POC,is_POFe,dum_i,dum_j)
          else
             bio_part_red(is_POC,is_POFe,dum_i,dum_j) = 1.0/par_part_red_FetoCmax
             bio_part_red(is_POFe,is_POC,dum_i,dum_j) = 1.0/bio_part_red(is_POC,is_POFe,dum_i,dum_j)
          end if
       end if
    end if   

    ! *** CALCULATE ISOTOPIC FRACTIONATION ***
    ! >>> NON-GENERIC ALGORITHM
    ! NOTE: implement isotopic fraction as a 'Redfield' field ratio (populate array <bio_part_red>)
    if (sed_select(is_POC_13C)) then
       ! calculate the 13C/12C fractionation between DIC and POC
       ! NOTE: calculate d14C anyway if d13C is done (as it doesn't involve much additional calculation)
       loc_r13C = ocn(io_DIC_13C,dum_i,dum_j,k)/ocn(io_DIC,dum_i,dum_j,k)
       loc_r14C = ocn(io_DIC_14C,dum_i,dum_j,k)/ocn(io_DIC,dum_i,dum_j,k)
       loc_alpha_1 = (0.99869 + loc_TC*4.9E-6)/(1.01078 - loc_TC*114.0E-6)
       loc_alpha_2 = (0.99869 + loc_TC*4.9E-6)/(1.00722 - loc_TC*52.0E-6)
       loc_alpha_sum = &
            & ( &
            &   carb(ic_conc_CO2,dum_i,dum_j,k) + &
            &   loc_alpha_1*carb(ic_conc_HCO3,dum_i,dum_j,k) + &
            &   loc_alpha_2*carb(ic_conc_CO3,dum_i,dum_j,k) &
            & ) &
            & /ocn(io_DIC,dum_i,dum_j,k)
       loc_Kq = const_d13C_DIC_Corg_Q2_c + const_d13C_DIC_Corg_Q2_x*ocn(io_T,dum_i,dum_j,k) + &
            & const_d13C_DIC_Corg_Q2_x2*ocn(io_T,dum_i,dum_j,k)**2
       loc_delta = -const_d13C_DIC_Corg_ef + &
            & (const_d13C_DIC_Corg_ef - const_d13C_DIC_Corg_ed)*loc_Kq/carb(ic_conc_CO2,dum_i,dum_j,k)
       loc_alpha = loc_alpha_sum*(1.0 + loc_delta/1000.0)
       bio_part_red(is_POC,is_POC_13C,dum_i,dum_j) = loc_alpha*loc_r13C
       ! calculate the 14C/12C fractionation between DIC and POC
       bio_part_red(is_POC,is_POC_14C,dum_i,dum_j) = (loc_alpha**2)*loc_r14C
       ! calculate 13C/12C fractionation between DIC and CaCO3
       ! NOTE: relate isotopic species to their dependent species, not to POC as is done for the type #1 solid sed tracers
       ! T-dependent fractionation for calcite following Mook [1986]
       ! NOTE: approximated by assuming that the d13C of HCO3- is the same as that of total DIC 
       loc_delta = 15.10 - 4232.0/ocn(io_T,dum_i,dum_j,k)
       loc_alpha = 1.0 + loc_delta/1000.0
       bio_part_red(is_CaCO3,is_CaCO3_13C,dum_i,dum_j) = loc_alpha*loc_r13C
       ! calculate 14C/12C fractionation between DIC and CaCO3
       bio_part_red(is_CaCO3,is_CaCO3_14C,dum_i,dum_j) = (loc_alpha**2)*loc_r14C
    end if

    if (sed_select(is_PON_15N)) then 
       ! calculate the 15N/14N fractionation between NO3 and PON 
       ! NOTE; check first for non-zero nitrate concentration to prevent potential unpleasantness ... 
       ! by M.Chikamoto 07-11-2006 considering 15N fraction of NO3_15N through nitrate assimilation 
       !                           PON_15N +5 permil (calculating above) 
       !                           NO3_15N -5 permil 
       loc_r15N = 0. 
       if (ocn(io_NO3,dum_i,dum_j,k) > const_real_nullsmall) then 
          loc_r15N = ocn(io_NO3_15N,dum_i,dum_j,k)/ocn(io_NO3,dum_i,dum_j,k) 
       end if 
       loc_delta = -5.0 !-5.0 ; km 8Mar07 set to 0 to debug
       loc_alpha = 1.0 + loc_delta/1000.0 
 
       bio_part_red(is_PON,is_PON_15N,dum_i,dum_j) = loc_alpha*loc_r15N 
    end if 
    
    if (sed_select(is_opal_30Si)) then
       ! by M. Chikamoto 07-06-2006 
       ! cauclulate the 30Si/28Si fractionation between SiO2 and opal 
       loc_r30Si = 0. 
       if (ocn(io_SiO2,dum_i,dum_j,k) > const_real_nullsmall)then 
          loc_r30Si = ocn(io_SiO2_30Si,dum_i,dum_j,k)/ocn(io_SiO2,dum_i,dum_j,k) 
       endif
       loc_alpha = 0.9989 ! by Wischmeyer et al. 2003 
       bio_part_red(is_Opal,is_Opal_30Si,dum_i,dum_j)  = loc_alpha*loc_r30Si 
    end if
    ! \/\/\/ INSERT CODE TO DEAL WITH ADDITIONAL ISOTOPES \/\/\/
    !
    ! /\/\/\ INSERT CODE TO DEAL WITH ADDITIONAL ISOTOPES /\/\/\

    ! *** CALCULATE BIOLOGICAL POC EXPORT ***
    ! NOTE: production is calculated as the concentration of newly-formed particulate material in the surface ocean layer
    !       that occurs within any single time step
! _MM is rate limiting
! _lim is concentration
    po4_MM = ocn(io_PO4,dum_i,dum_j,k)/(par_bio_c0_PO4 + ocn(io_PO4,dum_i,dum_j,k))  
    no3_MM = ocn(io_NO3,dum_i,dum_j,k)/(par_bio_c0_NO3 + ocn(io_NO3,dum_i,dum_j,k))  
    co2_MM = 10*carb(ic_conc_co2,dum_i,dum_j,k)/(par_bio_c0_CO2 + carb(ic_conc_co2,dum_i,dum_j,k))  
    
    !include stoichiometry here:                                                          1=PO4  2=NO3  3=CO2  
    !     note:  100* factor is added to co2_lim because par_bio_red_POP_POC is DIC/P, not [CO2aq]/P  (co2 ~ 1% of DIC)
    po4_lim = ocn(io_PO4,dum_i,dum_j,k)
    no3_lim = ocn(io_NO3,dum_i,dum_j,k)/par_bio_red_POP_PON
    co2_lim = 100.0*carb(ic_conc_co2,dum_i,dum_j,k)/par_bio_red_POP_POC
    
    SELECT CASE (par_bio_prodopt)
    CASE ('NONE')
       par_bio_k0_PO4 = c0                                   ! par_bio_k0_PO4 = 0.0 = no production
    CASE ('1N1T_PO4restore')
       ! calculate PO4 depletion and filter for positive productivity; indicated by a negative value of <dum_docn_restore>
       par_bio_k0_PO4 = c1                                   ! par_bio_k0_PO4 = 1.0 so growth rate doesn't affect production
       loc_ficefree = c1   !this I think, is not right.
       loc_biomass = c1
       if (.not. Llimit) loc_kI = c1                         
       if (.not. LMTBtau) dtime = c1

       print*,'RESTOREing to [NO3]'
!kst restore to prev_state of [NO3]                             
       if ( restore_prev_state ) then
          if ( dum_docn_restore(io_NO3,k) < -const_real_nullsmall ) then
             po4_MM =  (-dum_docn_restore(io_NO3,k)) /par_bio_red_POP_PON        !molP/kg
          else
             po4_MM = 0.0
          endif
          if ( LMTBtau ) then
             loc_biomass = po4_MM
             po4_MM = c1
          endif
        ! modify net ocean flux because restoring is achieved though negative particulate remieralization (i.e., production)--NOT FLUX
!kst dum_ocn_flux/locijk_focn is set = 0.0 so that it is not included twice when ocn is updated in b_main, only the remin part is included in the update
        dum_ocn_flux(io_PO4,k) = c0
        dum_ocn_flux(io_NO3,k) = c0
       else                          !restore to PO4 field if not prev. state
          if (dum_docn_restore(io_PO4,k) < -const_real_nullsmall) then
             po4_MM = -dum_docn_restore(io_PO4,k)
          else
             po4_MM = 0.0
          endif
       endif
    CASE ('1N1T_PO4MM','2N1T_PO4MM_SiO2')
       ! 1 x nutrient, 1 x 'taxa': PO4 Michaelis-Menton; Si not cosidered in org-c production
       if (ocn(io_PO4,dum_i,dum_j,k) > const_real_nullsmall) then 
          loc_biomass = ocn(io_PO4,dum_i,dum_j,k)
       end if 
    CASE ('3N1T_PNCMM','4N1T_PNCMM_SiO2') 
      ! 3 x nutrient, 1 x 'taxa': P/N/C colimitation Michaelis-Menton
       ! Chikamoto 05/12/2006  Michealis-Menten kinetics [PO4],[NO3] and [CO2]         
       !NOTE:  4N1T is really = 3N1T, since SiO2 is not really considered here....SiO2 "limitation" is via bio_part_red(is_POC,is_opal,:,:,:) calc'ed above

       if (ocn(io_PO4,dum_i,dum_j,k) > const_real_nullsmall .and. & 
            &     ocn(io_NO3,dum_i,dum_j,k) > const_real_nullsmall .and. & 
            &     carb(ic_conc_co2,dum_i,dum_j,k) > const_real_nullsmall ) then

          loc_biomass = min(po4_lim*po4_MM,no3_lim*no3_MM,co2_lim*co2_MM) 
          
          if ( abs( loc_biomass - po4_lim*po4_MM ) <= const_real_nullsmall ) MM_index(dum_i,dum_j,k) = 1.  
          if ( abs( loc_biomass - no3_lim*no3_MM ) <= const_real_nullsmall ) MM_index(dum_i,dum_j,k) = 2.
          if ( abs( loc_biomass - co2_lim*co2_MM ) <= const_real_nullsmall ) MM_index(dum_i,dum_j,k) = 3.
          
          po4_MM = c1

      endif

  CASE ('5N1T_PNCFeMM_SiO2' ) ! 
      ! 5 x nutrient, 1 x 'taxa': P/N/C colimitation Michaelis-Menton
       ! Chikamoto 05/12/2006  Michealis-Menten kinetics [PO4],[NO3] and [CO2], [SiO2], [Fe]
!                loc_FeT = ocn(io_Fe,dum_i,dum_j,n_kmax) + ocn(io_FeL,dum_i,dum_j,n_kmax)                        !already done above, besides, uses surface instead of k
!  see note for 4N1T re: SiO2
       if (ocn(io_PO4,dum_i,dum_j,k) > const_real_nullsmall .and. & 
            &     ocn(io_NO3,dum_i,dum_j,k) > const_real_nullsmall .and. & 
            &     carb(ic_conc_co2,dum_i,dum_j,k) > const_real_nullsmall .and. & 
            &     loc_FeT > const_real_nullsmall ) then 
          
          ! nutrient_MM is ratio, nutrient_lim is concentration[mol/kg] 
          
          fe_MM = loc_FeT/(par_bio_c0_Fe + loc_FeT)                  !one nutr.1/2 sat for all phytoplankton
          
!          loc_junk2 = min(po4_MM,no3_MM,co2_MM,fe_MM)       !limits by uptake rate
          
          fe_lim = (loc_FeT)/par_bio_red_POP_POC/bio_part_red(is_POC,is_POFe,dum_i,dum_j)
          
!          loc_junk1 = min(po4_lim,no3_lim,co2_lim,fe_lim)   !limits by concentration in P units

          loc_biomass = min(po4_lim*po4_MM,no3_lim*no3_MM,fe_lim*fe_MM,co2_lim*co2_MM)

          if ( abs( loc_biomass - po4_lim*po4_MM ) <= const_real_nullsmall ) MM_index(dum_i,dum_j,k) = 1.  
          if ( abs( loc_biomass - no3_lim*no3_MM ) <= const_real_nullsmall ) MM_index(dum_i,dum_j,k) = 2.
          if ( abs( loc_biomass - co2_lim*co2_MM ) <= const_real_nullsmall ) MM_index(dum_i,dum_j,k) = 3.
          if ( abs( loc_biomass - fe_lim*fe_MM   ) <= const_real_nullsmall ) MM_index(dum_i,dum_j,k) = 4.
          po4_MM = c1

       endif
       
       
    CASE ('5N2T_PNCFeMM_SiO2' ) ! 
       ! 5 x nutrient, 2 x 'taxa': P/N/C colimitation Michaelis-Menton
       ! Chikamoto 05/12/2006  Michealis-Menten kinetics [PO4],[NO3] and [CO2], [SiO2], [Fe]
       
       if (ocn(io_PO4,dum_i,dum_j,k) > const_real_nullsmall .and. & 
            &     ocn(io_NO3,dum_i,dum_j,k) > const_real_nullsmall .and. & 
            &     carb(ic_conc_co2,dum_i,dum_j,k) > const_real_nullsmall .and. & 
            &     loc_FeT > const_real_nullsmall ) then 
          
          
          ! nutrient_MM is ratio (rate), nutrient_lim is concentration: [mol/kg] AH/KT 7-16-09
 ! OK, now we have half sats varying with phytoplankton class:
         
          po4_MM_lg = ocn(io_PO4,dum_i,dum_j,k)/(par_bio_c0_PO4_lg + ocn(io_PO4,dum_i,dum_j,k))          
          po4_MM_sm = ocn(io_PO4,dum_i,dum_j,k)/(par_bio_c0_PO4_sm + ocn(io_PO4,dum_i,dum_j,k))
          
          no3_MM_lg = ocn(io_NO3,dum_i,dum_j,k)/(par_bio_c0_NO3_lg + ocn(io_NO3,dum_i,dum_j,k))
          no3_MM_sm = ocn(io_NO3,dum_i,dum_j,k)/(par_bio_c0_NO3_sm + ocn(io_NO3,dum_i,dum_j,k))
       
          co2_MM_lg = 10*carb(ic_conc_co2,dum_i,dum_j,k)/(par_bio_c0_CO2_lg + carb(ic_conc_co2,dum_i,dum_j,k)) 
          co2_MM_sm = 10*carb(ic_conc_co2,dum_i,dum_j,k)/(par_bio_c0_CO2_sm + carb(ic_conc_co2,dum_i,dum_j,k)) 
       
          fe_MM_lg = loc_FeT/(par_bio_c0_Fe_lg + loc_FeT)            !Hannah
          fe_MM_sm = loc_FeT/(par_bio_c0_Fe_sm + loc_FeT)

          if ( ocn(io_SiO2,dum_i,dum_j,k)  > const_real_nullsmall ) then
             sio2_MM = ocn(io_SiO2,dum_i,dum_j,k)/(par_bio_c0_SiO2 + ocn(io_SiO2,dum_i,dum_j,k) )
          else
             sio2_MM = c0
          endif
          loc_rate_lim_lg = min(po4_MM_lg,no3_MM_lg,co2_MM_lg,fe_MM_lg,sio2_MM)    
          loc_rate_lim_sm = min(po4_MM_sm,no3_MM_sm,co2_MM_sm,fe_MM_sm)    

!**********Hannah 2 July 2010****from above, L476********2 Taxa
          if (sed_select(is_POFe)) then
             
             !  for now:  hardwire in, then make these adjustable and read in biobio_5:
             par_bio_FetoC_C_lg = 15000.0
             par_bio_FetoC_C_sm = 20000.0
             par_bio_FetoC_K_lg = 115623.0
             par_bio_FetoC_K_sm = 31805.0
             par_bio_FetoC_pP = -0.65
             
             !NOTE:  for the following parameters, the POC/POFE ratio will not exceed ~340 000 regardless of the max in bio_bio_5, because
             !       of the 'usual' surface [FeT] - it doesn't reach less than about 
             !       par_bio_FetoC_C  = 0.
             !       par_bio_FetoC_K = 94500.
             !       par_bio_FetoC_K = 103684.
             !       par_bio_FetoC_pP = -0.4225
             if (.NOT. opt_bio(iopt_bio_Fe_fixedFetoC) ) then   !usually not fixed. ie: this if is done
                if (loc_FeT > par_part_red_FeTmin) then       
                   loc_bio_part_red_POC_POFe_lg = 1.0/ &
                        & MIN(par_part_red_FetoCmax, (par_bio_FetoC_C_lg + par_bio_FetoC_K_lg* (1.0E9*loc_FeT)**par_bio_FetoC_pP))    
                   loc_bio_part_red_POC_POFe_sm = 1.0/ &
                        & MIN(par_part_red_FetoCmax, (par_bio_FetoC_C_sm + par_bio_FetoC_K_sm* (1.0E9*loc_FeT)**par_bio_FetoC_pP))            
                else
                   loc_bio_part_red_POC_POFe_lg = 1.0/par_part_red_FetoCmax
                   loc_bio_part_red_POC_POFe_sm = 1.0/par_part_red_FetoCmax
                endif
             else
                loc_bio_part_red_POC_POFe_lg = bio_part_red(is_POC,is_POFe,dum_i,dum_j)
                loc_bio_part_red_POC_POFe_sm = bio_part_red(is_POC,is_POFe,dum_i,dum_j)
             end if
          end if

          !only fe has differnet lg and sm fe_conc_lim becuase of the C:Fe can vary
          fe_lim_lg = (loc_FeT)/par_bio_red_POP_POC/loc_bio_part_red_POC_POFe_lg 
          fe_lim_sm = (loc_FeT)/par_bio_red_POP_POC/loc_bio_part_red_POC_POFe_sm 
          
          
          !++++++++Hannah 2 July 2010 - Diatom biomass, [Fe] limitation, Si:N uptake ratio ++++++++++++++++variable Si:N available only in 5N2T
          !for now:  hardwire [Fe] and fit coefficients:
          loc_fet_max = 6.0e-10
          loc_fet_min = 1.0e-11 
          loc_a1 = -9.0909090909E9
          loc_a2 = 6.454545454545
          loc_a3 = 2e-7
          loc_a4 = -0.72
          loc_max_si2n = 5.2

          if ( ocn(io_SiO2,dum_i,dum_j,k) > const_real_nullsmall ) then
#ifdef powsi2n
             if (loc_FeT > loc_fet_min) then
                bio_part_red(is_PON,is_opal,dum_i,dum_j) = loc_a3*((loc_FeT)**loc_a4)
             else
                bio_part_red(is_PON,is_opal,dum_i,dum_j) = loc_max_si2n             !thick shells if Fe depleted--OK:  this bio_part_red(is_PON,is_opal,:,:) 
             end if

!             ! km 1/2012 - try to calibrate Si/N...but FeT needs to be more depleted perhaps
!             loc_a1 = ((loc_FeT*1e9)**(-1))/5.0 + 0.667
!             if (loc_a1 > 6.0) then
!                bio_part_red(is_PON,is_opal,dum_i,dum_j) = 6.0          ! FeT deplete, so thick shell
!             elseif (loc_a1 < 1.0) then
!                bio_part_red(is_PON,is_opal,dum_i,dum_j) = 1.0          ! FeT replete, so thin shell
!             else
!                bio_part_red(is_PON,is_opal,dum_i,dum_j) = loc_a1       ! power law dependence
!             end if

#else
!km 5/2012             if (loc_FeT < loc_fet_max) then
!km 5/2012                if (loc_FeT > loc_fet_min) then
!km 5/2012                   bio_part_red(is_PON,is_opal,dum_i,dum_j) = loc_a1*(loc_FeT) + loc_a2
!km 5/2012                else
!km 5/2012                   bio_part_red(is_PON,is_opal,dum_i,dum_j) = 6.0         !thick shells if Fe depleted--OK:  this bio_part_red(is_PON,is_opal,:,:) is only what will
!km 5/2012                endif                                                     !   go to POC, so  opal/PON(total) =opal/PON(POC)*(1-DOMfrac)
!km 5/2012             else
!km 5/2012                bio_part_red(is_PON,is_opal,dum_i,dum_j) = 1.0             !thin shells if plenty fe (FeT >= loc_fet_max)
!km 5/2012             endif

             ! km inverse formulation
             bio_part_red(is_PON,is_opal,dum_i,dum_j) = 2.0e-10/loc_FeT
             if ( bio_part_red(is_PON,is_opal,dum_i,dum_j)<0.3 ) then
                bio_part_red(is_PON,is_opal,dum_i,dum_j) = 0.3
             end if
!            ! km steep curve
!             bio_part_red(is_PON,is_opal,dum_i,dum_j) = (1.0e-10/loc_FeT)**2
!             if ( bio_part_red(is_PON,is_opal,dum_i,dum_j)<0.3 ) then
!                bio_part_red(is_PON,is_opal,dum_i,dum_j) = 0.3
!             else if ( bio_part_red(is_PON,is_opal,dum_i,dum_j)>5 ) then
!                bio_part_red(is_PON,is_opal,dum_i,dum_j) = 5.0
!             end if
#endif
#ifdef SitoNconst
             !seasonal Si:N mask  8/10/10
             if (biostep_test == 0 ) then
                print*,'error on seasonal previous state read'
                print*,'biostep_stop = ',biostep_test
                stop
             endif
             bio_part_red(is_PON,is_opal,dum_i,dum_j) = SitoN_season1(dum_i,dum_j,k,biostep)
#endif
             bio_part_red(is_POC,is_opal,dum_i,dum_j) = bio_part_red(is_PON,is_opal,dum_i,dum_j)* & 
                  & bio_part_red(is_POC,is_PON,dum_i,dum_j)*(1.0 - par_bio_red_DOMfrac)             
             sio2_lim = ocn(io_SiO2,dum_i,dum_j,k)/bio_part_red(is_POC,is_opal,dum_i,dum_j)/par_bio_red_POP_POC
          else
             bio_part_red(is_POC,is_opal,dum_i,dum_j) = 0.0  !if no SiO2, no opal produciton
             sio2_lim = 0.0
          endif

          !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
          bio_part_red_PON_opal(dum_i,dum_j,k) = bio_part_red(is_PON,is_opal,dum_i,dum_j)
!po4_lim = [PO4], no3_lim=[NO3]in PO4 units
          loc_limiter_lg = min(po4_lim*po4_MM_lg,no3_lim*no3_MM_lg,co2_lim*co2_MM_lg,fe_lim_lg*fe_MM_lg,sio2_lim*sio2_MM)
          loc_limiter_sm = min(po4_lim*po4_MM_sm,no3_lim*no3_MM_sm,co2_lim*co2_MM_sm,fe_lim_sm*fe_MM_sm)

#ifdef rclim  !lim(rate)*lim(conc)  , not necessarily of the same nutrient...
          loc_conc_lim_lg = min(po4_lim,no3_lim,co2_lim,fe_lim_lg,sio2_lim)
          loc_conc_lim_sm = min(po4_lim,no3_lim,co2_lim,fe_lim_sm)

          loc_biomass_lg = (loc_rate_lim_lg/(loc_rate_lim_lg + loc_rate_lim_sm))*loc_rate_lim_lg*loc_conc_lim_lg !ensure you don't use up more po4 than there is!
          loc_biomass_sm = (loc_rate_lim_sm/(loc_rate_lim_lg + loc_rate_lim_sm))*loc_rate_lim_sm*loc_conc_lim_sm !ensure you don't use up more po4 than there is!

          if ( abs( loc_rate_lim_lg - po4_MM_lg ) <= const_real_nullsmall ) then                 
             MM_index_lg(dum_i,dum_j,k) = 1.                                                                          !  1 = PO4 limits
          elseif ( abs(loc_rate_lim_lg - no3_MM_lg) <= const_real_nullsmall )  then                    
             MM_index_lg(dum_i,dum_j,k) = 2.   
          elseif ( abs(loc_rate_lim_lg - co2_MM_lg) <= const_real_nullsmall )  then                    
             MM_index_lg(dum_i,dum_j,k) = 3.   
          elseif ( abs(loc_rate_lim_lg - fe_MM_lg) <= const_real_nullsmall )  then                      
             MM_index_lg(dum_i,dum_j,k) = 4.   
          elseif ( abs(loc_rate_lim_lg - sio2_MM) <= const_real_nullsmall )  then                      
             MM_index_lg(dum_i,dum_j,k) = 5.   
          endif

          if ( abs( loc_rate_lim_sm - po4_MM_sm ) <= const_real_nullsmall ) then                 
             MM_index_sm(dum_i,dum_j,k) = 1.                                                                          !  1 = PO4 limits
          elseif ( abs(loc_rate_lim_sm - no3_MM_sm) <= const_real_nullsmall )  then                    
             MM_index_sm(dum_i,dum_j,k) = 2.   
          elseif ( abs(loc_rate_lim_sm - co2_MM_sm) <= const_real_nullsmall )  then                    
             MM_index_sm(dum_i,dum_j,k) = 3.   
          elseif ( abs(loc_rate_lim_sm - fe_MM_sm) <= const_real_nullsmall )  then                      
             MM_index_sm(dum_i,dum_j,k) = 4.   
          endif
#else  !limiter is min(rate*conc), as opposed to originally, which looked only at min(rate) and let biomass = min(rateNut)*concNut, which could result in concNut < 0...
!    this is default by before apr, 2011
!km 3/2012 MM_INDEX could give erroneous result in fields_biogem_3d.nc 
!km        e.g., N-limitation (=2) 6 months out of a year, and Fe limitation (=4) during the remaining 6 months
!km        then on annual mean, average of 2 and 4 gives 3 (CO2 limitation!)  This in fact occurs as the model does not 
!km        print the co2_sm limitation comments below in runlog.txt but shows the co2_sm limitation in 3d.nc
!km        Presumably this is less likely to happen when 3d output is monthly mean instead of annual mean.

          if ( abs( loc_limiter_lg - po4_lim*po4_MM_lg ) <= const_real_nullsmall ) then                 
             loc_biomass_lg = (loc_rate_lim_lg/(loc_rate_lim_lg + loc_rate_lim_sm))*(po4_lim)*po4_MM_lg 
             MM_index_lg(dum_i,dum_j,k) = 1.                                                                          !  1 = PO4 limits
          elseif ( abs(loc_limiter_lg - no3_lim*no3_MM_lg) <= const_real_nullsmall )  then                    
             loc_biomass_lg = (loc_rate_lim_lg/(loc_rate_lim_lg + loc_rate_lim_sm))*(no3_lim)*no3_MM_lg
             MM_index_lg(dum_i,dum_j,k) = 2.   
          elseif ( abs(loc_limiter_lg - co2_lim*co2_MM_lg) <= const_real_nullsmall )  then                    
             loc_biomass_lg = (loc_rate_lim_lg/(loc_rate_lim_lg + loc_rate_lim_sm))*(co2_lim)*co2_MM_lg
             MM_index_lg(dum_i,dum_j,k) = 3.   
          elseif ( abs(loc_limiter_lg - fe_lim_lg*fe_MM_lg) <= const_real_nullsmall )  then                      
             loc_biomass_lg = (loc_rate_lim_lg/(loc_rate_lim_lg + loc_rate_lim_sm))*(fe_lim_lg)*fe_MM_lg
             MM_index_lg(dum_i,dum_j,k) = 4.   
          elseif ( abs(loc_limiter_lg - sio2_lim*sio2_MM) <= const_real_nullsmall )  then                      
             loc_biomass_lg = (loc_rate_lim_lg/(loc_rate_lim_lg + loc_rate_lim_sm))*(sio2_lim)*sio2_MM
             MM_index_lg(dum_i,dum_j,k) = 5.   
          endif

          if ( abs( loc_limiter_sm - po4_lim*po4_MM_sm ) <= const_real_nullsmall ) then                 
             loc_biomass_sm = (loc_rate_lim_sm/(loc_rate_lim_lg + loc_rate_lim_sm))*(po4_lim)*po4_MM_sm
             MM_index_sm(dum_i,dum_j,k) = 1.    
          elseif ( abs(loc_limiter_sm - no3_lim*no3_MM_sm) <= const_real_nullsmall )  then                    
             loc_biomass_sm = (loc_rate_lim_sm/(loc_rate_lim_lg + loc_rate_lim_sm))*(no3_lim)*no3_MM_sm
             MM_index_sm(dum_i,dum_j,k) = 2.
          elseif ( abs(loc_limiter_sm - co2_lim*co2_MM_sm) <= const_real_nullsmall )  then                    
             loc_biomass_sm = (loc_rate_lim_sm/(loc_rate_lim_lg + loc_rate_lim_sm))*(co2_lim)*co2_MM_sm
             MM_index_sm(dum_i,dum_j,k) = 3.

!km             print*,'km-co2 sm limitation: ',loc_limiter_sm,co2_lim*co2_MM_sm, co2_lim, co2_MM_sm,co2_MM_sm
!km             print*,' p:  ',po4_lim*po4_MM_sm, po4_lim, po4_MM_sm
!km             print*,' n:  ',no3_lim*no3_MM_sm, no3_lim, no3_MM_sm
!km             print*,' fe: ',fe_lim_sm*fe_MM_sm, fe_lim_sm, fe_MM_sm

          elseif ( abs(loc_limiter_sm - fe_lim_sm*fe_MM_sm) <= const_real_nullsmall )  then                      
             loc_biomass_sm = (loc_rate_lim_sm/(loc_rate_lim_lg + loc_rate_lim_sm))*(fe_lim_sm)*fe_MM_sm
             MM_index_sm(dum_i,dum_j,k) = 4.
          endif
#endif
          po4_MM = c1
       endif
       
       !***************************************************
       
       par_bio_k0_PO4_lg = 1 / nuts_tau_lg * loc_temp * loc_mld * loc_biomass_lg         !nuts_tau = optimum nutrient uptake timescale in yr, loc_biomass=mol kg-1
       par_bio_k0_PO4_sm = 1 / nuts_tau_sm * loc_temp * loc_mld * loc_biomass_sm         !nuts_tau = optimum nutrient uptake timescale in yr, loc_biomass=mol kg-1
       
       
       loc_dPO4_lg = dtime* &                 !time step length (yr)
            & loc_ficefree* &               !km ice fraction
            & loc_kI* &                     !km light dependence (fraction)
            & po4_MM* &                     !km nutrient uptake in terms of PO4 (fraction)
            & par_bio_k0_PO4_lg             !km max growth rate   (mol kg-1 yr-1)
       loc_dPO4_sm = dtime* &                 !time step length (yr)
            & loc_ficefree* &               !km ice fraction
            & loc_kI* &                     !km light dependence (fraction)
            & po4_MM* &                     !km nutrient uptake in terms of PO4 (fraction)
            & par_bio_k0_PO4_sm             !km max growth rate   (mol kg-1 yr-1)
       loc_dPO4 = loc_dPO4_lg + loc_dPO4_sm
       
       !Hannah 6 July 2010  From Trunk
       if (sed_select(is_POFe)) then
          if (.NOT.opt_bio(iopt_bio_Fe_fixedFetoC)) then
             if (loc_dPO4 > const_real_nullsmall) then       !get a weighted average POFe:POC ratio for lg+sm 
                bio_part_red(is_POC,is_POFe,dum_i,dum_j) = (loc_dPO4_lg*loc_bio_part_red_POC_POFe_lg + loc_dPO4_sm*loc_bio_part_red_POC_POFe_sm)/loc_dPO4
                bio_part_red(is_POFe,is_POC,dum_i,dum_j) = 1.0/bio_part_red(is_POC,is_POFe,dum_i,dum_j)
             else
                bio_part_red(is_POC,is_POFe,dum_i,dum_j) = 0.
                bio_part_red(is_POFe,is_POC,dum_i,dum_j) = 0.
             endif
          end if
       end if
       
       bio_part_lg(dum_i,dum_j,k) = bio_part_red(is_POP,is_POC,dum_i,dum_j)*loc_dPO4_lg !note:this is POC only (not DOM frac like bio_part(is_POC)
       bio_part_sm(dum_i,dum_j,k) = bio_part_red(is_POP,is_POC,dum_i,dum_j)*loc_dPO4_sm
       
    CASE default
    end select

    SELECT CASE (par_bio_prodopt)
    CASE ('5N2T_PNCFeMM_SiO2' ) ! loc_dPO4 had to be calculated above, but for all else, see default:
    CASE default       
       if (sed_select(is_opal)) then
          ! allow opal production indedpendent of organic carbon production
          ! modify Si:C uptake according to H4SiO4 depletion
          ! NOTE: this is just to ensure that SiO2 does not drop below zero,                                          variable Si:N not available in anything but 5N2T
          !       and does not consistute a true 'nutrient limitation' of productivity
          if (ocn(io_SiO2,dum_i,dum_j,k) > const_real_nullsmall) then
             bio_part_red(is_POC,is_opal,dum_i,dum_j) = (1.0 - par_bio_red_DOMfrac)*par_bio_red_POC_opal*  &
                  & ocn(io_SiO2,dum_i,dum_j,k)/(par_bio_c0_SiO2 + ocn(io_SiO2,dum_i,dum_j,k))
          else
             bio_part_red(is_POC,is_opal,dum_i,dum_j) = 0.0
          end if
       end if
       par_bio_k0_PO4 = 1 / nuts_tau * loc_temp * loc_mld * loc_biomass         !nuts_tau = optimum nutrient uptake timescale in yr, loc_biomass=mol kg-1
       
       loc_dPO4 = dtime* &                      !time step length (yr)
            & loc_ficefree* &               !km ice fraction
            & loc_kI* &                     !km light dependence (fraction)
            & po4_MM* &                     !km nutrient uptake in terms of PO4 (fraction)
            & par_bio_k0_PO4                !km max growth rate   (mol kg-1 yr-1)
    end select       

    
#ifdef prodtotalconst
    !seasonal constant 5/12/09:
    if (biostep_test == 0 ) then
       print*,'error on seasonal previous state read'
       print*,'biostep_stop = ',biostep_test
       stop
    endif
    loc_dPO4 = dPO4_season1(dum_i,dum_j,k,biostep)
#endif
    dPO4_season(dum_i,dum_j,k,biostep) = loc_dPO4       !save seasonal dPO4 for restart
    
    ! *** CALCULATE ORGANIC CARBON PRODUCTION ***
    ! calculate export in currency of particulate carbon (rather than PO4)
    ! NOTE: put everything into particulate form initially, but re-scale later to account for DOM export


    bio_part(is_POC,dum_i,dum_j,k) = bio_part_red(is_POP,is_POC,dum_i,dum_j)*loc_dPO4

    ! *** CALCULATE ASSOCIATED BULK EXPORT ***
    ! >>> NON-GENERIC ALGORITHM (will be generic later)
    ! calculate any associated 'Redfield' and isotopic tracer components
    ! NOTE: treat <is_POM_recfrac> as a special case, and set to default value (as set in <biogem_config.par>)
    ! NOTE: sed_select(is) == 1,2 are the bulk biogenic particulates (POC,PON,POP,POFe,PO2,CaCO3,opal)
    ! NOTE: sed_select(is) == 11,12,13,14 are the isotopic particulates (POC_13C,POC_14C,PON_15N,CaCO3_13C,CaCO3_14C,CaCO3_18O)
    !       that can be related directly to the relevant bulk particulate by the defined dependency == sed_dep(is)
    !       (see gem_config_sed.par)
    DO l=1,n_ismax
       is = conv_iselected_is(l)
       select case (sed_type(is))
       case (par_sed_type_bio,par_sed_type_POM,par_sed_type_CaCO3,par_sed_type_opal) !par_sed_type_bio=1(POC,CaCO3,opal: bulk),par_sed_type_POM=3(POP,PON,POFe),
          SELECT CASE (par_bio_prodopt)                                              !!par_sed_type_CaCO3=4(CaCO3_Cd),par_sed_type_opal=5(opal_Ge ?!),(not scavenged)
          CASE ('5N2T_PNCFeMM_SiO2' ) !                    !NOTE:  radio isotopes haven't been dealt with yet for 2 taxa...
             if (is == is_opal ) then
                bio_part(is,dum_i,dum_j,k) = &                                             
                     & bio_part_red(is_POC,is,dum_i,dum_j)*bio_part_lg(dum_i,dum_j,k)  !(1-DOMfrac already accounted for in opal:POC var from opal:PON)
             elseif( is == is_CaCO3 ) then
                bio_part(is,dum_i,dum_j,k) = &                                             
                     & bio_part_red(is_POC,is,dum_i,dum_j)*bio_part_sm(dum_i,dum_j,k) !(1-DOMfrac already accounted for in CaCO3:POC)
             else
                bio_part(is,dum_i,dum_j,k) = &                                             
                     & bio_part_red(is_POC,is,dum_i,dum_j)*bio_part(is_POC,dum_i,dum_j,k)!(1-DOMfrac)will be adjusted for below,altho is_POC,is_opal isalreadyaccounted for
             endif
          CASE default
             bio_part(is,dum_i,dum_j,k) = &                                             
                  & bio_part_red(is_POC,is,dum_i,dum_j)*bio_part(is_POC,dum_i,dum_j,k)
          END SELECT
       case (11:20)                                                               !13C,14C 
          bio_part(is,dum_i,dum_j,k) = &
               & bio_part_red(sed_dep(is),is,dum_i,dum_j)*bio_part(sed_dep(is),dum_i,dum_j,k)
       end select
       
    end do
#ifdef CaCO3const
    !seasonal constant 5/12/09:
    if (biostep_test == 0 ) then
       print*,'error on seasonal previous state read'
       print*,'biostep_stop = ',biostep_test
       stop
    endif
    bio_part(is_CaCO3,dum_i,dum_j,k) = CaCO3_season1(dum_i,dum_j,k,biostep)
#endif
    CaCO3_season(dum_i,dum_j,k,biostep) = bio_part(is_CaCO3,dum_i,dum_j,k)
   
    ! \/\/\/ THIS CODE TO BE EVENTUALLY REPLACED BY GENERIC ARRAY \/\/\/            
    bio_part(is_POC_frac2,dum_i,dum_j,k)   = par_bio_remin_POC_frac2                      !biogem_bio_   .05455etc
    bio_part(is_CaCO3_frac2,dum_i,dum_j,k) = par_bio_remin_CaCO3_frac2                    !  "      "    .48829etc
    bio_part(is_opal_frac2,dum_i,dum_j,k)  = par_bio_remin_opal_frac2                     !biogem_bio_4  0.0, or 1.0 if not in biobio file...
    ! /\/\/\ THIS CODE TO BE EVENTUALLY REPLACED BY GENERIC ARRAY /\/\/\
    ! convert particulate sediment tracer indexed array concentrations to (dissolved) tracer indexed array
    ! >>> GENERIC ALGORITHM BUT CONTAINS USER-MODIFIABLE INFORMATION IN ARRAY <conv_sed_ocn>
!zB:  
!POC->DIC,O2
!PON ->NO3,O2
!POP->PO4
!POCd->Cd
!POFe->Fe
!POM_Fe->Fe, altho POM_Fe = 0.0 
!CaCO3->DIC
!CaCO3->ALK
!CaCO3->Ca
!opal->SiO2
    DO l=1,n_ismax
       is = conv_iselected_is(l)
       loc_tot_i = conv_sed_ocn_i(0,is)
       do loc_i=1,loc_tot_i
          io = conv_sed_ocn_i(loc_i,is)
          loc_bio_uptake(io) = loc_bio_uptake(io) + conv_sed_ocn(io,is)*bio_part(is,dum_i,dum_j,k)
       end do
    end DO

    ! *** RE-SCALE FOR DISSOLVED ORGANIC MATTER PRODUCTION ***
    ! >>> GENERIC ALGORITHM BUT CONTAINS USER-MODIFIABLE INFORMATION IN ARRAY <conv_POM_DOM>
    DO l=1,n_ismax
       is = conv_iselected_is(l)
       loc_tot_i = conv_POM_DOM_i(0,is)
       do loc_i=1,loc_tot_i
          io = conv_POM_DOM_i(loc_i,is)
 
          select case(sed_type(is)) 
          case(14) 
             ! by M.Chikamoto 07-13-2006 adding 15N fraction between PON_15N and DOM_N_15N 
             ! fraction factor 0 - -2 permil (Miyake and Wada 1971; Records of Oceanographic Works in Japan,11,1-6) 
             !                           PON_15N 0 permil or +1 
             !                           DON_15N 0 permil or -1 (Yoshikawa et al. 2005) 
             ! We assume no fractionation through nitrification without using ammonia.==> loc_delta = 0 
             loc_r15N = 0.0 
             if(bio_part(is_PON,dum_i,dum_j,k) > const_real_nullsmall )then 
                loc_r15N = bio_part(is_PON_15N,dum_i,dum_j,k)/bio_part(is_PON,dum_i,dum_j,k)  
             endif 
             loc_delta = 0. !-1.  
             loc_alpha = 1.0 + loc_delta/1000. 
 
             bio_remin(io,dum_i,dum_j,k) = bio_remin(io,dum_i,dum_j,k) + & 
                  & par_bio_red_DOMfrac*bio_part(sed_dep(is),dum_i,dum_j,k)*loc_r15N*loc_alpha 
 
             bio_part(is,dum_i,dum_j,k) = & 
                  & (1.0-par_bio_red_DOMfrac)*bio_part(sed_dep(is),dum_i,dum_j,k)*loc_r15N*loc_alpha 

    !         zB:     for io= 15,16,17,18,20,21,22 (DOM_C,DOM_C_13C,COM_C_14C,DOM_N,DOM_P,DOM_Cd,DOM_Fe) only POM/DOM_tracers opal and CaCO3 are NOT done here.
    !                and  is=  3, 4, 5, 6, 8, 9,10 (POC,    POC_13C,  POC_14C,  PON,  POP,  POCd, POFe)             they are taken care of elsewhere
          case(1:13,16)
             if (.not. restore_prev_state ) then        !kst added 4/20 
                bio_remin(io,dum_i,dum_j,k) = bio_remin(io,dum_i,dum_j,k) + & 
                     & par_bio_red_DOMfrac*bio_part(is,dum_i,dum_j,k) 
                bio_part(is,dum_i,dum_j,k) = & 
                     & (1.0 - par_bio_red_DOMfrac)*bio_part(is,dum_i,dum_j,k) 
             endif
          end select 

       end do
    end do
    SELECT CASE (par_bio_prodopt)                                             
    CASE ('5N2T_PNCFeMM_SiO2' ) !                    !NOTE bio_part_lg,sm is just for output bookkeeping now.
       bio_part_lg(dum_i,dum_j,k) = & 
            & (1.0 - par_bio_red_DOMfrac)*bio_part_lg(dum_i,dum_j,k) 
       bio_part_sm(dum_i,dum_j,k) = & 
            & (1.0 - par_bio_red_DOMfrac)*bio_part_sm(dum_i,dum_j,k) 
    END SELECT

!added Fe scavanging by newly created POM: (this means Fe is scavenged in the surface layer as well.)
!kst 4/7/2010
    ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

    ! *** Fe SCAVENGING ***
    ! calculate scavenging of Fe from water column by newly formed particulates
    if (ocn_select(io_Fe)) then
       if (ocn(io_Fe,dum_i,dum_j,k) > const_real_nullsmall) then
!          if (dum_i == 10 .and. dum_j == 10 ) print*,'3uptake: ',bio_part(is_POM_Fe,dum_i,dum_j,k),k,bio_remin(io_Fe,dum_i,dum_j,k)
             call sub_calc_scav_Fe( &
                  & dum_dt, &
                  & ocn(io_Fe,dum_i,dum_j,k), &
                  & bio_part(:,dum_i,dum_j,k), &
                  & bio_remin(:,dum_i,dum_j,k) &
                  & )
!          if (dum_i == 10 .and. dum_j == 10 ) print*,'4uptake ',bio_part(is_POM_Fe,dum_i,dum_j,k),k,bio_remin(io_Fe,dum_i,dum_j,k)
       end if
    end if

! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++=

    ! *** SET MODIFICATION OF TRACER CONCENTRATIONS ***
    ! >>> GENERIC ALGORITHM
    ! NOTE: depletion of dissolved species as a result of biological productivity is implimented as 'negative' remineralization
    !kst: O2 is negative, ALK can be either (?)

    DO l=3,n_iomax
       io = conv_iselected_io(l)
       bio_remin(io,dum_i,dum_j,k) = bio_remin(io,dum_i,dum_j,k) - loc_bio_uptake(io) 
    end do
    end DO levels


  end SUBROUTINE sub_calc_bio_uptake
!*************from trunk:
!kst  4/21/10

  SUBROUTINE sub_calc_geochem_Fe(dum_i,dum_j,dum_k1,dum_focnFe)
    ! dummy arguments
    INTEGER,INTENT(in)::dum_i,dum_j,dum_k1
    real,dimension(n_kmax),INTENT(in)::dum_focnFe
    ! local variables
    INTEGER::k
    real::loc_Fe,loc_FeL,loc_L
    real::loc_FeT,loc_LT
    real,DIMENSION(2)::loc_roots

    ! *** CALCULATE Fe SPECIATION THROUGHOUT THE WATER COLUMN***
    DO k=n_kmax,dum_k1,-1
       ! initialize variables
       loc_Fe  = ocn(io_Fe,dum_i,dum_j,k) + bio_remin(io_Fe,dum_i,dum_j,k) + dum_focnFe(k)
       loc_FeL = ocn(io_FeL,dum_i,dum_j,k) + bio_remin(io_FeL,dum_i,dum_j,k)
       loc_L   = ocn(io_Ligand,dum_i,dum_j,k) + bio_remin(io_Ligand,dum_i,dum_j,k)
       loc_FeT = loc_FeL + loc_Fe
       loc_LT  = loc_FeL + loc_L

       ! solve Fe speciation equation:
       ! K = FeL / (Fe*L) (e.g. see: Parekth et al. [2005])
       ! => FeL = Fe*L*K
       !    conservation relations:
       !    FeL + Fe = FeT => Fe = FeT - FeL
       !    FeL + L  = LT  => L  = LT  - FeL
       !    substitute:
       !    FeL = (FeT - FeL)*(LT - FeL)*K
       !    => FeL/K = FeT*LT + FeL^2 - LT*FeL - FeT*FeL
       !    => FeL/K = FeL^2 - (LT + FeT)*FeL + FeT*LT
       !    => 1.0*FeL^2 - (LT + FeT + 1.0/K)*FeL + FeT*LT = 0.0
       !       solve as: ax2 + bx + c = 0.0
       !                 where x = FeL
       loc_roots(:) = fun_quad_root(1.0,-(loc_LT + loc_FeT + 1.0/par_K_FeL),loc_FeT*loc_LT)
       ! filter returned roots
       if (maxval(loc_roots(:)) < const_real_nullsmall) then
          IF (opt_misc(iopt_misc_audit)) THEN
!          IF (ctrl_audit) THEN
             CALL sub_report_error( &
                  & 'biogem_box.f90','sub_calc_geochem_Fe', &
                  & 'No REAL root in Fe speciation calculation (or maybe zero ...).'// &
                  & ' / Data: dum_i,dum_j,k,loc_FeL(OLD),loc_Fe(OLD),loc_L(OLD),loc_FeT(OLD),loc_LT(OLD),', &
                  & 'SOD THIS FOR A GAME OF SOLDIERS: calculation abondoned ...', &
                  & (/real(dum_i),real(dum_j),real(k),loc_FeL,loc_Fe,loc_L,loc_FeT,loc_LT/),.false. &
                  & )
             error_stop = .FALSE.
          end IF
       elseif ((minval(loc_roots(:)) > loc_FeT) .AND. (minval(loc_roots(:)) > loc_LT)) then
          IF (opt_misc(iopt_misc_audit)) THEN
!          IF (ctrl_audit) THEN
             CALL sub_report_error( &
                  & 'biogem_box.f90','sub_calc_geochem_Fe', &
                  & 'No solution to Fe speciation calculation possible ... :('// &
                  & ' / Data: dum_i,dum_j,k,loc_FeL(OLD),loc_Fe(OLD),loc_L(OLD),loc_FeT(OLD),loc_LT(OLD),', &
                  & 'SOD THIS FOR A GAME OF SOLDIERS: calculation abondoned ...', &
                  & (/real(dum_i),real(dum_j),real(k),loc_FeL,loc_Fe,loc_L,loc_FeT,loc_LT/),.false. &
                  & )
             error_stop = .FALSE.
          end IF
       else
          if (minval(loc_roots(:)) < const_real_nullsmall) then
             loc_FeL = maxval(loc_roots(:))
          else
             loc_FeL = minval(loc_roots(:))
          end if
          loc_Fe  = loc_FeT - loc_FeL
          loc_L   = loc_LT - loc_FeL
       end if
       ! re-calculate reminerlization arrays to give rise to calculated Fe speciation
       ! NOTE: subtract <dum_focnFe> again because it is added subsequently in the main BIOGEM loop through <locijk_focn>
       bio_remin(io_Fe,dum_i,dum_j,k)  = loc_Fe - ocn(io_Fe,dum_i,dum_j,k) - dum_focnFe(k)
       bio_remin(io_FeL,dum_i,dum_j,k) = loc_FeL - ocn(io_FeL,dum_i,dum_j,k)
       bio_remin(io_Ligand,dum_i,dum_j,k)   = loc_L - ocn(io_Ligand,dum_i,dum_j,k)

    end DO

  end SUBROUTINE sub_calc_geochem_Fe
  ! ****************************************************************************************************************************** !


  ! *** calculate the oxidation H2S ***
  ! <<< CONTAINS NON-GENERIC ALGORITHMS >>>
  SUBROUTINE sub_calc_bio_remin_oxidize_H2S(dum_i,dum_j,dum_k1,dum_dtyr)
    ! dummy arguments
    INTEGER,INTENT(in)::dum_i,dum_j,dum_k1
    real,intent(in)::dum_dtyr
    ! local variables
    integer::l,io,k
    real::loc_potO2cap,loc_r34S
    real::loc_H2S_oxidation_const,loc_H2S_oxidation
    real,dimension(0:n_ocn,n_maxk)::loc_bio_remin

    ! *** INITIALIZE VARIABLES ***
    ! initialize local variables
    ! change units of H2S oxidation constant from mM-2 hr-1 to M-2 yr-1
    ! and convert from O2 consumption units to H2S units (i.e., divide by 2)
    loc_H2S_oxidation_const = 0.5*const_oxidation_coeff_H2S/conv_hr_yr/(conv_mmol_mol)**2
    ! initialize remineralization tracer arrays
    DO l=3,n_iomax
       io = conv_iselected_io(l)
       loc_bio_remin(io,:) = 0.0
    end do
    
    ! *** OXIDIZE H2S ***
    ! look for some H2S and see if it can be instantaneously oxidized (using O2; if there is any!)
    ! H2S + 2O2 -> SO4 + 2H
    DO k=n_kmax,dum_k1,-1
       ! calculate potential oxidation capacity
       loc_potO2cap = ocn(io_O2,dum_i,dum_j,k) + bio_remin(io_O2,dum_i,dum_j,k)
       if ((ocn(io_H2S,dum_i,dum_j,k) > const_real_nullsmall) .AND. (loc_potO2cap > const_real_nullsmall)) then
          ! calculate H2S oxidation, and cap value at H2S concentration if necessary
          loc_H2S_oxidation = dum_dtyr*loc_H2S_oxidation_const*ocn(io_H2S,dum_i,dum_j,k)*loc_potO2cap**2
          If (loc_H2S_oxidation > ocn(io_H2S,dum_i,dum_j,k)) loc_H2S_oxidation = ocn(io_H2S,dum_i,dum_j,k)
          ! calculate isotopic ratio
          loc_r34S = ocn(io_H2S_34S,dum_i,dum_j,k)/ocn(io_H2S,dum_i,dum_j,k)
          if (loc_H2S_oxidation <= 0.5*loc_potO2cap) then
             ! complete H2S oxidation (no S fractionation)
             loc_bio_remin(io_H2S,k)     = -loc_H2S_oxidation
             loc_bio_remin(io_SO4,k)     = loc_H2S_oxidation
             loc_bio_remin(io_O2,k)      = -2.0*loc_H2S_oxidation
             loc_bio_remin(io_H2S_34S,k) = -loc_r34S*loc_H2S_oxidation
             loc_bio_remin(io_SO4_34S,k) = loc_r34S*loc_H2S_oxidation
          else
             ! partial H2S oxidation (=> S isotope Rayleigh fractionation)
             loc_bio_remin(io_H2S,k) = -0.5*loc_potO2cap
             loc_bio_remin(io_SO4,k) = 0.5*loc_potO2cap
             loc_bio_remin(io_O2,k)  = -loc_potO2cap
!in trunk code 4/21/10: kst
!            loc_bio_remin(io_ALK,k) = -loc_potO2cap
             ! \/\/\/ INSERT ALTERNATIVE CODE FOR NON-ZERO S FRACTIONATION \/\/\/
             loc_bio_remin(io_H2S_34S,k) = -loc_r34S*0.5*loc_potO2cap
             loc_bio_remin(io_SO4_34S,k) = loc_r34S*0.5*loc_potO2cap
             ! /\/\/\ INSERT ALTERNATIVE CODE FOR NON-ZERO S FRACTIONATION /\/\/\
          end if
       end if
    end DO

    ! *** WRITE GLOBAL ARRAY DATA ***
    ! write ocean tracer remineralization field (global array)
    DO l=3,n_iomax
       io = conv_iselected_io(l)
       bio_remin(io,dum_i,dum_j,:) = bio_remin(io,dum_i,dum_j,:) + loc_bio_remin(io,:)
    end do
    
  end SUBROUTINE sub_calc_bio_remin_oxidize_H2S


  ! *** calculate the oxidation of CH4 ***
  ! <<< CONTAINS NON-GENERIC ALGORITHMS >>>
  SUBROUTINE sub_calc_bio_remin_oxidize_CH4(dum_i,dum_j,dum_k1)
    ! dummy arguments
    INTEGER,INTENT(in)::dum_i,dum_j,dum_k1
    ! local variables
    integer::l,io,k
    real::loc_potO2cap,loc_r13C,loc_r14C
    real,dimension(0:n_ocn,n_maxk)::loc_bio_remin

    ! *** INITIALIZE VARIABLES ***
    ! initialize local variables
    ! initialize remineralization tracer arrays
    DO l=3,n_iomax
       io = conv_iselected_io(l)
       loc_bio_remin(io,:) = 0.0
    end do
!in trunk 4/21/10 kst:  see below between insert code^^^^
!    loc_NH4_oxidation_const = 0.16666666 !(yr-1; Fennel et al, 2005)
!    loc_NH4_halfsat = 2.0E-05 !(mol kg-1; Fennel et al., 2005)

    ! *** OXIDIZE CH4 ***
    ! look for some CH4 and see if it can be instantaneously oxidized (using O2, NO3, or SO4; if there is any!)
    ! CH4 + O2 -> CO2 + 2H2
    DO k=n_kmax,dum_k1,-1
       ! calculate potential oxidation capacity
       loc_potO2cap = fun_potO2cap(ocn_select(:),ocn(:,dum_i,dum_j,k),bio_remin(:,dum_i,dum_j,k))
       ! \/\/\/ INSERT CODE \/\/\/
!in trunk 4/21/10 kst
!       if ((ocn(io_NH4,dum_i,dum_j,k) > const_real_nullsmall) .AND. (loc_potO2cap > const_real_nullsmall)) then
!          ! calculate NH4 oxidation, and cap value at NH4 concentration if necessary
!          loc_NH4_oxidation = dum_dtyr*loc_NH4_oxidation_const*ocn(io_NH4,dum_i,dum_j,k)*loc_potO2cap / &
!               & (loc_NH4_halfsat + loc_potO2cap)
!          If (loc_NH4_oxidation > ocn(io_NH4,dum_i,dum_j,k)) loc_NH4_oxidation = ocn(io_NH4,dum_i,dum_j,k)
!          ! calculate isotopic ratio
!          loc_r15N = ocn(io_NH4_15N,dum_i,dum_j,k)/ocn(io_NH4,dum_i,dum_j,k)
!          if (loc_NH4_oxidation <= 0.5*loc_potO2cap) then
!             ! complete NH4 oxidation (no N fractionation)
!             loc_bio_remin(io_NH4,k) = -loc_NH4_oxidation
!             loc_bio_remin(io_NO3,k) = loc_NH4_oxidation
!             loc_bio_remin(io_O2,k)  = -2.0*loc_NH4_oxidation
!             loc_bio_remin(io_ALK,k) = loc_bio_remin(io_NH4,k) - loc_bio_remin(io_NO3,k)
!             loc_bio_remin(io_NH4_15N,k) = -loc_r15N*loc_NH4_oxidation
!             loc_bio_remin(io_NO3_15N,k) = loc_r15N*loc_NH4_oxidation
!          else
!             ! partial NH4 oxidation (=> S isotope Rayleigh fractionation)
!             loc_bio_remin(io_NH4,k) = -0.5*loc_potO2cap
!             loc_bio_remin(io_NO3,k) = 0.5*loc_potO2cap
!             loc_bio_remin(io_O2,k)  = -loc_potO2cap
!             loc_bio_remin(io_ALK,k) = loc_bio_remin(io_NH4,k) - loc_bio_remin(io_NO3,k)
!             ! ### INSERT ALTERNATIVE CODE FOR NON-ZERO N FRACTIONATION ########################################################## !
!             loc_bio_remin(io_NH4_15N,k) = -loc_r15N*0.5*loc_potO2cap
!             loc_bio_remin(io_NO3_15N,k) = loc_r15N*0.5*loc_potO2cap
!             ! ################################################################################################################### !
!          end if
!       end if
!    end DO
       !
       ! /\/\/\ INSERT CODE /\/\/\
    end DO

    ! *** WRITE GLOBAL ARRAY DATA ***
    ! write ocean tracer remineralization field (global array)
    DO l=3,n_iomax
       io = conv_iselected_io(l)
       bio_remin(io,dum_i,dum_j,:) = bio_remin(io,dum_i,dum_j,:) + loc_bio_remin(io,:)
    end do
    
  end SUBROUTINE sub_calc_bio_remin_oxidize_CH4


  ! *** calculate water column remineralization of dissolved organic matter ***
  ! <<< CONTAINS NON-GENERIC ALGORITHMS >>>
  SUBROUTINE sub_calc_bio_remin_DOM(dum_i,dum_j,dum_k1,dum_dtyr)
    ! dummy arguments
    INTEGER,INTENT(in)::dum_i,dum_j,dum_k1
    real,intent(in)::dum_dtyr
    ! local variables
    integer::l,io,is,k
    integer::loc_i,loc_tot_i
    real::loc_potO2cap,loc_O2demand,loc_bio_remin_DOMratio
    real,dimension(0:n_ocn,n_maxk)::loc_bio_remin
    real,dimension(0:n_sed,n_maxk)::loc_bio_part
    real::loc_tanomaly, junkst
 
    ! by M.Chikamoto 07-13-2006  
    real::loc_r15N 
    real::loc_alpha, loc_delta 
    real::delta,ktemp,tv2 

    ! *** INITIALIZE VARIABLES ***
    ! initialize local variables
    ! initialize remineralization tracer arrays
    DO l=3,n_iomax
       io = conv_iselected_io(l)
       loc_bio_remin(io,:) = 0.0
    end do
    ! initialize particulate tracer arrays
    DO l=1,n_ismax
       is = conv_iselected_is(l)
       loc_bio_part(is,:) = 0.0
    end do

    ! *** REMINERALIZE DISSOLVED ORGANIC MATTER ***
    ! NOTE: the new algorithm converts the fraction of DOM marked to be remineralized first into POM before applying the
    !       'usual' generic conversion of sed -> ocn tracers, so as to avoid the need for 'special cases'
    !       (such as of the link between DON and ALK, or POP and ALK)
    DO k=n_kmax,dum_k1,-1
       ! calculate potential oxidation capacity
       loc_potO2cap = fun_potO2cap(ocn_select(:),ocn(:,dum_i,dum_j,k),bio_remin(:,dum_i,dum_j,k))
       If (ocn(io_DOM_C,dum_i,dum_j,k) > const_real_nullsmall) then
          ! check that the DOM lifetime is no shorter than the time-step and modify fraction remineralized accordingly
          if (par_bio_remin_DOMlifetime > dum_dtyr) then
             loc_bio_remin_DOMratio = dum_dtyr/par_bio_remin_DOMlifetime
          else
             loc_bio_remin_DOMratio = 1.0
          end if
#ifdef domrft
          ktemp = log(2.0)/10.            ! Eppley (1972) 0.03013
          loc_bio_remin_DOMratio = loc_bio_remin_DOMratio *exp(ktemp*(ocn(io_T,dum_i,dum_j,k)-273.15))
!in trunk 4/21/10 kst, but commented out...
!!$          loc_DOMlifetime = par_bio_remin_DOMlifetime/(0.59*exp(0.0633*(ocn(io_T,dum_i,dum_j,k) - const_zeroC)))
!  in essence, temperature effect = 0.59*exp(.0633(T))*ratio
          if (loc_bio_remin_DOMratio > const_real_one) loc_bio_remin_DOMratio = 1.0
#endif

          ! calculate potential oxidation capacity
          loc_O2demand = conv_sed_ocn(io_O2,is_POC)*conv_DOM_POM(is_POC,io_DOM_C)* &
               & loc_bio_remin_DOMratio*ocn(io_DOM_C,dum_i,dum_j,k)
          ! compare with potential oxygen availability and modify fraction remineralized accordingly
          if ((loc_O2demand > loc_potO2cap) .AND. (loc_O2demand > const_real_nullsmall)) then
             loc_bio_remin_DOMratio = (loc_potO2cap/loc_O2demand)*loc_bio_remin_DOMratio         
          end if
          ! remineralize dissolved organic matter and add released dissolved inorganic tracers to local remin array
          DO l=3,n_iomax
             io = conv_iselected_io(l)
             loc_tot_i = conv_DOM_POM_i(0,io)
             do loc_i=1,loc_tot_i
                is = conv_DOM_POM_i(loc_i,io)

                if(io == io_NO3_15N)then 
                   ! by M.Chikamoto 07-13-2006 adding 15N fractionation between DON_15N and NO3_15N through nitrification 
                   !                           DON_15N +17 permil  
                   !                           NO3_15N -17 permil  
                   !                               (about NH4+ fractionation Sigman et al., 2005, GBC,19) 
                   ! We assume no fractionation through nitrification without using ammonia.==> loc_delta = 0 
                   loc_r15N = 0. 
                   if(ocn(io_DOM_N,dum_i,dum_j,k) > const_real_nullsmall )then 
                      loc_r15N = ocn(io_DOM_N_15N,dum_i,dum_j,k)/ocn(io_DOM_N,dum_i,dum_j,k) 
                   endif 
                   loc_delta = 0. !-17.0 
                   loc_alpha = 1.0 + loc_delta/1000. 
 
                   delta = loc_bio_remin_DOMratio*conv_DOM_POM(is,io)*ocn(io_DOM_N,dum_i,dum_j,k) 
 
                   loc_bio_part(is,k) = loc_alpha * loc_r15N * delta   
 
                   loc_bio_remin(io,k) = -loc_bio_part(is,k)  
                else 
                   loc_bio_part(is,k)  = loc_bio_remin_DOMratio*conv_DOM_POM(is,io)*ocn(io,dum_i,dum_j,k) 
                   loc_bio_remin(io,k) = -loc_bio_part(is,k) 
                endif 
             end do
          end do

          DO l=1,n_ismax
             is = conv_iselected_is(l)
             loc_tot_i = conv_sed_ocn_i(0,is)
             do loc_i=1,loc_tot_i
                io = conv_sed_ocn_i(loc_i,is)
                loc_bio_remin(io,k) = loc_bio_remin(io,k) + conv_sed_ocn(io,is)*loc_bio_part(is,k)
             end do
          end DO
       end if
    end DO

    ! *** WRITE GLOBAL ARRAY DATA ***
    ! write ocean tracer remineralization field (global array)
    DO l=3,n_iomax
       io = conv_iselected_io(l)
       bio_remin(io,dum_i,dum_j,:) = bio_remin(io,dum_i,dum_j,:) + loc_bio_remin(io,:)
    end do

  end SUBROUTINE sub_calc_bio_remin_DOM


  ! *** calculate water column remineralization of pariculate organic matter ***
  ! <<< CONTAINS NON-GENERIC ALGORITHMS >>>
  SUBROUTINE sub_calc_bio_remin(dum_i,dum_j,dum_k1,dum_dtyr)
    ! dummy arguments
    INTEGER,INTENT(in)::dum_i,dum_j,dum_k1
    real,intent(in)::dum_dtyr
    ! local variables
    integer::l,io,is
    INTEGER::k,kk,loc_bio_remin_min_k
    integer::loc_io,loc_i,loc_tot_i
    real::loc_potO2cap,loc_O2demand
    real::loc_T,loc_SiO2                                               ! 
    real::loc_Si_eq,loc_u                                              ! 
    real::loc_bio_remin_dD
    real::loc_bio_remin_max_D
    real::loc_bio_remin_layerratio
    real::loc_bio_remin_sinkingrate                                    ! 
    real::loc_bio_remin_dt                                             ! layer residence time (in years)
    real::loc_bio_remin_POC_frac1,loc_bio_remin_POC_frac2
    real::loc_bio_part_POC_ratio
    real::loc_bio_remin_CaCO3_frac1,loc_bio_remin_CaCO3_frac2
    real::loc_bio_part_CaCO3_ratio
    real::loc_bio_remin_opal_frac1,loc_bio_remin_opal_frac2
    real::loc_bio_part_opal_ratio
    real,dimension(0:n_sed,n_maxk)::loc_bio_part_TMP
    real,dimension(0:n_sed,n_maxk)::loc_bio_part_OLD
    real,dimension(0:n_sed,n_maxk)::loc_bio_part
    real,dimension(0:n_ocn,n_maxk)::loc_bio_remin
    real,dimension(0:n_sed,n_maxk)::loc_bio_settle

    ! Chikamoto by 06-08-2006
!  kst 12/17/08, from 11/01/06   also added newtemp and loc_tc_anom, junkst
    real::ktemp, r0, dcdt, O2ratio, newtemp, loc_tc_anom, junkst,tv1,tv2,tv3,tv4
    real::loc_bio_part_TMPs, loc_bio_part_POC_ratios
    real::delta_PON, delta_NO3, delta_NO3_15N 
    real::loc_delta,loc_alpha 
    real::loc_r15N 
  

    ! *** USER-DEFINABLE OPTIONS ***
    ! NOTE: settings not included in the run-time configuration files for clarity
    ! ******************************
!    par_bio_remin_opal_K = 0.019/conv_d_yr ! opal particulate base dissolution rate (d-1 -> yr-1) [Ridgwell, 2001]

    ! turning for 4 m/day  by M. Chikamoto 06-07-2006
!    par_bio_remin_opal_K = 0.025/conv_d_yr 
!    if ( par_bio_prodopt == '5N2T_PNCFeMM_SiO2' ) par_bio_remin_opal_K = nuts_tau/conv_d_yr    !only used for testing,
    ! ******************************

    ! *** INITIALIZE VARIABLES ***
    ! initialize local variables
    loc_bio_remin_sinkingrate = par_bio_remin_sinkingrate

    ! initialize particulate tracer arrays
    DO l=1,n_ismax
       is = conv_iselected_is(l)
       ! copy particulate tracer field to temporary array and reset value of particulate tracer field
       loc_bio_part_OLD(is,:) = bio_part(is,dum_i,dum_j,:)
       bio_part(is,dum_i,dum_j,:)   = 0.0
       ! initialize local particulate tracer field and settling flux arrays
       loc_bio_part(is,:) = 0.0
       loc_bio_settle(is,:) = 0.0
    end do
    ! initialize remineralization tracer arrays
    DO l=3,n_iomax
       io = conv_iselected_io(l)
       loc_bio_remin(io,:) = 0.0
    end do

    ! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
    ! *** k WATER-COLUMN LOOP START ***
    ! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>


    kloop : DO k=n_kmax,dum_k1,-1      ! 16(surface) : bottom(dum_k1)
       ! find some particulates in the water column
       ! NOTE: assume that there will always be some particulate organic carbon if there is going to be anything at all
       ! km : w/ nlayer_prod=2, the following is True for k=16 and k=15

       If (loc_bio_part_OLD(is_POC,k) > const_real_nullsmall) then
          ! if the identified particulate material is already residing in the bottom-most ocean layer, flag as sediment flux
          If (k == dum_k1) then
             loc_bio_remin_min_k = dum_k1 - 1
          else
             ! determine the deepest layer that sinking material can reach within the current time-step
             ! NOTE: do this regardless of whether a fixed remineralization profile is selected, or whether 
             ! remineralization is calculated as a function of residence time in an ocena layer (and ambient environmental conditions)
             ! NOTE: trap the situation where the depth of the sediment surface is surpassed
             ! NOTE: start loop from the layer lying below the one in which the identified particulate material resides
             ! phys_ocn(ipo_Dbot,i,j,k) = depth of the bottom of layer k 
             loc_bio_remin_min_k = dum_k1 - 1
             loc_bio_remin_max_D = phys_ocn(ipo_Dbot,dum_i,dum_j,k) + dum_dtyr*loc_bio_remin_sinkingrate
             do kk=k-1,dum_k1,-1
                If (phys_ocn(ipo_Dbot,dum_i,dum_j,kk) > loc_bio_remin_max_D) then
                   loc_bio_remin_min_k = kk
                   exit
                end if
             end do
          end if
          
          ! zero local (temporary) particulate field array, and seed value at location in water column identified
          DO l=1,n_ismax
             is = conv_iselected_is(l)
             loc_bio_part_TMP(is,:) = 0.0
             loc_bio_part_TMP(is,k) = loc_bio_part_OLD(is,k)
          end do
          ! >>>>>>>>>>>>>>>>>>>>>>>>>
          ! *** kk SUB-LOOP START ***
          ! >>>>>>>>>>>>>>>>>>>>>>>>>

          ! for each of the three (POC, CaCO3, and opal) primary remineralizable species (if selected),
          ! loop down remineralization column identified previously;
          ! (1) calculating the fractional remineralization in each layer, moving the particulate remainder to the layer below
          ! (2) calculate tracer remineralization from particulate supply from layer above
          ! (3) update particulate tracer field for current layer
          ! then, if the sediments are reached, calculate sediment flux
          ! NOTE: the particulate tracer field is in units of mol kg-1, following the (dissolved) ocean tracers, and as a result,
          !       corrections must be made for changes in ocean layer thickness
          junkst = 0.0
          k2loop : do kk=k-1,loc_bio_remin_min_k,-1
             ! test to see whether the ocean bottom has been reached
             If (kk >= dum_k1) then                
                ! calculate ratio of layer thicknesses
                ! (used to convert particulate matter concentrations as particulates settle through the water column
                !  comprising layers of non-uniform thickness)
                loc_bio_remin_layerratio = phys_ocn(ipo_dD,dum_i,dum_j,kk+1)/phys_ocn(ipo_dD,dum_i,dum_j,kk)     ! layer thickness above/thickness below (ratio)
                loc_bio_remin_dD = phys_ocn(ipo_dD,dum_i,dum_j,kk)                                               !layer thickness(kk) in meters
                ! calculate residence time (yr) of particulates in ocean layer (from layer thickness and sinking speed)
                if (loc_bio_remin_sinkingrate > const_real_nullsmall) loc_bio_remin_dt = loc_bio_remin_dD/loc_bio_remin_sinkingrate     !!residence time in layer kk (years)
                junkst = junkst + loc_bio_remin_dt
                if (junkst > dum_dtyr) then
                   loc_bio_remin_dt = loc_bio_remin_dt - ( junkst - dum_dtyr )                                          ! residence time in partial layer (yrs)
                   !                   loc_bio_remin_dD = phys_ocn(ipo_dD,dum_i,dum_j,kk)-( (junkst-dum_dtyr)*loc_bio_remin_sinkingrate )               
                   !                  print*,'dD1=',loc_bio_remin_dD
                   if( loc_bio_remin_dt < const_real_nullsmall) then
                      print*,'now dt>0: ',loc_bio_remin_dt
                      loc_bio_remin_dt = 0.0
                   endif
                   loc_bio_remin_dD = loc_bio_remin_dt*loc_bio_remin_sinkingrate                       ! layer thickness above/partial thickness below
                   !                   loc_bio_remin_layerratio = phys_ocn(ipo_dD,dum_i,dum_j,kk+1)/loc_bio_remin_dD     
                   !                print*,'now restime ',loc_bio_remin_dt,loc_bio_remin_dD
                endif
                ! *** Calculate fractional change in particulate fluxes ***
                ! particulate organic matter
                if (sed_select(is_POC)) then
                   If (.NOT. opt_bio(iopt_bio_remin_POC_fixed)) then
                      ! by M.Chikamoto 2006-06-30 
                      ! calculate residence time in each ocean layer

                      r0    = 0.1 * conv_yr_d         ! Yamanaka et al. (2004)     r0 = decay rate(based on ...)
                      ktemp = log(2.0)/(10.)            ! Eppley (1972)

                      loc_bio_remin_POC_frac1 = loc_bio_remin_dt*r0 &
                           & *en_tempwt*exp(ktemp*(ocn(io_T,dum_i,dum_j,kk)-273.15))   !this is done if not messing with the temperature (default)

#ifdef reminTconst
!seasonal constant 5/12/09
                      if (biostep_test == 0 ) then
                         print*,'error on seasonal previous state read'
                         print*,'biostep_stop = ',biostep_test
                         stop
                      endif
                      loc_bio_remin_POC_frac1 = loc_bio_remin_dt*r0 &
                           & *exp(ktemp*(ocn_T_season1(dum_i,dum_j,kk,biostep)-273.15))
#endif                      
                      if (loc_bio_remin_POC_frac1 > const_real_one) loc_bio_remin_POC_frac1 = 1.0

                      loc_bio_remin_POC_frac2 = 0. ! refractory

                   !km remin based on profile
                   else
                      loc_bio_remin_POC_frac1 = (1.0 - EXP(-loc_bio_remin_dD/par_bio_remin_POC_eL1))
                      loc_bio_remin_POC_frac2 = (1.0 - EXP(-loc_bio_remin_dD/par_bio_remin_POC_eL2))
                   end if
                   ! calculate the ratio of particulate tracer between layers
                      loc_bio_part_POC_ratio = 1.0 - &
                           & ( &
                           &   (1.0 - loc_bio_part_TMP(is_POC_frac2,kk+1))*loc_bio_remin_POC_frac1 + &
                           &          loc_bio_part_TMP(is_POC_frac2,kk+1) *loc_bio_remin_POC_frac2 &
                           & )
                   ! calculate potential oxidation capacity
                      loc_potO2cap = fun_potO2cap(ocn_select(:),ocn(:,dum_i,dum_j,kk),bio_remin(:,dum_i,dum_j,kk))
                   ! compare with potential oxygen availability and modify fraction remineralized accordingly
                      loc_O2demand = conv_sed_ocn(io_O2,is_POC)* &
                           & loc_bio_remin_layerratio*loc_bio_part_POC_ratio*loc_bio_part_TMP(is_POC,kk+1)
                   ! decrease effective remin fraction to take into account depletion of electon acceptors
                      if ((loc_O2demand > loc_potO2cap) .AND. (loc_O2demand > const_real_nullsmall)) then
                         loc_bio_part_POC_ratio = (loc_potO2cap/loc_O2demand)*loc_bio_part_POC_ratio
                      end if

                   ! calculate change in partitioning between different fractions
                      is = is_POC_frac2
                      if (loc_bio_part_TMP(is,kk+1) > const_real_nullsmall.and.loc_bio_part_POC_ratio > const_real_nullsmall) then
                         loc_bio_part_TMP(is,kk) =  &
                              & (1.0 - loc_bio_remin_POC_frac2)*loc_bio_part_TMP(is,kk+1)/loc_bio_part_POC_ratio
                      else
                         loc_bio_part_TMP(is,kk) = 0.0
                      end if
                end if
                ! carbonate
                if (sed_select(is_CaCO3)) then
                   If (.NOT. opt_bio(iopt_bio_remin_CaCO3_fixed)) then                           !usually done
                      ! by M.Chiamoto 2006-07-06
                      ktemp = log(2.0)/(10.)             ! Eppley (1972) =.069
                      r0    = dec_rate * conv_yr_d        ! 0.05 /day Yamanaka et al. (2004) {= decomposition rate of CaCO3 at 0C}
!org                      r0    = 0.2 * conv_yr_d        ! 0.05 /day Yamanaka et al. (2004) {= decomposition rate of CaCO3 at 0C}
                      loc_bio_remin_CaCO3_frac1 = loc_bio_remin_dt*r0 &
                           & * en_tempwt_remin*exp(ktemp*(ocn(io_T,dum_i,dum_j,kk)-273.15))
#ifdef reminTconst
!seasonal constant 5/12/09
                      if (biostep_test == 0 ) then
                         print*,'error on seasonal previous state read'
                         print*,'biostep_stop = ',biostep_test
                         stop
                      endif
                      loc_bio_remin_CaCO3_frac1 = loc_bio_remin_dt*r0 &
                           & *exp(ktemp*(ocn_T_season1(dum_i,dum_j,kk,biostep)-273.15))
#endif                      
                      if (loc_bio_remin_CaCO3_frac1 > const_real_one) loc_bio_remin_CaCO3_frac1 = 1.0 !kst added 8/08 NOTE:  MESMO (descr.paper)does NOT have this...see notes10/9/08
                      loc_bio_remin_CaCO3_frac2 = 0.   ! refractory

                   else                                                                         !for fixed-profile CaCO3 remineralization:
                      loc_bio_remin_CaCO3_frac1 = (1.0 - EXP(-loc_bio_remin_dD/par_bio_remin_CaCO3_eL1))
                      loc_bio_remin_CaCO3_frac2 = (1.0 - EXP(-loc_bio_remin_dD/par_bio_remin_CaCO3_eL2))
                   endif

                   ! calculate the ratio of particulate tracer between layers
                   loc_bio_part_CaCO3_ratio = 1.0 - &
                        & ( &
                        &   (1.0 - loc_bio_part_TMP(is_CaCO3_frac2,kk+1))*loc_bio_remin_CaCO3_frac1 + &
                        &   loc_bio_part_TMP(is_CaCO3_frac2,kk+1)*loc_bio_remin_CaCO3_frac2 &
                        & )
                   ! calculate change in partitioning between different fractions
                   is = is_CaCO3_frac2
                   if (loc_bio_part_TMP(is,kk+1) > const_real_nullsmall .and. loc_bio_part_CaCO3_ratio > const_real_nullsmall) then
!                         print*,'caco3_frac2(kk+1) > 0)',kk
                      loc_bio_part_TMP(is,kk) = &
                           & (1.0 - loc_bio_remin_CaCO3_frac2)*loc_bio_part_TMP(is,kk+1)/loc_bio_part_CaCO3_ratio
                   else
                      loc_bio_part_TMP(is,kk) = 0.0
                   end if
                end if

                ! opal
                if (sed_select(is_opal)) then
                   If (.NOT. opt_bio(iopt_bio_remin_opal_fixed)) then
                      ! set local variables - temperature (K) and silicic acid concentration (mol kg-1)
                      loc_T     = ocn(io_T,dum_i,dum_j,kk)
#ifdef reminTconst
!seasonal constant 5/12/09
                      if (biostep_test == 0 ) then
                         print*,'error on seasonal previous state read'
                         print*,'biostep_stop = ',biostep_test
                         stop
                      endif
                      loc_T     = ocn_T_season1(dum_i,dum_j,kk,biostep)
#endif
                      loc_SiO2  = ocn(io_SiO2,dum_i,dum_j,kk)
                      ! calculate opal equilibrium H4SiO4 saturation concentration
                      loc_Si_eq = conv_umol_mol*10.0**(6.44 - 968.0/loc_T)
                      ! calculate degree of opal undersatruation
                      loc_u     = (loc_Si_eq - loc_SiO2)/loc_Si_eq
                      IF (loc_u > const_real_one)       loc_u = 1.0
                      IF (loc_u < const_real_nullsmall) loc_u = 0.0

                      ! calculate opal fractional dissolution
                      ! NOTE: for now, assume that both opal 'fractionas' behave identically
!km                      print*,'km what is par_bio_remin_opal_K? ',par_bio_remin_opal_K
                      loc_bio_remin_opal_frac1 =                                             &
!km                           & loc_bio_remin_dt*par_bio_remin_opal_K*                          &
                           & loc_bio_remin_dt*(par_bio_remin_opal_K/conv_d_yr)*                          &
                           & (1.0/0.71)*                                                     &
                           & (                                                               &
                           &   (0.16*(1.0 + (loc_T - const_zeroC)/15.0)*loc_u) +             &
                           &   (0.55*((1.0 + (loc_T - const_zeroC)/400.0)**4.0*loc_u)**9.25) &
                           & )

                      if (loc_bio_remin_opal_frac1 > const_real_one) loc_bio_remin_opal_frac1 = 1.0
                      loc_bio_remin_opal_frac2 = loc_bio_remin_opal_frac1

                   else
                      loc_bio_remin_opal_frac1 = (1.0 - EXP(-loc_bio_remin_dD/par_bio_remin_opal_eL1))
                      loc_bio_remin_opal_frac2 = (1.0 - EXP(-loc_bio_remin_dD/par_bio_remin_opal_eL2))
                   endif


                   ! calculate the ratio of particulate tracer between layers
                   loc_bio_part_opal_ratio = 1.0 - &
                        & ( &
                        &   (1.0 - loc_bio_part_TMP(is_opal_frac2,kk+1))*loc_bio_remin_opal_frac1 + &
                        &   loc_bio_part_TMP(is_opal_frac2,kk+1)*loc_bio_remin_opal_frac2 &
                        & )
                   ! calculate change in partitioning between different fractions
                   is = is_opal_frac2
                   if (loc_bio_part_TMP(is,kk+1) > const_real_nullsmall .and. loc_bio_part_opal_ratio > const_real_nullsmall) then
                      loc_bio_part_TMP(is,kk) = &
                           & (1.0 - loc_bio_remin_opal_frac2)*loc_bio_part_TMP(is,kk+1)/loc_bio_part_opal_ratio
                   else
                      loc_bio_part_TMP(is,kk) = 0.0
                   end if

                end if
                

                ! *** Calculate particle concentrations in layer below ***
                ! calculate local (temporary) particulate tracer concentration;
                ! (a) take the particulate concentration in the layer above
                ! (b) modify it by the remineralization ratio to take into account dissolution loss, and
                ! (c) convert to the equivalent particulate concentration in the new (lower) layer
                ! NOTE: NO fractionation (in elemental composition or isotopic ratio) is currently assumed
                !       => additional CASE statements are required to deal with specific fractionation cases
                ! NOTE: adjust fraction of scavenged material that can be returned (par_scav_fremin)
                !       => e.g., par_scav_Fe_remin = 1.0 will result in
                !          scavanged material being returned in the same proportion as remineralization of the scavenger
                DO l=1,n_ismax
                   is = conv_iselected_is(l)
                   ! particulate organic matter (plus elemental components, and particle-reactive scavenged elements)
                   !by M.Chikamoto 07-26-2007 15N 
 
!                   if ((sed_dep(is) == is_POC) .or. (sed_type(is) == par_sed_type_POM) .OR. (sed_dep(is) == is_PON)) then 
                   if ((sed_dep(is) == is_POC) .or. (sed_type(is) == par_sed_type_POM)) then 
!                      loc_bio_part_TMP(is,kk) = loc_bio_remin_layerratio*loc_bio_part_POC_ratio*loc_bio_part_TMP(is,kk+1) 
!didn't have this until 6/24/10...
                      if (sed_type(is) == par_sed_type_scavenged) then
                         loc_bio_part_TMP(is,kk) = loc_bio_part_TMP(is,kk+1)* &
                              & loc_bio_remin_layerratio*(1.0 - par_scav_Fe_remin*(1.0 - loc_bio_part_POC_ratio))
                      else
                         loc_bio_part_TMP(is,kk) = loc_bio_part_TMP(is,kk+1)* &
                              & loc_bio_remin_layerratio*loc_bio_part_POC_ratio
                      end if
!didn't have this until 6/24/10...
                   end if 
                   if(is == is_PON_15N)then 
                     ! by M.Chikamoto 07-13-2006 adding 15N fractionation between PON_15N and NO3_15N through nitrification 
                     !                           PON_15N +17 permil  
                     !                           NO3_15N -17 permil  
                     !                               (Sigman et al., 2005, GBC,19) 
                     !  loc_bio_part_TMP(is_PON_15N,kk) ---- d[PON_15N] (= dissolution + isotope fraction) 
                     !  note: loc_bio_part_TMP[is_PON]      \= loc_bio_remin[io_NO3]   
                     !        loc_bio_part_TMP[is_PON_15N]  \= loc_bio_remin[io_NO3_15N]  why? 
                     ! 
                     ! We assume no fractionation through nitrification without using ammonia.==> loc_delta = 0 
                     loc_delta = 0.0 !+17. 
                     loc_alpha = 1.0 + loc_delta/1000. 
 
                     loc_bio_part_TMP(is,kk) = loc_bio_remin_layerratio*loc_bio_part_POC_ratio*loc_bio_part_TMP(is,kk+1) * loc_alpha 
                  endif 
                   ! carbonate (plus elemental components, and particle-reactive scavenged elements)
                   if ((sed_dep(is) == is_CaCO3) .OR. (sed_type(is) == par_sed_type_CaCO3)) then
                      loc_bio_part_TMP(is,kk) = loc_bio_remin_layerratio*loc_bio_part_CaCO3_ratio*loc_bio_part_TMP(is,kk+1)
                   end if
                   ! opal (plus elemental components, and particle-reactive scavenged elements)
                   if ((sed_dep(is) == is_opal) .OR. (sed_type(is) == par_sed_type_opal)) then                      
                      loc_bio_part_TMP(is,kk) = loc_bio_remin_layerratio*loc_bio_part_opal_ratio*loc_bio_part_TMP(is,kk+1)
                   end if
                   if ((sed_dep(is) == is_det) .OR. (sed_type(is) == par_sed_type_det)) then
                      loc_bio_part_TMP(is,kk) = loc_bio_remin_layerratio*loc_bio_part_TMP(is,kk+1)
                   end if
                end do
                
                ! *** Calculate increase in tracer concentrations due to particle remineralization ***
                ! >>> GENERIC ALGORITHM
                ! add 'missing' (remineralized) particulate sediment tracers to respective remineralization array components
                ! NOTE: ensure the particulate concentration in the upper layer is scaled w.r.t. 
                !       the diference in relative layer thickness
                !                    (DIC,ALK,Ca from CaCO3)
                DO l=1,n_ismax
                   is = conv_iselected_is(l)
                   loc_tot_i = conv_sed_ocn_i(0,is)
                   do loc_i=1,loc_tot_i
                      io = conv_sed_ocn_i(loc_i,is)
                      ! km - this is where the negative value comes in
                      loc_bio_remin(io,kk) = loc_bio_remin(io,kk) + conv_sed_ocn(io,is)* &
                           & (loc_bio_remin_layerratio*loc_bio_part_TMP(is,kk+1) - loc_bio_part_TMP(is,kk))
                        ! 2007/01 Chikamoto
                      if(io.eq.io_DIC.and.is.eq.is_POC)then
                         res_ocn(dum_i,dum_j,kk) = conv_sed_ocn(io,is)* & 
                           & (loc_bio_remin_layerratio*loc_bio_part_TMP(is,kk+1) - loc_bio_part_TMP(is,kk)) &
                           & /dum_dtyr*phys_ocn(ipo_M,dum_i,dum_j,k)*par_bio_red_POP_PON/par_bio_red_POP_POC ! molN yr-1  kst
                      endif
                   end do
                end DO

                ! *** Scavenge Fe from water column *** Copy from Andy's codes  Sun 2009/03/03
                ! NOTE: Fe scavenging must be called AFTER particulates have been 'moved' to the next layer down
                !       - they are assumed to start at the BASE of the originating layer,
                !         which is why they are not scavenged from level (kk+1)
!kst 4/29/10   restrict scavenging to non-production layers, as production layers are scavenged during uptake
                if (kk < n_kmax+1-nlayer_prod) then
                   if (ocn_select(io_Fe)) then
                      if (ocn(io_Fe,dum_i,dum_j,kk) > const_real_nullsmall) then
                            call sub_calc_scav_Fe( &
                                 & loc_bio_remin_dt, &
                                 & ocn(io_Fe,dum_i,dum_j,kk), &
                                 & loc_bio_part_TMP(:,kk), &
                                 & loc_bio_remin(:,kk) &
                                 & )
                      end if
                   end if
                endif
             else
                If (.NOT. opt_misc(iopt_misc_sed_select)) then
                   ! >>> GENERIC ALGORITHM BUT CONTAINS USER-MODIFIABLE INFORMATION IN ARRAY <conv_sed_ocn>
                   ! convert all 'missing' particulate sediment tracer reaching the sediments to (dissolved) tracer form 
                   ! and add to remineralization array
                   DO l=1,n_ismax
                      is = conv_iselected_is(l)
                      loc_tot_i = conv_sed_ocn_i(0,is)
                      do loc_i=1,loc_tot_i
                         io = conv_sed_ocn_i(loc_i,is)
                         loc_bio_remin(io,kk+1) = loc_bio_remin(io,kk+1) + conv_sed_ocn(io,is)* &
                              & loc_bio_part_TMP(is,kk+1)
                      end do
                   end DO
                else
                end If
             end If
          end do k2loop


          ! <<<<<<<<<<<<<<<<<<<<<<<
          ! *** kk SUB-LOOP END ***
          ! <<<<<<<<<<<<<<<<<<<<<<<

          ! *** UPDATE PARTICULATE MATTER INFORMATION ***
          ! >>> GENERIC ALGORITHM BUT CONTAINS USER-MODIFIABLE INFORMATION IN ARRAY <conv_sed_ocn>
          ! update local ocean particulate tracer field - store residual particulate tracer at the point of 
          ! the deepest level reached
          ! NOTE: do not store if the sediment surface is reached
          If (loc_bio_remin_min_k >= dum_k1) then
             DO l=1,n_ismax
                is = conv_iselected_is(l)
                ! km - frac2 from nlayer_prod=2 are added to give erroneously > 1
                loc_bio_part(is,loc_bio_remin_min_k) = loc_bio_part(is,loc_bio_remin_min_k) + &
                     & loc_bio_part_TMP(is,loc_bio_remin_min_k)
             end do
          end if

          ! >>> NON-GENERIC ALGORITHM
          ! record particulate fluxes at base of each layer (units of; mol per time-step)
          ! NOTE: implicitly includes sedimentation flux (kk=dum_k1)\
          ! NOTE: exclude meaningless 
          do kk=k,loc_bio_remin_min_k+1,-1
             DO l=1,n_ismax
                is = conv_iselected_is(l)
                SELECT CASE (sed_type(is))
                case (par_sed_type_frac)
                   loc_bio_settle(is,kk) = loc_bio_settle(is,kk) + &
                        & loc_bio_part_TMP(is,kk)
                case default
                   loc_bio_settle(is,kk) = loc_bio_settle(is,kk) + &
                        & phys_ocn(ipo_M,dum_i,dum_j,kk)*loc_bio_part_TMP(is,kk)
                end SELECT
             end do
          end do
       else
       end IF
    end do kloop

    ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
    ! *** k WATER-COLUMN LOOP END ***
    ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

    ! *** WRITE GLOBAL ARRAY DATA ***
    ! NOTE: because sub_calc_bio_reminis the first call in the sequence of events (in biogem_main),
    !       data arrays are over-written rather than incremented
    ! write ocean tracer field and settling flux arrays (global array)
    DO l=1,n_ismax
       is = conv_iselected_is(l)

       ! km rectify the fractions for nlayer_prod>1
       select case (sed_type(is))
       case (par_sed_type_frac)       ! this happens only for POC_frac2, CaCO3_frac2, and opal_frac2    
          loc_bio_part(is,:)=loc_bio_part(is,:)/nlayer_prod     
!   kst:  this takes the average of frac2 for production layers>  
!           I don't think it matters when not using prescribed remin profile because w/o prescribed profile, frac2 = 0 anyway.
       end select

       bio_part(is,dum_i,dum_j,:)   = loc_bio_part(is,:)
       bio_settle(is,dum_i,dum_j,:) = loc_bio_settle(is,:)
    end do       

    tv1 = SUM(bio_part_lg(dum_i,dum_j,15:16))            !diatom orgC
    tv2 = SUM(bio_part_sm(dum_i,dum_j,15:16))            !sm phyt C
    if (tv1 > const_real_nullsmall .or. tv2 > const_real_nullsmall) then
       bio_settle_lg(dum_i,dum_j,:) = loc_bio_settle(is_POC,:)*tv1/(tv1+tv2)
       bio_settle_sm(dum_i,dum_j,:) = loc_bio_settle(is_POC,:)*tv2/(tv1+tv2)
    else
       bio_settle_lg(dum_i,dum_j,:) = 0.0
       bio_settle_sm(dum_i,dum_j,:) = 0.0
    endif

  ! write ocean tracer remineralization field (global array)
    DO l=3,n_iomax
       io = conv_iselected_io(l)
       bio_remin(io,dum_i,dum_j,:) = bio_remin(io,dum_i,dum_j,:) + loc_bio_remin(io,:)
    end do

  END SUBROUTINE sub_calc_bio_remin


  ! *** calculate nitrate reduction arrising from anarobic organic matter remineralization *** 
  ! <<< CONTAINS NON-GENERIC ALGORITHMS >>> 
  SUBROUTINE sub_calc_bio_remin_reduce_NO3(dum_i,dum_j,dum_k1,loc_dtyr) 
    ! dummy arguments 
    INTEGER,INTENT(in)::dum_i,dum_j,dum_k1 
    ! local variables 
    integer::l,io,k 
    real::loc_potO2,loc_potO2def
    real::loc_NO3 
    real,dimension(0:n_ocn,n_maxk)::loc_bio_remin 
 
    ! by M.Chikamoto 07-13-2006  
    real::loc_r15N
    real::loc_delta, loc_alpha 
    real::loc_dtyr
    ! Chikamoto 2007-01
    real::loc_den
 
    ! *** INITIALIZE VARIABLES *** 
    ! initialize local variables 
    ! initialize remineralization tracer arrays 
    DO l=3,n_iomax 
       io = conv_iselected_io(l) 
       loc_bio_remin(io,:) = 0.0 
    end do 
 
    ! *** NITRATE REDUCTION *** 
    ! look for oxygen deficits; carry out nitrate reduction to remove the O2 deficit 
    ! NOTE: 2 moles of NO3 are required for every 3 moles of O2 deficit (hence the factor 1.5) 
    !       2NO3 -> N2 + 3O2 
    DO k=n_kmax,dum_k1,-1 
       ! calculate potential oxygen availability 
       loc_potO2 = ocn(io_O2,dum_i,dum_j,k) + bio_remin(io_O2,dum_i,dum_j,k) 
       ! if oxygen deficit exists, go do something useful about it 
       loc_den = 0.
       if ((loc_potO2 < -const_real_nullsmall) .AND. (ocn(io_NO3,dum_i,dum_j,k) > const_real_nullsmall)) then 
 
          ! set local NO3 concentration 
          loc_NO3 = ocn(io_NO3,dum_i,dum_j,k) 
          ! calculate potential oxygen  deficit 
          loc_potO2def = abs(loc_potO2) 

          ! by M.Chikamoto 07-13-2006 calculating 15N fractionation between NO3_15N and N2_15N 
          !                           NO3_15N +25 permil  
          !                           N2_15N  -25 permil  
          !                               (Sigman et al., 2005, GBC,19) 
          loc_delta = -25.0 !-25. ; km 8Mar07 set to 0 to debug
          loc_alpha = 1.0 + loc_delta/1000. 
 
          if (loc_potO2def < 1.5*loc_NO3) then 
             ! partial NO3 utilization (=> N isotope Rayleigh fractionation) 
             loc_bio_remin(io_NO3,k) =  -(2.0/3.0)*loc_potO2def 
             loc_bio_remin(io_N2,k)  = 0.5*(2.0/3.0)*loc_potO2def 
             loc_bio_remin(io_O2,k)  = -loc_potO2 

!km             loc_den = loc_bio_remin(io_O2,k) * par_red_PC / par_red_PO ! molC kg-1
             loc_den = loc_bio_remin(io_O2,k) * par_bio_red_POP_POC / (-1.*par_bio_red_POP_PO2) ! molN kg-1
          else 
             ! complete NO3 utilization (no N fractionation) 
             loc_bio_remin(io_O2,k)      = 1.5*loc_NO3 
             loc_bio_remin(io_NO3,k)     = -loc_NO3 
             loc_bio_remin(io_N2,k)      = 0.5*loc_NO3 
 
!km             loc_den = - loc_bio_remin(io_NO3,k) * 3./2. * par_red_PC / par_red_PO ! mol kg-1
             loc_den = - loc_bio_remin(io_NO3,k) * 3./2. * par_bio_red_POP_PON / (-1.*par_bio_red_POP_PO2) ! molN kg-1
          end if

          ! by Chikamoto 07-13-2006 ; km Mar2007 - taken outside above if clause
          if(ocn_select(io_NO3_15N))then 
             loc_r15N = 0. 
             if(ocn(io_NO3,dum_i,dum_j,k) > const_real_nullsmall)then 
                loc_r15N = ocn(io_NO3_15N,dum_i,dum_j,k)/ocn(io_NO3,dum_i,dum_j,k) 
             endif 

             loc_bio_remin(io_NO3_15N,k) = loc_bio_remin(io_NO3,k)*loc_r15N*loc_alpha
             loc_bio_remin(io_N2_15N,k)  = -1./2. * loc_bio_remin(io_NO3_15N,k)
          endif 

       end if
       den_ocn(dum_i,dum_j,k) = loc_den /loc_dtyr*phys_ocn(ipo_M,dum_i,dum_j,k) ! molN yr-1
    end do
 
    ! *** WRITE GLOBAL ARRAY DATA *** 
    ! write ocean tracer remineralization field (global array) 
    DO l=3,n_iomax 
       io = conv_iselected_io(l) 
       bio_remin(io,dum_i,dum_j,:) = bio_remin(io,dum_i,dum_j,:) + loc_bio_remin(io,:) 
    end do 
 
  end SUBROUTINE sub_calc_bio_remin_reduce_NO3 
 


  ! *** calculate sulphate reduction arrising from anarobic organic matter remineralization ***
  ! <<< CONTAINS NON-GENERIC ALGORITHMS >>>
  SUBROUTINE sub_calc_bio_remin_reduce_SO4(dum_i,dum_j,dum_k1)
    ! dummy arguments
    INTEGER,INTENT(in)::dum_i,dum_j
    integer,intent(in)::dum_k1
    ! local variables
    integer::l,io,is
    INTEGER::k
    real::loc_potO2,loc_potO2def,loc_r34S
    real::loc_SO4
    real,dimension(0:n_ocn,n_maxk)::loc_bio_remin

    ! *** INITIALIZE VARIABLES ***
    ! initialize local variables
    ! initialize remineralization tracer arrays
    DO l=3,n_iomax
       io = conv_iselected_io(l)
       loc_bio_remin(io,:) = 0.0
    end do

    ! *** SULPHATE REDUCTION ***
    ! >>> NON-GENERIC ALGORITHM
    ! look for oxygen deficits; carry out sulphate reduction to remove the O2 deficit
    ! *** sulphate reduction ***
    ! NOTE: 1 mole of SO4 are required for every 2 moles of O2 deficit (hence the factor 2.0)
    !       2H+ + SO4 -> H2S + 2O2
    DO k=n_kmax,dum_k1,-1
       ! calculate potential oxygen availability
       loc_potO2 = ocn(io_O2,dum_i,dum_j,k) + bio_remin(io_O2,dum_i,dum_j,k)
       ! if oxygen deficit exists, go do something useful about it
       If ((loc_potO2 < -const_real_nullsmall) .AND. (ocn(io_SO4,dum_i,dum_j,k) > const_real_nullsmall)) then
          ! set local NO3 concentration
          loc_SO4 = ocn(io_SO4,dum_i,dum_j,k)
          ! calculate potential oxygen deficit
          loc_potO2def = abs(loc_potO2)
          ! calculate isotopic ratio
          loc_r34S = ocn(io_SO4_34S,dum_i,dum_j,k)/ocn(io_SO4,dum_i,dum_j,k)
          if (loc_potO2def < 2.0*loc_SO4) then
             ! partial SO4 utilization (=> S isotope Rayleigh fractionation)
             loc_bio_remin(io_SO4,k) = -0.5*loc_potO2def
             loc_bio_remin(io_H2S,k) = 0.5*loc_potO2def
             loc_bio_remin(io_O2,k)  = -loc_potO2
             ! \/\/\/ INSERT ALTERNATIVE CODE FOR NON-ZERO S FRACTIONATION \/\/\/
             ! [currently there is no fractionation occurring]
             loc_bio_remin(io_SO4_34S,k) = -loc_r34S*0.5*loc_potO2def
             loc_bio_remin(io_H2S_34S,k) = loc_r34S*0.5*loc_potO2def
             ! /\/\/\ INSERT ALTERNATIVE CODE FOR NON-ZERO S FRACTIONATION /\/\/\
          else
             ! complete SO4 utilization (no S fractionation)
             loc_bio_remin(io_O2,k)      = 2.0*loc_SO4
             loc_bio_remin(io_SO4,k)     = -loc_SO4
             loc_bio_remin(io_H2S,k)     = loc_SO4
             loc_bio_remin(io_SO4_34S,k) = -loc_r34S*loc_SO4
             loc_bio_remin(io_H2S_34S,k) = loc_r34S*loc_SO4
          end if
       end if
    end do

    ! *** WRITE GLOBAL ARRAY DATA ***
    ! write ocean tracer remineralization field (global array)
    DO l=3,n_iomax
       io = conv_iselected_io(l)
       bio_remin(io,dum_i,dum_j,:) = bio_remin(io,dum_i,dum_j,:) + loc_bio_remin(io,:)
    end do

  end SUBROUTINE sub_calc_bio_remin_reduce_SO4


 ! ****************************************************************************************************************************** !
 ! Copy from Andy's Sun 2009/03/03
 ! ****************************************************************************************************************************** !
  ! Calculate Fe scavenging
  SUBROUTINE sub_calc_scav_Fe(dum_dtyr,dum_ocn_Fe,dum_bio_part,dum_bio_remin)
    ! dummy arguments
    REAL,INTENT(in)::dum_dtyr
    REAL,INTENT(in)::dum_ocn_Fe
    real,dimension(0:n_sed),INTENT(inout)::dum_bio_part    !5/6/10  ok, 0:n_sed instead of just n_sed straightened out array dimensions...
    real,dimension(0:n_ocn),INTENT(inout)::dum_bio_remin
    ! local variables
    real::loc_scav_Fe_k_POC,loc_scav_Fe_k_CaCO3,loc_scav_Fe_k_opal,loc_scav_Fe_k_det
    real::loc_scav_Fe_k_tot
    real::loc_scav_dFe_tot
    real::loc_part_den_POC,loc_part_den_CaCO3,loc_part_den_opal,loc_part_den_det
    real::loc_part_den_tot

    ! *** Calculate Fe scavenging ***
    ! NOTE: residence time in each ocean layer must be estimated for the Parekh et al. [2005] model
    ! NOTE: Dutkiewicz et al. [2005] scavenging rate par_scav_Fe_Ks has been converted to units of (yr-1)
    ! NOTE: scavening from surface ocean is NOT calculated here ...kst:  well, actually, it is if you called this at the surface (as in bio_uptake)
    if (sed_select(is_POM_Fe)) then 
       loc_part_den_POC = conv_g_mg*conv_POC_mol_g*dum_bio_part(is_POC)/conv_kg_l
    else
       loc_part_den_POC = 0.0
    end if
    if (sed_select(is_CaCO3_Fe)) then 
       loc_part_den_CaCO3 = conv_g_mg*conv_cal_mol_g*dum_bio_part(is_CaCO3)/conv_kg_l
    else
       loc_part_den_CaCO3 = 0.0
    end if
    if (sed_select(is_opal_Fe)) then 
       loc_part_den_opal  = conv_g_mg*conv_opal_mol_g*dum_bio_part(is_opal)/conv_kg_l
    else
       loc_part_den_opal  = 0.0
    end if
    if (sed_select(is_det_Fe)) then 
       loc_part_den_det   = conv_g_mg*conv_det_mol_g*dum_bio_part(is_det)/conv_kg_l
    else
       loc_part_den_det   = 0.0
    endif
    loc_part_den_tot = loc_part_den_POC + loc_part_den_CaCO3 + loc_part_den_opal + loc_part_den_det 
    if (loc_part_den_tot > const_real_nullsmall) then
       if (opt_bio(iopt_bio_Fe_fixedKscav)) then
          ! calculate scavenging following Dutkiewicz et al. [2005]
          loc_scav_Fe_k_tot = par_scav_Fe_ks
          loc_scav_Fe_k_POC = (loc_part_den_POC/loc_part_den_tot)*loc_scav_Fe_k_tot
          loc_scav_Fe_k_CaCO3 = (loc_part_den_CaCO3/loc_part_den_tot)*loc_scav_Fe_k_tot
          loc_scav_Fe_k_opal = (loc_part_den_opal/loc_part_den_tot)*loc_scav_Fe_k_tot
          loc_scav_Fe_k_det = (loc_part_den_det/loc_part_den_tot)*loc_scav_Fe_k_tot
          ! calculate total Fe scavenged  in mol/kg
          loc_scav_dFe_tot = dum_dtyr*loc_scav_Fe_k_tot*dum_ocn_Fe
       else
          ! calculate scavenging following Parekh et al. [2005]
          loc_scav_Fe_k_POC = par_scav_Fe_sf_POC*par_scav_Fe_k0*loc_part_den_POC**par_scav_Fe_exp
          loc_scav_Fe_k_CaCO3 = par_scav_Fe_sf_CaCO3*par_scav_Fe_k0*loc_part_den_CaCO3**par_scav_Fe_exp
          loc_scav_Fe_k_opal = par_scav_Fe_sf_opal*par_scav_Fe_k0*loc_part_den_opal**par_scav_Fe_exp
          loc_scav_Fe_k_det = par_scav_Fe_sf_det*par_scav_Fe_k0*loc_part_den_det**par_scav_Fe_exp
          loc_scav_Fe_k_tot = loc_scav_Fe_k_POC + loc_scav_Fe_k_CaCO3 + loc_scav_Fe_k_opal + loc_scav_Fe_k_det
          ! calculate total Fe scavenged in mol/kg
          loc_scav_dFe_tot = dum_dtyr*loc_scav_Fe_k_tot*dum_ocn_Fe
       end if
       ! calculate Fe scavenged by particulates
       ! and update local remineralization array to take into account the removal of Fe from solution
       dum_bio_part(is_POM_Fe)   = dum_bio_part(is_POM_Fe) + (loc_scav_Fe_k_POC/loc_scav_Fe_k_tot)*loc_scav_dFe_tot
       dum_bio_part(is_CaCO3_Fe) = dum_bio_part(is_CaCO3_Fe) + (loc_scav_Fe_k_CaCO3/loc_scav_Fe_k_tot)*loc_scav_dFe_tot
       dum_bio_part(is_opal_Fe)  = dum_bio_part(is_opal_Fe) + (loc_scav_Fe_k_opal/loc_scav_Fe_k_tot)*loc_scav_dFe_tot
       dum_bio_part(is_det_Fe)   = dum_bio_part(is_det_Fe) + (loc_scav_Fe_k_det/loc_scav_Fe_k_tot)*loc_scav_dFe_tot
       dum_bio_remin(io_Fe) = dum_bio_remin(io_Fe) - loc_scav_dFe_tot
    end if

  end SUBROUTINE sub_calc_scav_Fe
  ! ****************************************************************************************************************************** !



  ! *** correct for spurious negative H2S concentrations ***
  ! <<< CONTAINS NON-GENERIC ALGORITHMS >>>
  SUBROUTINE sub_calc_bio_remin_correct_H2S(dum_i,dum_j,dum_k1)
    ! dummy arguments
    INTEGER,INTENT(in)::dum_i,dum_j,dum_k1
    ! local variables
    integer::l,io,k
    real::loc_potH2S,loc_r34S
    real::loc_H2S_oxidation_const,loc_H2S_oxidation
    real,dimension(0:n_ocn,n_maxk)::loc_bio_remin

    ! *** INITIALIZE VARIABLES ***
    ! initialize local variables
    ! initialize remineralization tracer arrays
    DO l=3,n_iomax
       io = conv_iselected_io(l)
       loc_bio_remin(io,:) = 0.0
    end do
    
    ! *** 'FIX' H2S ... ***
    ! 'fix' to prevent 'excessive' negative H2S developing
    ! NOTE: this is not really consistent with the philosophy of the numerical time-stepping scheme, but so what
    ! NOTE: after LRK
    DO k=n_kmax,dum_k1,-1
       ! calculate potential H2S
       loc_potH2S = ocn(io_H2S,dum_i,dum_j,k)
       if (loc_potH2S < -const_real_nullsmall) then
          ! calculate isotopic ratio
          loc_r34S = ocn(io_H2S_34S,dum_i,dum_j,k)/ocn(io_H2S,dum_i,dum_j,k)
          ! just crudely remove all negative H2S if there is sufficient SO4 to 'neutralize' it
          if (-loc_potH2S < ocn(io_SO4,dum_i,dum_j,k)) then
             loc_bio_remin(io_H2S,k)     = -loc_potH2S
             loc_bio_remin(io_SO4,k)     = loc_potH2S
             loc_bio_remin(io_O2,k)      = -2.0*loc_potH2S
             loc_bio_remin(io_H2S_34S,k) = -loc_r34S*loc_potH2S
             loc_bio_remin(io_SO4_34S,k) = loc_r34S*loc_potH2S
          end if
       end if
    end DO

    ! *** WRITE GLOBAL ARRAY DATA ***
    ! write ocean tracer remineralization field (global array)
    DO l=3,n_iomax
       io = conv_iselected_io(l)
       bio_remin(io,dum_i,dum_j,:) = bio_remin(io,dum_i,dum_j,:) + loc_bio_remin(io,:)
    end do
    
  end SUBROUTINE sub_calc_bio_remin_correct_H2S


  ! *********************************
  ! *** FORCING FUNCTION ROUTINES ***
  ! *********************************


  ! *** updating environment at the current (BioGeM) model time w.r.t. a defined signal function ***
  ! <<< GENERIC >>>
  SUBROUTINE sub_update_sig(dum_t,dum_sig,dum_sig_i,dum_x)
    ! dummy arguments
    REAL, INTENT(in)::dum_t
    REAL,INTENT(in),DIMENSION(n_data_max)::dum_sig
    INTEGER,INTENT(inout),DIMENSION(2)::dum_sig_i
    REAL, INTENT(out)::dum_x
    ! update forcing signal indices (if required) and carry put linear interpolation
    ! NOTE: t(1) is the lower age bounding point, t(2) is the upper age bounding point
    IF (dum_sig_i(1) > 1) THEN
       IF (dum_t < dum_sig(dum_sig_i(1))) THEN
          DO
             dum_sig_i(1) = dum_sig_i(1) - 1
             IF (dum_t > dum_sig(dum_sig_i(1))) THEN
                ! found correct index - exit loop
                EXIT
             ELSEIF (dum_sig_i(1) == 1) THEN
                EXIT
             END IF
          END DO
       END IF
    END IF
    IF (dum_sig_i(2) > 1) THEN
       IF (dum_t < dum_sig(dum_sig_i(2))) THEN
          DO
             dum_sig_i(2) = dum_sig_i(2) - 1
             IF (dum_t > dum_sig(dum_sig_i(2))) THEN
!trunk 4/21/10 uses >=   kst
!             IF (dum_t >= dum_sig(dum_sig_i(2))) THEN
                ! come too far - add one back to index(2) and exit loop
                dum_sig_i(2) = dum_sig_i(2) + 1
                EXIT
             ELSEIF (dum_sig_i(2) == 1) THEN
                EXIT
             END IF
          END DO
       END IF
    END IF
    ! calculate relative position of current time w.r.t. upper and lower bounding points of the signal function
    ! NOTE: if upper and lower bounding points are identical 
    !       (i.e., if current time is outside of maximum or minimum signal time values)
    !       avoid divide-by-zero problems and assume a value of 0.5
    IF (ABS(dum_sig(dum_sig_i(2)) - dum_sig(dum_sig_i(1))) > const_real_nullsmall) THEN
       dum_x = (dum_sig(dum_sig_i(2)) - dum_t)/(dum_sig(dum_sig_i(2)) - dum_sig(dum_sig_i(1)))
    ELSE
       dum_x = 0.5
    ENDIF
  END SUBROUTINE sub_update_sig


  ! *** update ocean restoring forcing function value ***
  ! <<< GENERIC >>>
  SUBROUTINE sub_update_force_restore_ocn(dum_t,dum_io)
    ! dummy arguments
    REAL,INTENT(in)::dum_t
    INTEGER,INTENT(in)::dum_io
    ! local variables
    INTEGER::i,j,k
    REAL::loc_x
    REAL::loc_force_restore_ocn
    real::loc_tot,loc_standard
    ! calculate new forcing time series values
    CALL sub_update_sig(dum_t,force_restore_ocn_sig(dum_io,1,:),force_restore_ocn_sig_i(dum_io,:),loc_x)
    force_restore_ocn_sig_x(dum_io) = &
         & (1 - loc_x)*force_restore_ocn_sig(dum_io,2,force_restore_ocn_sig_i(dum_io,2)) + &
         & loc_x*force_restore_ocn_sig(dum_io,2,force_restore_ocn_sig_i(dum_io,1))
    ! *** update prescribed (restoring) boundary conditions ***
    ! NOTE: use different <k> limits for the ocean restoring forcing loop (to enable surface-only forcing to be implemented)
    ! NOTE: flux forcings are in units of mol a-1
    DO i=1,n_imax
       DO j=1,n_jmax
          DO k=force_restore_ocn_k1(dum_io,i,j),n_kmax
             loc_force_restore_ocn = &
                  & force_restore_ocn_I(dum_io,i,j,k) + &
                  & force_restore_ocn_sig_x(dum_io)*(force_restore_ocn_II(dum_io,i,j,k) - force_restore_ocn_I(dum_io,i,j,k))
             SELECT CASE (ocn_type(dum_io))
             CASE (0,1)
                force_restore_ocn(dum_io,i,j,k) = loc_force_restore_ocn
!             case (11,12,13,14) 
                ! by M. Chikamoto 07-11-2006 adding 30Si 
              case (11,12,13,14,16) 
                loc_tot  = force_restore_ocn(ocn_dep(dum_io),i,j,k)
                loc_standard = const_standards(ocn_type(dum_io))
                force_restore_ocn(dum_io,i,j,k) = fun_calc_isotope_fraction(loc_force_restore_ocn,loc_standard)*loc_tot
             END SELECT
          END DO
       END DO
    END DO
  END SUBROUTINE sub_update_force_restore_ocn


  ! *** update ocean flux forcing function value ***
  ! <<< GENERIC >>>
  SUBROUTINE sub_update_force_flux_ocn(dum_t,dum_io)
    ! dummy arguments
    REAL,INTENT(in)::dum_t
    INTEGER,INTENT(in)::dum_io
    ! local variables
    INTEGER::i,j,k
    REAL::loc_x
    REAL::loc_force_flux_ocn_tot
    REAL::loc_force_flux_ocn_rtot
    REAL::loc_force_flux_ocn
    real::loc_tot,loc_standard
    ! calculate new forcing time series values
    CALL sub_update_sig(dum_t,force_flux_ocn_sig(dum_io,1,:),force_flux_ocn_sig_i(dum_io,:),loc_x)
    force_flux_ocn_sig_x(dum_io) = &
         & (1 - loc_x)*force_flux_ocn_sig(dum_io,2,force_flux_ocn_sig_i(dum_io,2)) + &
         & loc_x*force_flux_ocn_sig(dum_io,2,force_flux_ocn_sig_i(dum_io,1))
    ! *** update flux boundary conditions ***
    ! NOTE: use different <k> limits for the ocean restoring forcing loop (to enable surface-only forcing to be implemented)
    ! NOTE: flux forcings are in units of mol yr-1
    DO i=1,n_imax
       DO j=1,n_jmax
          DO k=goldstein_k1(i,j),n_kmax
             loc_force_flux_ocn = &
                  & force_flux_ocn_I(dum_io,i,j,k) + &
                  & force_flux_ocn_sig_x(dum_io)*(force_flux_ocn_II(dum_io,i,j,k) - force_flux_ocn_I(dum_io,i,j,k))
             SELECT CASE (ocn_type(dum_io))
             CASE (0,1)
                force_flux_ocn(dum_io,i,j,k) = loc_force_flux_ocn
!             case (11,12,13,14) 
                ! by M. Chikamoto 07-11-2006 adding 30Si 
             case (11,12,13,14,16) 
                loc_tot  = force_flux_ocn(ocn_dep(dum_io),i,j,k)
                loc_standard = const_standards(ocn_type(dum_io))
                force_flux_ocn(dum_io,i,j,k) = fun_calc_isotope_fraction(loc_force_flux_ocn,loc_standard)*loc_tot
             END SELECT
          END DO
       END DO
    END DO
    ! normalize flux forcings (if selected) so that the total flux is equal to the magnitude (at the current time step) 
    ! defined in the forcing signal file
    IF (force_flux_ocn_scale(dum_io)) THEN
       loc_force_flux_ocn_tot = SUM(force_flux_ocn(dum_io,:,:,:))
       if (abs(loc_force_flux_ocn_tot) > const_real_nullsmall) then
          loc_force_flux_ocn_rtot = 1.0/loc_force_flux_ocn_tot
       else
          loc_force_flux_ocn_rtot = 0.0
       end if
       DO i=1,n_imax
          DO j=1,n_jmax
             DO k=goldstein_k1(i,j),n_kmax
                SELECT CASE (ocn_type(dum_io))
                CASE (0,1)
                   force_flux_ocn(dum_io,i,j,k) = force_flux_ocn_sig_x(dum_io)*force_flux_ocn(dum_io,i,j,k)*loc_force_flux_ocn_rtot
                end SELECT
             END DO
          END DO
       END DO
    END IF
  END SUBROUTINE sub_update_force_flux_ocn


  ! *** update atmosphere tracer restoring forcing function value ***
  ! km 12CO2 must be restored, when 13CO2 and 14CO2 are restored
  ! <<< GENERIC >>>
  SUBROUTINE sub_update_force_restore_atm(dum_t,dum_ia)
!km  SUBROUTINE sub_update_force_restore_atm(dum_t,dum_ia,dum_atm)
    ! dummy arguments
    REAL,INTENT(in)::dum_t
    INTEGER,INTENT(in)::dum_ia
    ! local variables
    INTEGER::i,j
    REAL::loc_x
    REAL::loc_force_restore_atm
    real::loc_tot,loc_standard
    ! km 07 Mar 28 -local vars to restore 14C correctly (ie., BIG-D14C, not little-d14C)
    real::loc_littled14C
    real::loc_d13C, loc_frac13
    ! calculate new atmosphere forcing time series values
    CALL sub_update_sig(dum_t,force_restore_atm_sig(dum_ia,1,:),force_restore_atm_sig_i(dum_ia,:),loc_x)
    force_restore_atm_sig_x(dum_ia) = &
         & (1 - loc_x)*force_restore_atm_sig(dum_ia,2,force_restore_atm_sig_i(dum_ia,2)) + &
         & loc_x*force_restore_atm_sig(dum_ia,2,force_restore_atm_sig_i(dum_ia,1))
    ! update prescribed (restoring) boundary conditions
    DO i=1,n_imax
       DO j=1,n_jmax
          loc_force_restore_atm =  &
               & force_restore_atm_I(dum_ia,i,j) + &
               & force_restore_atm_sig_x(dum_ia)*(force_restore_atm_II(dum_ia,i,j) - force_restore_atm_I(dum_ia,i,j))
          SELECT CASE (atm_type(dum_ia))
          CASE (1)
             force_restore_atm(dum_ia,i,j) = loc_force_restore_atm
          !km case (11,12,13,14)
          case (11,13,14)
             loc_tot  = force_restore_atm(atm_dep(dum_ia),i,j)      ! main tracer has to be restored as well
             loc_standard = const_standards(atm_type(dum_ia))
             force_restore_atm(dum_ia,i,j) = fun_calc_isotope_fraction(loc_force_restore_atm,loc_standard)*loc_tot
          !km 28-Mar-08 restore atm 14C to BIG-D14C, not little-d14C
          case (12)
             loc_tot  = force_restore_atm(atm_dep(dum_ia),i,j)      ! main tracer has to be restored as well
             loc_frac13 = force_restore_atm(dum_ia-1,i,j)                       !  13C must also be restored, mol-13C/kg
             loc_standard = const_standards(atm_type(dum_ia-1))                 !  13C standard=0.011237  
             loc_d13c = fun_calc_isotope_delta(loc_tot,loc_frac13,loc_standard) ! d13C

             loc_littled14C = fun_convert_D14Ctodelta14C(loc_d13c,loc_force_restore_atm)
             loc_standard = const_standards(atm_type(dum_ia))                   !  14C standard=1.117d-12

             force_restore_atm(dum_ia,i,j) = fun_calc_isotope_fraction(loc_littled14C,loc_standard)*loc_tot
          END SELECT
       END DO
    END DO
  END SUBROUTINE sub_update_force_restore_atm


  ! *** update atmosphere tracer flux forcing function value ***
  ! <<< GENERIC >>>
  SUBROUTINE sub_update_force_flux_atm(dum_t,dum_ia)
    ! dummy arguments
    REAL,INTENT(in)::dum_t
    INTEGER,INTENT(in)::dum_ia
    ! local variables
    INTEGER::i,j
    REAL::loc_x
    REAL::loc_force_flux_atm_tot
    REAL::loc_force_flux_atm_rtot
    REAL::loc_force_flux_atm
    real::loc_tot,loc_standard
    ! calculate new atmosphere forcing time series values
    CALL sub_update_sig(dum_t,force_flux_atm_sig(dum_ia,1,:),force_flux_atm_sig_i(dum_ia,:),loc_x)
    force_flux_atm_sig_x(dum_ia) = &
         & (1 - loc_x)*force_flux_atm_sig(dum_ia,2,force_flux_atm_sig_i(dum_ia,2)) + &
         & loc_x*force_flux_atm_sig(dum_ia,2,force_flux_atm_sig_i(dum_ia,1))
    ! update flux boundary conditions
    ! NOTE: flux forcings are in units of mol yr-1
    DO i=1,n_imax
       DO j=1,n_jmax
          loc_force_flux_atm = &
               & force_flux_atm_I(dum_ia,i,j) + &
               & force_flux_atm_sig_x(dum_ia)*(force_flux_atm_II(dum_ia,i,j) - force_flux_atm_I(dum_ia,i,j))
          SELECT CASE (atm_type(dum_ia))
          CASE (1)
             force_flux_atm(dum_ia,i,j) = loc_force_flux_atm
          case (11,12,13,14)
             loc_tot  = force_flux_atm(atm_dep(dum_ia),i,j)
             loc_standard = const_standards(atm_type(dum_ia))
             force_flux_atm(dum_ia,i,j) = fun_calc_isotope_fraction(loc_force_flux_atm,loc_standard)*loc_tot
          END SELECT
       END DO
    END DO
    ! normalize flux forcings (if selected) so that the total flux is equal to the magnitude (at the current time step) 
    ! defined in the forcing signal file
    ! NOTE: only re-scale type 1 atmosphere tracers -
    !       the isotopic tracers will be automatically normalized because they are related directly to the total flux,
    !       and when this subroutine call is made for an isotopic tracer,
    !       it has already been called to deal with the related bulk tracer where the normalization is done
    IF (force_flux_atm_scale(dum_ia)) THEN
       loc_force_flux_atm_tot = SUM(force_flux_atm(dum_ia,:,:))
       if (abs(loc_force_flux_atm_tot) > const_real_nullsmall) then
          loc_force_flux_atm_rtot = 1.0/loc_force_flux_atm_tot
       else
          loc_force_flux_atm_rtot = 0.0
       end if
       DO i=1,n_imax
          DO j=1,n_jmax
             SELECT CASE (atm_type(dum_ia))
             CASE (1)
                force_flux_atm(dum_ia,i,j) = force_flux_atm(dum_ia,i,j)*force_flux_atm_sig_x(dum_ia)*loc_force_flux_atm_rtot
             END SELECT
          END DO
       END DO
    END IF
  END SUBROUTINE sub_update_force_flux_atm

  ! *** update sediment tracer flux forcing function value ***
  ! <<< GENERIC >>>
  SUBROUTINE sub_update_force_flux_sed(dum_t,dum_is)
    ! dummy arguments
    REAL,INTENT(in)::dum_t
    INTEGER,INTENT(in)::dum_is
    ! local variables
    INTEGER::i,j
    REAL::loc_x
    REAL::loc_force_flux_sed_tot
    REAL::loc_force_flux_sed_rtot
    REAL::loc_force_flux_sed
    real::loc_tot,loc_standard
    ! calculate new sediment tracer forcing time series values
    CALL sub_update_sig(dum_t,force_flux_sed_sig(dum_is,1,:),force_flux_sed_sig_i(dum_is,:),loc_x)
    force_flux_sed_sig_x(dum_is) = &
         & (1 - loc_x)*force_flux_sed_sig(dum_is,2,force_flux_sed_sig_i(dum_is,2)) + &
         & loc_x*force_flux_sed_sig(dum_is,2,force_flux_sed_sig_i(dum_is,1))
    ! update flux boundary conditions
    ! NOTE: flux forcings are in units of mol yr-1
    DO i=1,n_imax
       DO j=1,n_jmax
          loc_force_flux_sed = &
               & force_flux_sed_I(dum_is,i,j) + &
               & force_flux_sed_sig_x(dum_is)*(force_flux_sed_II(dum_is,i,j) - force_flux_sed_I(dum_is,i,j))
          SELECT CASE (sed_type(dum_is))
          CASE (par_sed_type_bio,par_sed_type_det,par_sed_type_POM,par_sed_type_CaCO3,par_sed_type_opal) !bio=1,det=2,POM=3,CaCO3=4,opal=5 so det_fe=6 (as well as all the 
             !                                                                                                                  other scavenged species) is not done here
             force_flux_sed(dum_is,i,j) = loc_force_flux_sed
          case (11:20)
             loc_tot  = force_flux_sed(sed_dep(dum_is),i,j)
             loc_standard = const_standards(sed_type(dum_is))
             force_flux_sed(dum_is,i,j) = fun_calc_isotope_fraction(loc_force_flux_sed,loc_standard)*loc_tot
          END SELECT
       END DO
    END DO
    ! normalize flux forcings (if selected) so that the total flux is equal to the magnitude (at the current time step) 
    ! defined in the forcing signal file
    IF (force_flux_sed_scale(dum_is)) THEN
       loc_force_flux_sed_tot = SUM(force_flux_sed(dum_is,:,:))
       if (abs(loc_force_flux_sed_tot) > const_real_nullsmall) then
          loc_force_flux_sed_rtot = 1.0/loc_force_flux_sed_tot
       else
          loc_force_flux_sed_rtot = 0.0
       end if
       DO i=1,n_imax
          DO j=1,n_jmax
             SELECT CASE (sed_type(dum_is))
             CASE (par_sed_type_bio)
                force_flux_sed(dum_is,i,j) = force_flux_sed(dum_is,i,j)*force_flux_sed_sig_x(dum_is)*loc_force_flux_sed_rtot
             end SELECT
          END DO
       END DO
    END IF
  END SUBROUTINE sub_update_force_flux_sed


! ********************************
! *** INVENTORY AUDIT ROUTINES ***
! ********************************


  ! *** calculate ocean tracer inventory ***
  ! <<< CONTAINS NON-GENERIC ALGORITHMS >>>
  FUNCTION fun_calc_ocn_tot()
    ! result variable
    REAL,dimension(0:n_ocn)::fun_calc_ocn_tot
    ! local variables
    INTEGER::l,i,j,k,io,is
    real,dimension(0:n_sed,n_maxi,n_maxj,n_maxk)::loc_bio_part
    real,dimension(0:n_ocn,n_maxi,n_maxj,n_maxk)::loc_bio_part_ocn
    real,dimension(0:n_ocn,n_maxi,n_maxj,n_maxk)::loc_ocn
    real,dimension(0:n_ocn,n_maxi,n_maxj,n_maxk)::loc_ocn_tot
    ! set local variables
    loc_bio_part(:,:,:,:)     = 0.0
    loc_bio_part_ocn(:,:,:,:) = 0.0
    loc_ocn(:,:,:,:)          = 0.0
    loc_ocn_tot(:,:,:,:)      = 0.0
    ! set default result
    fun_calc_ocn_tot(:) = 0.0
    ! convert particulate sediment and dissolved organic matter tracer concentrations to (dissolved) tracers
    ! NOTE: the 'old' algorithm for converting between tracers is retained even though it is much more numerically expensive
    !       as it is probably more robust to have the audit checking carried out by an independent method
    DO i=1,n_imax
       DO j=1,n_jmax
          DO k=goldstein_k1(i,j),n_kmax
             loc_ocn(:,i,j,k) = ocn(:,i,j,k)
             loc_bio_part(:,i,j,k) = bio_part(:,i,j,k)
             DO l=3,n_iomax
                io = conv_iselected_io(l)
                is = maxval(maxloc(abs(conv_DOM_POM(:,io))))-1
                if (is /= 0) then
                   loc_bio_part(is,i,j,k) = loc_bio_part(is,i,j,k) + loc_ocn(io,i,j,k)
                   loc_ocn(io,i,j,k) = 0.0
                end if
             end do
             loc_bio_part_ocn(:,i,j,k) = matmul(conv_sed_ocn(:,:),loc_bio_part(:,i,j,k))
          end do
       end do
    END DO
    ! determine ocean tracer inventory (mol)
    DO i=1,n_imax
       DO j=1,n_jmax
          DO k=goldstein_k1(i,j),n_kmax
             loc_ocn(:,i,j,k) = fun_audit_combinetracer(loc_ocn(:,i,j,k))
          end do
       end do
    END DO
!!$    ! NOTE: combine transformable tracer pairs when testing for drift;
!!$    !       NO3 + N2 (N2 -> NO3 during nitrogen fixation, and NO3 -> N2 during denitrification))
!!$    !       CO2 + CH4
!!$    !       SO4 + H2S
!!$    ! NOTE: ignore O2 component in oxidizing CH4->CO2 for now ...
!!$    loc_ocn(io_DIC,:,:,:) = loc_ocn(io_DIC,:,:,:) + loc_ocn(io_CH4,:,:,:)
!!$    loc_ocn(io_CH4,:,:,:) = 0.0
!!$    loc_ocn(io_O2,:,:,:) = loc_ocn(io_O2,:,:,:) + 1.5*loc_ocn(io_NO3,:,:,:) + 2.0*loc_ocn(io_SO4,:,:,:)
!!$    loc_ocn(io_NO3,:,:,:) = loc_ocn(io_NO3,:,:,:) + 2.0*loc_ocn(io_N2,:,:,:) + 2.0*loc_ocn(io_N2O,:,:,:)
!!$    loc_ocn(io_N2,:,:,:) = 0.0
!!$    loc_ocn(io_N2O,:,:,:) = 0.0
!!$    loc_ocn(io_SO4,:,:,:) = loc_ocn(io_SO4,:,:,:) + loc_ocn(io_H2S,:,:,:)
!!$    loc_ocn(io_H2S,:,:,:) = 0.0
!!$    ! ***********************************************
!!$    ! *** <INSERT CODE FOR FURTHER SPECIAL CASES> ***
!!$    ! ***********************************************
    DO l=3,n_iomax
       io = conv_iselected_io(l)
       loc_ocn_tot(io,:,:,:) = loc_ocn(io,:,:,:) + loc_bio_part_ocn(io,:,:,:)
       fun_calc_ocn_tot(io) = sum(loc_ocn_tot(io,:,:,:)*phys_ocn(ipo_M,:,:,:))
    end do
  END function fun_calc_ocn_tot


  ! *** carry out updated tracer audit ***
  ! <<< GENERIC FOR n_ocn TRACERS >>>
  SUBROUTINE sub_audit_update()
    ! local variables
    INTEGER::l,io
    REAL,dimension(0:n_ocn)::loc_audit_ocn_relerr
    ! set local variables
    loc_audit_ocn_relerr(:)   = 0.0
    ! calculate inventory drift
    audit_ocn_new(:) = fun_calc_ocn_tot()
    ! adjust ocean tracer inventory change (audit_ocn_delta) to combine different forms of the same element
    audit_ocn_delta(:) = fun_audit_combinetracer(audit_ocn_delta(:))
    ! 
    DO l=3,n_iomax
       io = conv_iselected_io(l)
       if (abs(audit_ocn_new(io)) > const_real_nullsmall) then
          loc_audit_ocn_relerr(io) = (audit_ocn_new(io) - (audit_ocn_old(io) + audit_ocn_delta(io)))/audit_ocn_new(io)
       else
          loc_audit_ocn_relerr(io) = 0.0
       end if
    end DO
    ! compare current (ocean) tracer inventory with estimate since last audit,
    ! upate maximum error encountered, and report error if relative change exceeds pre-defined threshold
    ! NOTE: do not report 14C (ocn 'type' 12), because it decays and is 'lost' in mass terms ...
    DO l=3,n_iomax
       io = conv_iselected_io(l)
       SELECT CASE (ocn_type(io))
       CASE (1,11,13,14)
          IF (ABS(loc_audit_ocn_relerr(io)) > par_misc_audit_relerr) THEN
             CALL sub_report_error('biogem_box','audit_update', &
                  & '(ocean) tracer inventory drift: new(top)/old(middle)/expected(bottom)): '//TRIM(string_ocn(io)), &
                  & 'n/a', &
                  & (/ &
                  &   audit_ocn_new(io), &
                  &   audit_ocn_old(io), &
                  &   (audit_ocn_old(io) + audit_ocn_delta(io)) &
                  & /), &
                  & opt_misc(iopt_misc_audit_fatal))
          ENDIF
          audit_ocn_old(io) = audit_ocn_new(io)
          audit_ocn_delta(io) = 0.0
       end SELECT
    END DO
  END SUBROUTINE sub_audit_update


  ! *** combine different forms of the same element ***
  function fun_audit_combinetracer(dum_ocn)
    ! result variable
    real,dimension(0:n_ocn)::fun_audit_combinetracer
    ! dummy arguments
    REAL,dimension(0:n_ocn),INTENT(in)::dum_ocn
    ! initialze result variable
    fun_audit_combinetracer(:) = dum_ocn(:)
    ! adjust ocean tracer inventory change (audit_ocn_delta) to combine different forms of the same element
    ! NOTE: combine transformable tracer pairs when testing for drift;
    !       NO3 + N2 (N2 -> NO3 during nitrogen fixation, and NO3 -> N2 during denitrification))
    !       CO2 + CH4
    !       SO4 + H2S
    ! NOTE: ignore O2 component in oxidizing CH4->CO2 for now ...
    fun_audit_combinetracer(io_DIC) = fun_audit_combinetracer(io_DIC) + dum_ocn(io_CH4)
    fun_audit_combinetracer(io_CH4) = 0.0
    fun_audit_combinetracer(io_O2)  = fun_audit_combinetracer(io_O2)  + 1.5*dum_ocn(io_NO3) + 2.0*dum_ocn(io_SO4)
    fun_audit_combinetracer(io_NO3) = fun_audit_combinetracer(io_NO3) + 2.0*dum_ocn(io_N2)  + 2.0*dum_ocn(io_N2O)
    fun_audit_combinetracer(io_N2)  = 0.0
    fun_audit_combinetracer(io_N2O) = 0.0
    fun_audit_combinetracer(io_SO4) = fun_audit_combinetracer(io_SO4) + dum_ocn(io_H2S)
    fun_audit_combinetracer(io_H2S) = 0.0
!COPY from Andy's codes Sun 2009/03/03
    if (ocn_select(io_Fe)) then
       fun_audit_combinetracer(io_Fe)  = fun_audit_combinetracer(io_Fe)  + dum_ocn(io_FeL)
    end if
    if (ocn_select(io_Ligand)) then
       fun_audit_combinetracer(io_Ligand)   = fun_audit_combinetracer(io_Ligand)   + dum_ocn(io_FeL)
    end if
    fun_audit_combinetracer(io_FeL) = 0.0

    ! ***********************************************
    ! *** <INSERT CODE FOR FURTHER SPECIAL CASES> ***
    ! ***********************************************
  end function fun_audit_combinetracer


  ! ******************************
  ! *** MISCELLANEOUS ROUTINES ***
  ! ******************************


  ! *** copy tracer array ***
  ! <<< GENERIC FOR n_ocn TRACERS >>>
  SUBROUTINE sub_biogem_copy_ocntots(dum_ts,dum_ts1)
    USE biogem_lib
    ! dummy arguments
    REAL,DIMENSION(1:n_maxl,1:n_maxi,1:n_maxj,1:n_maxk),INTENT(inout)::dum_ts
    REAL,DIMENSION(1:n_maxl,1:n_maxi,1:n_maxj,1:n_maxk),INTENT(inout)::dum_ts1
    ! local variables
    INTEGER::i,j,k,l,io
    real::loc_ocn_mean_S,loc_ocn_tot_M
#ifdef Snorm                                  !yes
    ! [SALINITY NORMALIZED SCHEME]
    ! calculate total ocean mass
    loc_ocn_tot_M = sum(phys_ocn(ipo_M,:,:,:))
    ! calculate mean ocean salinity
    loc_ocn_mean_S = SUM(ocn(io_S,:,:,:)*phys_ocn(ipo_M,:,:,:))/loc_ocn_tot_M
    ! copy GOLDSTEIn <ts> array values from the relevant <ocn> array of BioGeM
    ! NOTE: leave T (index 1) and S (index 2) well alone ;-)
    ! NOTE: no offset (array: <tstoocn_offset()>) required for biogeochem-only tracers
    ! NOTE: normalize by relative salinity deviation from ocean mean
    DO i=1,n_imax
       DO j=1,n_jmax
          DO k=goldstein_k1(i,j),n_kmax
             DO l=3,n_iomax
                io = conv_iselected_io(l)
                dum_ts(l,i,j,k)  = ocn(io,i,j,k)*(loc_ocn_mean_S/ocn(io_S,i,j,k))
                dum_ts1(l,i,j,k) = dum_ts(l,i,j,k)
             end do
          END DO
       END DO
    END DO
#else
    ! [NON-SALINITY NORMALIZED SCHEME]
    DO i=1,n_imax
       DO j=1,n_jmax
          DO k=goldstein_k1(i,j),n_kmax
             DO l=3,n_iomax
                io = conv_iselected_io(l)
                dum_ts(l,i,j,k)  = ocn(io,i,j,k)
                dum_ts1(l,i,j,k) = dum_ts(l,i,j,k)
             end do
          END DO
       END DO
    END DO
#endif
  END SUBROUTINE sub_biogem_copy_ocntots
  
  
  ! *** copy tracer array ***
  ! <<< GENERIC FOR n_ocn TRACERS >>>
  SUBROUTINE sub_biogem_copy_tstoocn(dum_ts)
    USE biogem_lib
    ! dummy arguments
    REAL,DIMENSION(1:n_maxl,1:n_maxi,1:n_maxj,1:n_maxk),INTENT(in)::dum_ts
    ! local variables
    INTEGER::i,j,k
    ! copy BioGeM <ocn> array values from the relevant <ts> (or <ts1>) array of GOLDSTEIn
    ! NOTE: restrict to T and S
    DO i=1,n_imax
       DO j=1,n_jmax
          DO k=1,n_kmax
             IF (k >= goldstein_k1(i,j)) THEN
                ocn(io_T,i,j,k) = dum_ts(1,i,j,k) + tstoocn_offset(1)
                ocn(io_S,i,j,k) = dum_ts(2,i,j,k) + tstoocn_offset(2)
             END IF
          END DO
       END DO
    END DO
  END SUBROUTINE sub_biogem_copy_tstoocn


  ! *** copy of GOLDSTEIn overturning streamfunction calculation ***
  ! <<< CONTAINS NON-GENERIC ALGORITHMS >>>
  SUBROUTINE sub_calc_psi(dum_u,dum_opsi,dum_opsia,dum_opsip,dum_zpsi)
    ! dummy arguments
    REAL,INTENT(in),DIMENSION(3,1:n_maxi,1:n_maxj,1:n_maxk)::dum_u
    REAL,INTENT(inout),DIMENSION(0:n_maxj,0:n_maxk)::dum_opsi,dum_opsia,dum_opsip,dum_zpsi
    ! local variables
    INTEGER::i,j,k
    REAL::loc_ominp,loc_omaxp
    REAL::loc_omina,loc_omaxa
    REAL,DIMENSION(n_maxj,n_maxk)::loc_ou,loc_zu
    REAL,DIMENSION(0:n_maxj,0:n_maxk)::loc_opsi,loc_opsia,loc_opsip,loc_zpsi
    ! c Calculate meridional overturning streamfunction opsi on C grid only
    loc_opsi(:,:) = 0.0
    loc_opsia(:,:) = 0.0
    loc_opsip(:,:) = 0.0
    DO j=1,n_jmax-1
       DO k=1,n_kmax-1
          loc_ou(j,k) = 0.0
          DO i=1,n_imax
             loc_ou(j,k) = loc_ou(j,k) + goldstein_cv(j)*dum_u(2,i,j,k)*goldstein_dphi
          END DO
          loc_opsi(j,k) = loc_opsi(j,k-1) - goldstein_dz(k)*loc_ou(j,k)
       END DO
    END DO
    ! c Pacific and Atlantic overturning streamfunctions
    loc_ominp = 0.0
    loc_omaxp = 0.0
    DO j=goldstein_jsf+1,n_jmax-1
       DO k=1,n_kmax-1
          loc_ou(j,k) = 0.0
          DO i=goldstein_ips(j),goldstein_ipf(j)
             loc_ou(j,k) = loc_ou(j,k) + goldstein_cv(j)*dum_u(2,i,j,k)*goldstein_dphi
          ENDDO
          loc_opsip(j,k) = loc_opsip(j,k-1) - goldstein_dz(k)*loc_ou(j,k)
          IF((loc_opsip(j,k) < loc_ominp) .AND. (z_at_k(k) >= 300.0)) loc_ominp = loc_opsip(j,k)
          IF((loc_opsip(j,k) > loc_omaxp) .AND. (z_at_k(k) >= 300.0)) loc_omaxp = loc_opsip(j,k)
        ENDDO
    ENDDO
    loc_omina = 0.0
    loc_omaxa = 0.0
    DO j=goldstein_jsf+1,n_jmax-1
       DO k=1,n_kmax-1
          loc_ou(j,k) = 0.0
          DO i=goldstein_ias(j),goldstein_iaf(j)
             loc_ou(j,k) = loc_ou(j,k) + goldstein_cv(j)*dum_u(2,i,j,k)*goldstein_dphi
          ENDDO
          loc_opsia(j,k) = loc_opsia(j,k-1) - goldstein_dz(k)*loc_ou(j,k)
          IF((loc_opsia(j,k) < loc_omina) .AND. (z_at_k(k) >= 300.0)) loc_omina = loc_opsia(j,k)
          IF((loc_opsia(j,k) > loc_omaxa) .AND. (z_at_k(k) >= 300.0)) loc_omaxa = loc_opsia(j,k)
       ENDDO
    ENDDO
    loc_zpsi(:,:) = 0.0
    DO i=1,n_imax-1
       DO k=1,n_kmax-1
          loc_zu(i,k) = 0
          DO j=1,n_jmax
             loc_zu(i,k) = loc_zu(i,k) + dum_u(1,i,j,k)/goldstein_c(j)*goldstein_ds
          ENDDO
          loc_zpsi(i,k) = loc_zpsi(i,k-1) - goldstein_dz(k)*loc_zu(i,k)
       ENDDO
    ENDDO
    ! set results arrays
    dum_opsi(1:n_maxj,1:n_maxk)  = loc_opsi(1:n_maxj,1:n_maxk)
    dum_opsia(1:n_maxj,1:n_maxk) = loc_opsia(1:n_maxj,1:n_maxk)
    dum_opsip(1:n_maxj,1:n_maxk) = loc_opsip(1:n_maxj,1:n_maxk)
    dum_zpsi(1:n_maxj,1:n_maxk)  = loc_zpsi(1:n_maxj,1:n_maxk)
  END SUBROUTINE sub_calc_psi

  ! km 29-Mar-07; flux (Tg-N/yr) in equal proportions in forms of DIN and DON (literature)
  !   BUT put all in NO3 and none in DON, because N-fixation to compensate for NO3 river flux
  !    is easy but to compensate for DON river flux is difficult; also DON->DIN conversion is "fast"
  !   ALSO, no need to multiply forcing flux (Tg-N/yr) by dt; this subroutine is called once ayear
  !   LAST, the addition of NO3- adds negative charge, so reduces ALK
  SUBROUTINE sub_anthN_river(loc_anthn_flux) 

    USE biogem_lib 

    integer::i,j,l
    real,dimension(indexr)::loc_anthn_flux
    real::loc_delta,loc_standard,loc_frac

    do i = 1, n_maxi
       do j = 1, n_maxj
          if(riv_f(i,j) > 0.)then
             if(anth_area(i,j).eq.1.or.anth_area(i,j).eq.2.)then      ! Atlantic
                l = 2   
             elseif(anth_area(i,j).eq.3.or.anth_area(i,j).eq.4.)then  ! Pacific
                l = 5 
             elseif(anth_area(i,j).eq.5.or.anth_area(i,j).eq.6.)then  ! Indian
                l = 3
             elseif(anth_area(i,j).eq.7.)then                         ! Mediterranean
                l = 4
             elseif(anth_area(i,j).eq.8.)then                         ! Arctic
                l = 1
             endif
                                                        
             flux_rivN(i,j) = riv_f(i,j)*loc_anthn_flux(l)*1.e12/14.*phys_ocn(ipo_rM,i,j,n_maxk)   ! mol-N/kg
             ocn(io_NO3,i,j,n_maxk) = ocn(io_NO3,i,j,n_maxk) + flux_rivN(i,j)                      ! flux is /yr
#ifndef noalk
             ocn(io_ALK,i,j,n_maxk) = ocn(io_ALK,i,j,n_maxk) - flux_rivN(i,j)                      ! flux is /yr             
#endif
!add po4, too, so it won't be limiting
             ocn(io_PO4,i,j,n_maxk) = ocn(io_PO4,i,j,n_maxk) + flux_rivN(i,j)/par_bio_red_POP_PON  ! flux is /yr
             if(ocn_select(io_NO3_15N))then 
                loc_delta = 0.d0                                      ! assume 0 permil for river N15, ok?
                loc_standard = const_standards(ocn_type(io_NO3_15N)) 
                loc_frac     = fun_calc_isotope_fraction(loc_delta,loc_standard) 
                flux_riv15N(i,j) = loc_frac * flux_rivN(i,j)                           ! mol/kg
 
                ocn(io_NO3_15N,i,j,n_maxk) = ocn(io_NO3_15N,i,j,n_maxk) + flux_riv15N(i,j)
             endif 
          endif
       enddo
    enddo

  END SUBROUTINE sub_anthN_river
  

  ! km Alkalinity river flux distributed uniformly to coastal grids w/ runoff>0
  !    There is DIC associated with ALK (assumes HCO3-; so the same moles as ALK)
  SUBROUTINE sub_anthA_river(loc_antha_flux) 

    USE biogem_lib 

    integer::i,j
    real::loc_antha_flux            ! 1e13 moles-Alk/year

    do i = 1, n_maxi
       do j = 1, n_maxj
          if(riv_f(i,j) > 0.)then
             flux_rivA(i,j) = loc_antha_flux*1.e13/ngrid_tot*phys_ocn(ipo_rM,i,j,n_maxk)   ! mol-A/kg

             ocn(io_ALK,i,j,n_maxk) = ocn(io_ALK,i,j,n_maxk) + flux_rivA(i,j)                 ! flux is /yr
             ocn(io_DIC,i,j,n_maxk) = ocn(io_DIC,i,j,n_maxk) + flux_rivA(i,j)                 ! flux is /yr
          endif
       enddo
    enddo

  END SUBROUTINE sub_anthA_river


  ! km Carbon (DOC+POC) river flux distributed uniformly to coastal grids w/ runoff>0
  !    Kempe says POC is erosional and mostly inert; so forget this component (45%)
  !          remaining 55% is DOC - add this; so assume that this fraction remains the same
  SUBROUTINE sub_anthC_river(loc_anthc_flux) 

    USE biogem_lib 

    integer::i,j
    real::loc_anthc_flux           ! 1e14 grams-C/year

    do i = 1, n_maxi
       do j = 1, n_maxj
          if(riv_f(i,j) > 0.)then
             flux_rivC(i,j) = loc_anthc_flux*1.e14/ngrid_tot/12.*phys_ocn(ipo_rM,i,j,n_maxk)   ! mol-C/kg
             ocn(io_DOM_C,i,j,n_maxk) = ocn(io_DOM_C,i,j,n_maxk) + flux_rivC(i,j)*0.55                ! flux is /yr
          endif
       enddo
    enddo

  END SUBROUTINE sub_anthC_river


END MODULE biogem_box


  !        *******
  !    ***************
  !  ********************
  ! **********************
  ! *** CODE GRAVEYARD ***
  ! **********************
  ! **********************
  ! **      **  **      **
  ! **  **  **  **  **  **
  ! **  **  **  **  **  **
  ! **  *  ***  **      **
  ! **  **  **  **  ******
  ! **  **  **  **  ******
  ! **  **  **  **  ******
  ! **********************
  ! **********************
  ! **********************
  ! **********************
