123     &                        LBi, UBi, LBj, UBj, UBk, UBt,             &
 
  124     &                        IminS, ImaxS, JminS, JmaxS,               &
 
  128# if defined WET_DRY && defined DIAGNOSTICS_BIO
 
  134#ifdef DIAGNOSTICS_BIO
 
  149      integer, 
intent(in) :: ng, tile
 
  150      integer, 
intent(in) :: LBi, UBi, LBj, UBj, UBk, UBt
 
  151      integer, 
intent(in) :: IminS, ImaxS, JminS, JmaxS
 
  152      integer, 
intent(in) :: nstp, nnew
 
  156      real(r8), 
intent(in) :: rmask(LBi:,LBj:)
 
  157#  if defined WET_DRY && defined DIAGNOSTICS_BIO 
  158      real(r8), 
intent(in) :: rmask_full(LBi:,LBj:)
 
  161      real(r8), 
intent(in) :: Hz(LBi:,LBj:,:)
 
  162      real(r8), 
intent(in) :: z_r(LBi:,LBj:,:)
 
  163      real(r8), 
intent(in) :: z_w(LBi:,LBj:,0:)
 
  164      real(r8), 
intent(in) :: SpecIr(LBi:,LBj:,:)
 
  165      real(r8), 
intent(in) :: avcos(LBi:,LBj:,:)
 
  166# ifdef DIAGNOSTICS_BIO 
  167      real(r8), 
intent(inout) :: DiaBio3d(LBi:,LBj:,:,:)
 
  168      real(r8), 
intent(inout) :: DiaBio4d(LBi:,LBj:,:,:,:)
 
  170      real(r8), 
intent(inout) :: t(LBi:,LBj:,:,:,:)
 
  173      real(r8), 
intent(in) :: rmask(LBi:UBi,LBj:UBj)
 
  174#  if defined WET_DRY && defined DIAGNOSTICS_BIO 
  175      real(r8), 
intent(in) :: rmask_full(LBi:UBi,LBj:UBj)
 
  178      real(r8), 
intent(in) :: Hz(LBi:UBi,LBj:UBj,UBk)
 
  179      real(r8), 
intent(in) :: z_r(LBi:UBi,LBj:UBj,UBk)
 
  180      real(r8), 
intent(in) :: z_w(LBi:UBi,LBj:UBj,0:UBk)
 
  181      real(r8), 
intent(in) :: SpecIr(LBi:UBi,LBj:UBj,NBands)
 
  182      real(r8), 
intent(in) :: avcos(LBi:UBi,LBj:UBj,NBands)
 
  183# ifdef DIAGNOSTICS_BIO 
  184      real(r8), 
intent(inout) :: DiaBio3d(LBi:UBi,LBj:UBj,              &
 
  186      real(r8), 
intent(inout) :: DiaBio4d(LBi:UBi,LBj:UBj,N(ng),        &
 
  189      real(r8), 
intent(inout) :: t(LBi:UBi,LBj:UBj,UBk,3,UBt)
 
  194      integer, 
parameter :: Msink = 30
 
  196      integer :: i, j, k, ks
 
  197      integer :: Iter, Tindex, ic, isink, ibio, id, itrc, ivar
 
  198      integer :: ibac, iband, idom, ifec, iphy, ipig
 
  201      integer, 
dimension(Msink) :: idsink
 
  202      integer, 
dimension(IminS:ImaxS,N(ng)) :: ksource
 
  204      real(r8), 
parameter :: MinVal = 0.0_r8
 
  206      real(r8) :: FV1, FV2, FV3, FV4, FV5, FV6, FV7, dtbio
 
  207      real(r8) :: DOC_lab, Ed_tot, Nup_max, aph442, aPHYN_wa
 
  208      real(r8) :: avgcos_min, par_b, par_s, photo_DIC, photo_DOC
 
  209      real(r8) :: photo_decay, slope_AC, tChl, theta_m, total_photo
 
  210      real(r8) :: WLE, factint
 
  213      real(r8) :: N_quota, RelDOC1, RelDON1, RelDOP1, RelFe
 
  214      real(r8) :: cff, cff1, cff2, cffL, cffR, cu, dltL, dltR
 
  215#ifdef DIAGNOSTICS_BIO 
  219      real(r8), 
dimension(Msink) :: Wbio
 
  221      real(r8), 
dimension(4) :: Bac_G
 
  223      real(r8), 
dimension(NBands) :: dATT_sum
 
  225      real(r8), 
dimension(N(ng),NBands) :: avgcos, dATT
 
  226      real(r8), 
dimension(N(ng),NBands) :: specir_d
 
  227      real(r8), 
dimension(N(ng),NBands) :: tot_ab, tot_b, tot_s
 
  229      real(r8), 
dimension(0:N(ng),NBands) :: specir_w
 
  231      real(r8), 
dimension(N(ng),Nphy) :: C2CHL, C2CHL_w
 
  232      real(r8), 
dimension(N(ng),Nphy) :: Gt_fl, Gt_ll, Gt_nl
 
  233      real(r8), 
dimension(N(ng),Nphy) :: Gt_sl, Gt_pl
 
  234      real(r8), 
dimension(N(ng),Nphy) :: alfa
 
  235      real(r8), 
dimension(N(ng),Nphy) :: pac_eff
 
  237      real(r8), 
dimension(N(ng),Nphy,Npig) :: Pigs_w
 
  239      integer, 
dimension(IminS:ImaxS) :: Keuphotic
 
  241      real(r8), 
dimension(IminS:ImaxS,N(ng)) :: E0_nz
 
  242      real(r8), 
dimension(IminS:ImaxS,N(ng)) :: Ed_nz
 
  243      real(r8), 
dimension(IminS:ImaxS,N(ng)) :: DOC_frac
 
  244      real(r8), 
dimension(IminS:ImaxS,N(ng)) :: NitrBAC
 
  245      real(r8), 
dimension(IminS:ImaxS,N(ng)) :: NH4toNO3
 
  246      real(r8), 
dimension(IminS:ImaxS,N(ng)) :: NtoNBAC
 
  247      real(r8), 
dimension(IminS:ImaxS,N(ng)) :: NtoPBAC
 
  248      real(r8), 
dimension(IminS:ImaxS,N(ng)) :: NtoFeBAC
 
  249      real(r8), 
dimension(IminS:ImaxS,N(ng)) :: totDOC_d
 
  250      real(r8), 
dimension(IminS:ImaxS,N(ng)) :: totDON_d
 
  251      real(r8), 
dimension(IminS:ImaxS,N(ng)) :: totDOP_d
 
  252      real(r8), 
dimension(IminS:ImaxS,N(ng)) :: totFe_d
 
  253      real(r8), 
dimension(IminS:ImaxS,N(ng)) :: totNH4_d
 
  254      real(r8), 
dimension(IminS:ImaxS,N(ng)) :: totNO3_d
 
  255      real(r8), 
dimension(IminS:ImaxS,N(ng)) :: totPO4_d
 
  256      real(r8), 
dimension(IminS:ImaxS,N(ng)) :: totSiO_d
 
  258      real(r8), 
dimension(IminS:ImaxS,N(ng),Nbac) :: GtBAC
 
  259      real(r8), 
dimension(IminS:ImaxS,N(ng),Nbac) :: NupDOC_ba
 
  260      real(r8), 
dimension(IminS:ImaxS,N(ng),Nbac) :: NupDON_ba
 
  261      real(r8), 
dimension(IminS:ImaxS,N(ng),Nbac) :: NupDOP_ba
 
  262      real(r8), 
dimension(IminS:ImaxS,N(ng),Nbac) :: NupFe_ba
 
  263      real(r8), 
dimension(IminS:ImaxS,N(ng),Nbac) :: NupNH4_ba
 
  264      real(r8), 
dimension(IminS:ImaxS,N(ng),Nbac) :: NupPO4_ba
 
  266      real(r8), 
dimension(IminS:ImaxS,N(ng),Nphy) :: C2fALG
 
  267      real(r8), 
dimension(IminS:ImaxS,N(ng),Nphy) :: C2nALG
 
  268      real(r8), 
dimension(IminS:ImaxS,N(ng),Nphy) :: C2pALG
 
  269      real(r8), 
dimension(IminS:ImaxS,N(ng),Nphy) :: C2sALG
 
  270      real(r8), 
dimension(IminS:ImaxS,N(ng),Nphy) :: GtALG
 
  271      real(r8), 
dimension(IminS:ImaxS,N(ng),Nphy) :: GtALG_r
 
  272      real(r8), 
dimension(IminS:ImaxS,N(ng),Nphy) :: NupDOP
 
  273      real(r8), 
dimension(IminS:ImaxS,N(ng),Nphy) :: NupDON
 
  274      real(r8), 
dimension(IminS:ImaxS,N(ng),Nphy) :: NupFe
 
  275      real(r8), 
dimension(IminS:ImaxS,N(ng),Nphy) :: NupNH4
 
  276      real(r8), 
dimension(IminS:ImaxS,N(ng),Nphy) :: NupNO3
 
  277      real(r8), 
dimension(IminS:ImaxS,N(ng),Nphy) :: NupPO4
 
  278      real(r8), 
dimension(IminS:ImaxS,N(ng),Nphy) :: NupSiO
 
  279      real(r8), 
dimension(IminS:ImaxS,N(ng),Nphy) :: graz_act
 
  280      real(r8), 
dimension(IminS:ImaxS,N(ng),Nphy) :: mu_bar_f
 
  281      real(r8), 
dimension(IminS:ImaxS,N(ng),Nphy) :: mu_bar_n
 
  282      real(r8), 
dimension(IminS:ImaxS,N(ng),Nphy) :: mu_bar_p
 
  283      real(r8), 
dimension(IminS:ImaxS,N(ng),Nphy) :: mu_bar_s
 
  284      real(r8), 
dimension(IminS:ImaxS,N(ng),Nphy) :: refuge
 
  286      real(r8), 
dimension(IminS:ImaxS,N(ng),Nfec) :: Regen_C
 
  287      real(r8), 
dimension(IminS:ImaxS,N(ng),Nfec) :: Regen_F
 
  288      real(r8), 
dimension(IminS:ImaxS,N(ng),Nfec) :: Regen_N
 
  289      real(r8), 
dimension(IminS:ImaxS,N(ng),Nfec) :: Regen_P
 
  290      real(r8), 
dimension(IminS:ImaxS,N(ng),Nfec) :: Regen_S
 
  292      real(r8), 
dimension(IminS:ImaxS,N(ng),NBands) :: specir_scal
 
  293      real(r8), 
dimension(IminS:ImaxS,N(ng),Nphy,NBands) :: aPHYN_al
 
  294      real(r8), 
dimension(IminS:ImaxS,N(ng),Nphy,NBands) :: aPHYN_at
 
  295      real(r8), 
dimension(IminS:ImaxS,N(ng),NBands) :: aDET
 
  296      real(r8), 
dimension(IminS:ImaxS,N(ng),NBands) :: aCDC
 
  297      real(r8), 
dimension(IminS:ImaxS,N(ng),NBands) :: b_phy
 
  298      real(r8), 
dimension(IminS:ImaxS,N(ng),NBands) :: s_phy
 
  299      real(r8), 
dimension(IminS:ImaxS,N(ng),NBands) :: b_tot
 
  300      real(r8), 
dimension(IminS:ImaxS,N(ng),NBands) :: s_tot
 
  302      real(r8), 
dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio
 
  303      real(r8), 
dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio_old
 
  304      real(r8), 
dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio_new
 
  306      real(r8), 
dimension(IminS:ImaxS,0:N(ng)) :: FC
 
  307      real(r8), 
dimension(IminS:ImaxS,N(ng)) :: Hz_inv
 
  308      real(r8), 
dimension(IminS:ImaxS,N(ng)) :: Hz_inv2
 
  309      real(r8), 
dimension(IminS:ImaxS,N(ng)) :: Hz_inv3
 
  310      real(r8), 
dimension(IminS:ImaxS,N(ng)) :: WL
 
  311      real(r8), 
dimension(IminS:ImaxS,N(ng)) :: WR
 
  312      real(r8), 
dimension(IminS:ImaxS,N(ng)) :: bL
 
  313      real(r8), 
dimension(IminS:ImaxS,N(ng)) :: bR
 
  314      real(r8), 
dimension(IminS:ImaxS,N(ng)) :: qc
 
  316#include "set_bounds.h" 
  317#ifdef DIAGNOSTICS_BIO 
  324     &     (mod(
iic(ng),
ndia(ng)).eq.1)).or.                            &
 
  331                diabio3d(i,j,k,ivar)=0.0_r8
 
  341                  diabio4d(i,j,k,iband,ivar)=0.0_r8
 
  357#ifdef DIAGNOSTICS_BIO 
  362      fiter=1.0_r8/real(
bioiter(ng),r8)
 
  370        idsink(ic)=
ifecn(ifec)
 
  373        idsink(ic)=
ifecc(ifec)
 
  376        idsink(ic)=
ifecp(ifec)
 
  379        idsink(ic)=
ifecs(ifec)
 
  382        idsink(ic)=
ifecf(ifec)
 
  388        idsink(ic)=
iphyn(iphy)
 
  391        idsink(ic)=
iphyc(iphy)
 
  394        idsink(ic)=
iphyp(iphy)
 
  396        IF (
iphys(iphy).ne.0) 
THEN 
  398          idsink(ic)=
iphys(iphy)
 
  402        idsink(ic)=
iphyf(iphy)
 
  412      j_loop : 
DO j=jstr,jend
 
  415            hz_inv(i,k)=1.0_r8/hz(i,j,k)
 
  420            hz_inv2(i,k)=1.0_r8/(hz(i,j,k)+hz(i,j,k+1))
 
  425            hz_inv3(i,k)=1.0_r8/(hz(i,j,k-1)+hz(i,j,k)+hz(i,j,k+1))
 
  443              bio(i,k,itrc)=max(minval,t(i,j,k,nstp,itrc))
 
  444              bio_old(i,k,itrc)=bio(i,k,itrc)
 
  453              bio_new(i,k,itrc)=0.0_r8
 
  491                regen_c(i,k,ifec)=
regcr(ifec,ng)*fv1
 
  492                regen_n(i,k,ifec)=
regnr(ifec,ng)*fv1
 
  493                regen_p(i,k,ifec)=
regpr(ifec,ng)*fv1
 
  494                regen_f(i,k,ifec)=
regfr(ifec,ng)*fv1
 
  495                regen_s(i,k,ifec)=
regsr(ifec,ng)*fv1
 
  512              fv1=
maxc2nalg(iphy,ng)*(1.0_r8+gtalg(i,k,iphy))
 
  513              mu_bar_n(i,k,iphy)=gtalg(i,k,iphy)*                       &
 
  516                fv1=
maxc2sialg(iphy,ng)*(1.0_r8+gtalg(i,k,iphy))
 
  517                mu_bar_s(i,k,iphy)=gtalg(i,k,iphy)*                     &
 
  523                fv1=
maxc2palg(iphy,ng)*(1.0_r8+gtalg(i,k,iphy))
 
  524                mu_bar_p(i,k,iphy)=gtalg(i,k,iphy)*                     &
 
  530                fv1=
maxc2fealg(iphy,ng)*(1.0_r8+gtalg(i,k,iphy))
 
  531                mu_bar_f(i,k,iphy)=gtalg(i,k,iphy)*                     &
 
  562              fv1=max(1.0_r8,(bio(i,k,
iphyc(iphy))/refuge(i,k,iphy)))
 
  563              graz_act(i,k,iphy)=
hsgrz(iphy,ng)*log(fv1)
 
  572        iter_loop : 
DO iter=1,
bioiter(ng)
 
  589                nupnh4(i,k,iphy)=0.0_r8
 
  590                nupno3(i,k,iphy)=0.0_r8
 
  591                nuppo4(i,k,iphy)=0.0_r8
 
  592                nupsio(i,k,iphy)=0.0_r8
 
  593                nupfe(i,k,iphy)=0.0_r8
 
  594                nupdon(i,k,iphy)=0.0_r8
 
  595                nupdop(i,k,iphy)=0.0_r8
 
  606                c2nalg(i,k,iphy)=0.0_r8
 
  607                IF (bio(i,k,
iphyn(iphy)).gt.0.0_r8) 
THEN 
  608                  c2nalg(i,k,iphy)=bio(i,k,
iphyc(iphy))/                &
 
  609     &                             bio(i,k,
iphyn(iphy))
 
  611                c2palg(i,k,iphy)=0.0_r8
 
  612                IF (bio(i,k,
iphyp(iphy)).gt.0.0_r8) 
THEN 
  613                  c2palg(i,k,iphy)=bio(i,k,
iphyc(iphy))/                &
 
  614     &                             bio(i,k,
iphyp(iphy))
 
  616                c2salg(i,k,iphy)=0.0_r8
 
  617                IF (
iphys(iphy).gt.0) 
THEN 
  618                  IF (bio(i,k,
iphys(iphy)).gt.0.0_r8) 
THEN 
  619                    c2salg(i,k,iphy)=bio(i,k,
iphyc(iphy))/              &
 
  620     &                               bio(i,k,
iphys(iphy))
 
  623                c2falg(i,k,iphy)=0.0_r8
 
  624                IF (bio(i,k,
iphyf(iphy)).gt.0.0_r8) 
THEN 
  625                  c2falg(i,k,iphy)=bio(i,k,
iphyc(iphy))/                &
 
  626     &                             bio(i,k,
iphyf(iphy))
 
  639            ed_nz(i,n(ng))=0.0_r8
 
  640            e0_nz(i,n(ng))=0.0_r8
 
  642            IF (specir(i,j,21).gt.
vsmall) 
THEN 
  648                datt_sum(iband)=0.0_r8
 
  650                  avgcos(k,iband)=0.0_r8
 
  652#ifdef DIAGNOSTICS_BIO 
  653                  adet(i,k,iband)=0.0_r8
 
  654                  acdc(i,k,iband)=0.0_r8
 
  655                  b_phy(i,k,iband)=0.0_r8
 
  656                  s_phy(i,k,iband)=0.0_r8
 
  657                  b_tot(i,k,iband)=0.0_r8
 
  658                  s_tot(i,k,iband)=0.0_r8
 
  663                    aphyn_at(i,k,iphy,iband)=0.0_r8
 
  664                    aphyn_al(i,k,iphy,iband)=0.0_r8
 
  674                ed_tot=ed_tot+specir(i,j,iband)*
dlam 
  675                avgcos(n(ng),iband)=avcos(i,j,iband)
 
  684                IF (ed_tot.ge.1.0_r8) 
THEN 
  688                    IF (bio(i,k,
iphyc(iphy)).gt.0.0_r8) 
THEN 
  690                      pac_eff(k,iphy)=1.0_r8
 
  693     &                      (bio(i,k,
iphyc(iphy))*12.0_r8)
 
  694                        pac_eff(k,iphy)=max(0.5_r8,                     &
 
  703                        IF (
ipigs(iphy,ipig).gt.0) 
THEN 
  705     &                           bio(i,k,
ipigs(iphy,ipig))*             &
 
  706     &                           apigs(ipig,iband)*pac_eff(k,iphy)
 
  721                        IF (
ipigs(iphy,ipig).gt.0) 
THEN 
  722                          aphyn_at(i,k,iphy,iband)=                     &
 
  723     &                                      aphyn_at(i,k,iphy,iband)+   &
 
  724     &                                      bio(i,k,
ipigs(iphy,ipig))*  &
 
  725     &                                      apigs(ipig,iband)*          &
 
  729                      tot_ab(k,iband)=tot_ab(k,iband)+                  &
 
  730     &                                aphyn_at(i,k,iphy,iband)
 
  731#ifdef DIAGNOSTICS_BIO 
  732                      diabio4d(i,j,k,iband,
idaphy)=diabio4d(i,j,k,iband,&
 
  734     &                                             aphyn_at(i,k,iphy,   &
 
  741                      IF (
ipigs(iphy,ipig).gt.0) 
THEN 
  742                        aphyn_al(i,k,iphy,iband)=                       &
 
  743     &                                    aphyn_at(i,k,iphy,iband)-     &
 
  744     &                                    bio(i,k,
ipigs(iphy,ipig))*    &
 
  745     &                                    apigs(ipig,iband)*            &
 
  752                    cff=aph442*exp(0.011_r8*                            &
 
  754     &                              (397.0_r8+real(iband,r8)*
dlam)))
 
  755                    tot_ab(k,iband)=tot_ab(k,iband)+cff
 
  756#ifdef DIAGNOSTICS_BIO 
  757                    adet(i,k,iband)=adet(i,k,iband)+cff
 
  758                    diabio4d(i,j,k,iband,
idadet)=diabio4d(i,j,k,iband,  &
 
  770     &                            adoc(
ilab,iband)+                     &
 
  773                    tot_ab(k,iband)=tot_ab(k,iband)+cff+awater(iband)
 
  774#ifdef DIAGNOSTICS_BIO 
  775                    acdc(i,k,iband)=acdc(i,k,iband)+cff
 
  776                    diabio4d(i,j,k,iband,
idacdc)=diabio4d(i,j,k,iband,  &
 
  787                    par_s=0.3_r8*(tchl**0.62_r8)       
 
  789                    IF (tchl.gt.0.0_r8) 
THEN 
  790                      par_b=par_s*(0.002_r8+                            &
 
  792     &                             (0.5_r8-0.25_r8*log10(tchl))*        &
 
  795                    par_b=max(par_b,0.0_r8)
 
  796#ifdef DIAGNOSTICS_BIO 
  797                    s_phy(i,k,iband)=s_phy(i,k,iband)+par_s
 
  798                    b_phy(i,k,iband)=b_phy(i,k,iband)+par_b
 
  799                    diabio4d(i,j,k,iband,
idsphy)=diabio4d(i,j,k,iband,  &
 
  802                    diabio4d(i,j,k,iband,
idbphy)=diabio4d(i,j,k,iband,  &
 
  810                    tot_s(k,iband)=bwater(iband)+par_s*
wavedp(iband)
 
  811#ifdef DIAGNOSTICS_BIO 
  812                    s_tot(i,k,iband)=s_tot(i,k,iband)+tot_s(k,iband)
 
  817                    tot_b(k,iband)=0.5_r8*bwater(iband)+par_b
 
  818#ifdef DIAGNOSTICS_BIO 
  819                    b_tot(i,k,iband)=b_tot(i,k,iband)+tot_b(k,iband)
 
  821                    diabio4d(i,j,k,iband,
idstot)=diabio4d(i,j,k,iband,  &
 
  824                    diabio4d(i,j,k,iband,
idbtot)=diabio4d(i,j,k,iband,  &
 
  835     &                   0.005_r8*acos(avgcos(k,iband))*
rad2deg 
  836                    cff2=4.18_r8*(1.0_r8-0.52_r8*                       &
 
  837     &                            exp(-10.8_r8*tot_ab(k,iband)))
 
  838                    datt(k,iband)=cff1*tot_ab(k,iband)+                 &
 
  839     &                            cff2*tot_b(k,iband)
 
  847                    datt(k,iband)=(tot_ab(k,iband)+                     &
 
  848     &                             tot_b(k,iband))/avgcos(k,iband)
 
  853                    avgcos_min=avgcos(k,iband)+                         &
 
  854     &                         (0.5_r8-avgcos(k,iband))*                &
 
  856     &                          (tot_ab(k,iband)+tot_s(k,iband)))
 
  863     &                      7.0_r8-datt(k,iband)*abs(z_r(i,j,k)))
 
  864                    slope_ac =min(0.0_r8,                               &
 
  865     &                            (avgcos_min-avgcos(k,iband))/fv1)
 
  866                    avgcos(k,iband)=avgcos(k,iband)+                    &
 
  867     &                              slope_ac*datt(k,iband)*hz(i,j,k)
 
  874     &                   0.005_r8*acos(avgcos(k,iband))*
rad2deg 
  875                    datt(k,iband)=cff1*tot_ab(k,iband)+                 &
 
  876     &                            cff2*tot_b(k,iband)
 
  878                    datt(k,iband)=(tot_ab(k,iband)+                     &
 
  879     &                             tot_b(k,iband))/avgcos(k,iband)
 
  885                      avgcos(k-1,iband)=avgcos(k,iband)
 
  890                    fv1=datt(k,iband)*hz(i,j,k)
 
  891                    fv2=datt_sum(iband)+0.5_r8*fv1
 
  892                    datt_sum(iband)=datt_sum(iband)+fv1
 
  893                    specir_d(k,iband)=specir(i,j,iband)*                &
 
  898                    specir_scal(i,k,iband)=specir_d(k,iband)*           &
 
  901                    e0_nz(i,k)=e0_nz(i,k)+specir_scal(i,k,iband)
 
  905                    ed_nz(i,k)=ed_nz(i,k)+specir_d(k,iband)
 
  906#ifdef DIAGNOSTICS_BIO 
  907                    diabio3d(i,j,iband,
idspir)=diabio3d(i,j,iband,      &
 
  910                    diabio4d(i,j,k,iband,
iddirr)=diabio4d(i,j,k,iband,  &
 
  913                    diabio4d(i,j,k,iband,
idsirr)=diabio4d(i,j,k,iband,  &
 
  915     &                                           specir_scal(i,k,iband)
 
  916                    diabio4d(i,j,k,iband,
idlatt)=diabio4d(i,j,k,iband,  &
 
  919                    diabio4d(i,j,k,iband,
idacos)=diabio4d(i,j,k,iband,  &
 
  944                IF ((bio(i,k,
idomc(
ilab)).gt.0.0_r8).and.               &
 
  945     &              (bio(i,k,
idomn(
ilab)).gt.0.0_r8).and.               &
 
  947                  nupdoc_ba(i,k,ibac)=gtbac(i,k,ibac)*                  &
 
  948     &                                bio(i,k,
ibacc(ibac))*             &
 
  953                  nupdon_ba(i,k,ibac)=nupdoc_ba(i,k,ibac)*              &
 
  956                  nupdop_ba(i,k,ibac)=nupdoc_ba(i,k,ibac)*              &
 
  960                  nupdoc_ba(i,k,ibac)=0.0_r8
 
  961                  nupdon_ba(i,k,ibac)=0.0_r8
 
  962                  nupdop_ba(i,k,ibac)=0.0_r8
 
  964                totdoc_d(i,k)=totdoc_d(i,k)+nupdoc_ba(i,k,ibac)
 
  965                totdon_d(i,k)=totdon_d(i,k)+nupdon_ba(i,k,ibac)
 
  966                totdop_d(i,k)=totdop_d(i,k)+nupdop_ba(i,k,ibac)
 
  970                nupnh4_ba(i,k,ibac)=gtbac(i,k,ibac)*                    &
 
  971     &                              bio(i,k,
ibacn(ibac))*               &
 
  974                totnh4_d(i,k)=totnh4_d(i,k)+nupnh4_ba(i,k,ibac)
 
  978                nuppo4_ba(i,k,ibac)=gtbac(i,k,ibac)*                    &
 
  979     &                              bio(i,k,
ibacp(ibac))*               &
 
  982                totpo4_d(i,k)=totpo4_d(i,k)+nuppo4_ba(i,k,ibac)
 
  986                nupfe_ba(i,k,ibac)=gtbac(i,k,ibac)*                     &
 
  987     &                             bio(i,k,
ibacf(ibac))*                &
 
  990                totfe_d(i,k)=totfe_d(i,k)+nupfe_ba(i,k,ibac)
 
 1009                  nup_max=gtalg(i,k,iphy)
 
 1010                  nupno3(i,k,iphy)=(bio(i,k,
ino3_)/                     &
 
 1013                  nupnh4(i,k,iphy)=bio(i,k,
inh4_)/                      &
 
 1018                  fv1=nupno3(i,k,iphy)+nupnh4(i,k,iphy)
 
 1019                  IF (fv1.gt.1.0_r8) 
THEN 
 1021                    nupno3(i,k,iphy)=nupno3(i,k,iphy)*fv1
 
 1022                    nupnh4(i,k,iphy)=nupnh4(i,k,iphy)*fv1
 
 1027                  fv1=nup_max*bio(i,k,
iphyn(iphy))
 
 1028                  nupno3(i,k,iphy)=nupno3(i,k,iphy)*fv1
 
 1029                  nupnh4(i,k,iphy)=nupnh4(i,k,iphy)*fv1
 
 1033                  IF (c2nalg(i,k,iphy).gt.
c2nnupdon(iphy,ng)) 
THEN 
 1034                    nupdon(i,k,iphy)=fv1*                               &
 
 1036     &                               (
hsdon(iphy,ng)+                   &
 
 1042                  totno3_d(i,k)=totno3_d(i,k)+nupno3(i,k,iphy)
 
 1043                  totnh4_d(i,k)=totnh4_d(i,k)+nupnh4(i,k,iphy)
 
 1044                  totdon_d(i,k)=totdon_d(i,k)+nupdon(i,k,iphy)
 
 1052                    nup_max=gtalg(i,k,iphy)
 
 1053                    nupsio(i,k,iphy)=bio(i,k,
isio_)/                    &
 
 1058                    IF (
iphys(iphy).gt.0) 
THEN 
 1059                      fv1=nup_max*bio(i,k,
iphys(iphy))
 
 1060                      nupsio(i,k,iphy)=nupsio(i,k,iphy)*fv1
 
 1062                      nupsio(i,k,iphy)=0.0_r8
 
 1067                    totsio_d(i,k)=totsio_d(i,k)+nupsio(i,k,iphy)
 
 1076                    nup_max=gtalg(i,k,iphy)
 
 1077                    nuppo4(i,k,iphy)=bio(i,k,
ipo4_)/                    &
 
 1082                    fv1=nup_max*bio(i,k,
iphyp(iphy))
 
 1083                    nuppo4(i,k,iphy)=nuppo4(i,k,iphy)*fv1
 
 1087                    IF (c2palg(i,k,iphy).gt.
c2palkphos(iphy,ng)) 
THEN 
 1088                      nupdop(i,k,iphy)=fv1*                             &
 
 1090     &                                 (
hsdop(iphy,ng)+                 &
 
 1096                    totpo4_d(i,k)=totpo4_d(i,k)+nuppo4(i,k,iphy)
 
 1097                    totdop_d(i,k)=totdop_d(i,k)+nupdop(i,k,iphy)
 
 1106                    nup_max=gtalg(i,k,iphy)
 
 1107                    nupfe(i,k,iphy)=bio(i,k,
ifeo_)/                     &
 
 1112                    fv1=nup_max*bio(i,k,
iphyf(iphy))
 
 1113                    nupfe(i,k,iphy)=nupfe(i,k,iphy)*fv1
 
 1117                    totfe_d(i,k)=totfe_d(i,k)+nupfe(i,k,iphy)
 
 1131              nh4tono3(i,k)=0.0_r8
 
 1134              ntofebac(i,k)=0.0_r8
 
 1135              IF (k.lt.keuphotic(i)) 
THEN 
 1136                nh4tono3(i,k)=
rtnit(ng)*                                &
 
 1143                nitrbac(i,k)=nh4tono3(i,k)/7.0_r8
 
 1144                ntonbac(i,k)=nitrbac(i,k)*
n2cbac(ng)
 
 1145                ntopbac(i,k)=nitrbac(i,k)*
p2cbac(ng)
 
 1146                ntofebac(i,k)=nitrbac(i,k)*
fe2cbac(ng)
 
 1147                totnh4_d(i,k)=totnh4_d(i,k)+nh4tono3(i,k)+ntonbac(i,k)
 
 1148                totpo4_d(i,k)=totpo4_d(i,k)+ntopbac(i,k)
 
 1149                totfe_d(i,k)=totfe_d(i,k)+ntofebac(i,k)
 
 1162              fv2=totno3_d(i,k)*dtbio
 
 1163              IF (fv2.gt.bio(i,k,
ino3_)) 
THEN 
 1166                  nupno3(i,k,iphy)=nupno3(i,k,iphy)*fv1
 
 1170              fv2=totnh4_d(i,k)*dtbio
 
 1171              IF (fv2.gt.bio(i,k,
inh4_)) 
THEN 
 1174                  nupnh4(i,k,iphy)=nupnh4(i,k,iphy)*fv1
 
 1177                  nupnh4_ba(i,k,ibac)=nupnh4_ba(i,k,ibac)*fv1
 
 1179                nh4tono3(i,k)=nh4tono3(i,k)*fv1
 
 1180                ntonbac(i,k)=ntonbac(i,k)*fv1
 
 1183              fv2=totsio_d(i,k)*dtbio
 
 1184              IF (fv2.gt.bio(i,k,
isio_)) 
THEN 
 1187                  nupsio(i,k,iphy)=nupsio(i,k,iphy)*fv1
 
 1191              fv2=totpo4_d(i,k)*dtbio
 
 1192              IF (fv2.gt.bio(i,k,
ipo4_)) 
THEN 
 1195                  nuppo4(i,k,iphy)=nuppo4(i,k,iphy)*fv1
 
 1198                  nuppo4_ba(i,k,ibac)=nuppo4_ba(i,k,ibac)*fv1
 
 1200                ntopbac(i,k)=ntopbac(i,k)*fv1
 
 1203              fv2=totfe_d(i,k)*dtbio
 
 1204              IF (fv2.gt.bio(i,k,
ifeo_)) 
THEN 
 1207                  nupfe(i,k,iphy)=nupfe(i,k,iphy)*fv1
 
 1210                  nupfe_ba(i,k,ibac)=nupfe_ba(i,k,ibac)*fv1
 
 1212                ntofebac(i,k)=ntofebac(i,k)*fv1
 
 1218              fv2=totdoc_d(i,k)*dtbio
 
 1221                totdoc_d(i,k)=totdoc_d(i,k)*fv1
 
 1223                  nupdoc_ba(i,k,ibac)=nupdoc_ba(i,k,ibac)*fv1
 
 1224                  totdon_d(i,k)=totdon_d(i,k)-nupdon_ba(i,k,ibac)
 
 1225                  nupdon_ba(i,k,ibac)=nupdon_ba(i,k,ibac)*fv1
 
 1226                  totdon_d(i,k)=totdon_d(i,k)+nupdon_ba(i,k,ibac)
 
 1227                  totdop_d(i,k)=totdop_d(i,k)-nupdop_ba(i,k,ibac)
 
 1228                  nupdop_ba(i,k,ibac)=nupdop_ba(i,k,ibac)*fv1
 
 1229                  totdop_d(i,k)=totdop_d(i,k)+nupdop_ba(i,k,ibac)
 
 1235              fv2=totdon_d(i,k)*dtbio
 
 1238                totdon_d(i,k)=totdon_d(i,k)*fv1
 
 1239                totdoc_d(i,k)=totdoc_d(i,k)*fv1
 
 1241                  nupdon(i,k,iphy)=nupdon(i,k,iphy)*fv1
 
 1244                  nupdon_ba(i,k,ibac)=nupdon_ba(i,k,ibac)*fv1
 
 1245                  nupdoc_ba(i,k,ibac)=nupdoc_ba(i,k,ibac)*fv1
 
 1246                  totdop_d(i,k)=totdop_d(i,k)-nupdop_ba(i,k,ibac)
 
 1247                  nupdop_ba(i,k,ibac)=nupdop_ba(i,k,ibac)*fv1
 
 1248                  totdop_d(i,k)=totdop_d(i,k)+nupdop_ba(i,k,ibac)
 
 1254              fv2=totdop_d(i,k)*dtbio
 
 1257                totdop_d(i,k)=totdop_d(i,k)*fv1
 
 1258                totdoc_d(i,k)=totdoc_d(i,k)*fv1
 
 1260                  nupdop(i,k,iphy)=nupdop(i,k,iphy)*fv1
 
 1263                  nupdop_ba(i,k,ibac)=nupdop_ba(i,k,ibac)*fv1
 
 1264                  totdon_d(i,k)=totdon_d(i,k)-nupdon_ba(i,k,ibac)
 
 1265                  nupdon_ba(i,k,ibac)=nupdon_ba(i,k,ibac)*fv1
 
 1266                  totdon_d(i,k)=totdon_d(i,k)+nupdon_ba(i,k,ibac)
 
 1267                  nupdoc_ba(i,k,ibac)=nupdoc_ba(i,k,ibac)*fv1
 
 1278                bio_new(i,k,
iphyn(iphy))=bio_new(i,k,
iphyn(iphy))+      &
 
 1279     &                                   nupno3(i,k,iphy)+              &
 
 1280     &                                   nupnh4(i,k,iphy)+              &
 
 1282                bio_new(i,k,
iphyp(iphy))=bio_new(i,k,
iphyp(iphy))+      &
 
 1283     &                                   nuppo4(i,k,iphy)+              &
 
 1285                bio_new(i,k,
iphyf(iphy))=bio_new(i,k,
iphyf(iphy))+      &
 
 1287                IF (
iphys(iphy).gt.0) 
THEN 
 1288                  bio_new(i,k,
iphys(iphy))=bio_new(i,k,
iphys(iphy))+    &
 
 1328                bio_new(i,k,
ibacc(ibac))=bio_new(i,k,
ibacc(ibac))+      &
 
 1330                bio_new(i,k,
ibacn(ibac))=bio_new(i,k,
ibacn(ibac))+      &
 
 1332                bio_new(i,k,
ibacp(ibac))=bio_new(i,k,
ibacp(ibac))+      &
 
 1334                bio_new(i,k,
ibacf(ibac))=bio_new(i,k,
ibacf(ibac))+      &
 
 1347     &                           (nh4tono3(i,k)+ntonbac(i,k))
 
 1360            DO k=n(ng),keuphotic(i),-1
 
 1362                IF (bio(i,k,
iphyc(iphy)).gt.0.0_r8) 
THEN 
 1368                    aphyn_wa=aphyn_wa+(aphyn_al(i,k,iphy,iband)*        &
 
 1369     &                                 specir_scal(i,k,iband))
 
 1376                  aphyn_wa=aphyn_wa/e0_nz(i,k)
 
 1381                  alfa(k,iphy)=(aphyn_wa/bio(i,k,
iphyc(iphy)))*         &
 
 1382     &                          
qu_yld(iphy,ng)*0.001_r8
 
 1386                  fv1=max(0.0_r8,e0_nz(i,k)-
e0_comp(iphy,ng))
 
 1388                  IF (fv2.gt.0.0_r8) 
THEN 
 1389                    gt_ll(k,iphy)=gtalg(i,k,iphy)*                      &
 
 1390     &                            tanh(alfa(k,iphy)*fv1/                &
 
 1391     &                            gtalg(i,k,iphy))*                     &
 
 1394                    gt_ll(k,iphy)=gtalg(i,k,iphy)*                      &
 
 1395     &                            tanh(alfa(k,iphy)*fv1/                &
 
 1407                  IF (bio(i,k,
iphyn(iphy)).gt.0.0_r8) 
THEN 
 1408                    fv1=bio(i,k,
iphyc(iphy))/                           &
 
 1409     &                  (bio(i,k,
iphyn(iphy))+bio_new(i,k,
iphyn(iphy)))
 
 1410                    gt_nl(k,iphy)=mu_bar_n(i,k,iphy)*                   &
 
 1412                    gt_nl(k,iphy)=max(0.0_r8,                           &
 
 1413     &                                min(gt_nl(k,iphy),                &
 
 1420                  IF (
iphys(iphy).gt.0) 
THEN 
 1422     &                  (bio(i,k,
iphys(iphy)).gt.0.0_r8)) 
THEN 
 1423                      fv1=bio(i,k,
iphyc(iphy))/                         &
 
 1424     &                    (bio(i,k,
iphys(iphy))+                        &
 
 1425     &                     bio_new(i,k,
iphys(iphy)))
 
 1426                      gt_sl(k,iphy)=mu_bar_s(i,k,iphy)*                 &
 
 1428                      gt_sl(k,iphy)=max(0.0_r8,                         &
 
 1429     &                                  min(gt_sl(k,iphy),              &
 
 1441     &                (bio(i,k,
iphyp(iphy)).gt.0.0_r8)) 
THEN 
 1442                    fv1=bio(i,k,
iphyc(iphy))/                           &
 
 1443     &                  (bio(i,k,
iphyp(iphy))+bio_new(i,k,
iphyp(iphy)))
 
 1444                    gt_pl(k,iphy)=mu_bar_p(i,k,iphy)*                   &
 
 1446                    gt_pl(k,iphy)=max(0.0_r8,                           &
 
 1447     &                                min(gt_pl(k,iphy),                &
 
 1456     &                (bio(i,k,
iphyf(iphy)).gt.0.0_r8)) 
THEN 
 1457                    fv1=bio(i,k,
iphyc(iphy))/                           &
 
 1458     &                  (bio(i,k,
iphyf(iphy))+bio_new(i,k,
iphyf(iphy)))
 
 1459                    gt_fl(k,iphy)=mu_bar_f(i,k,iphy)*                   &
 
 1461                    gt_fl(k,iphy)=max(0.0_r8,                           &
 
 1462     &                                min(gt_fl(k,iphy),                &
 
 1470                  gtalg_r(i,k,iphy)=min(gt_ll(k,iphy),gt_nl(k,iphy),    &
 
 1471     &                                  gt_sl(k,iphy),gt_pl(k,iphy),    &
 
 1473                  IF (gtalg_r(i,k,iphy).ge.
larger) 
THEN 
 1474                    gtalg_r(i,k,iphy)=0.0_r8
 
 1479                  fv1=bio(i,k,
iphyc(iphy))*gtalg_r(i,k,iphy)
 
 1480                  bio_new(i,k,
iphyc(iphy))=bio_new(i,k,
iphyc(iphy))+    &
 
 1488                    IF (
ipigs(iphy,ipig).gt.0) 
THEN 
 1489                      itrc=
ipigs(iphy,ipig)
 
 1490                      IF (bio(i,k,
iphyc(iphy)).gt.0.0_r8) 
THEN 
 1491                        fv1=bio(i,k,itrc)*gtalg_r(i,k,iphy)
 
 1492                        bio_new(i,k,itrc)=bio_new(i,k,itrc)+fv1
 
 1520                fv1=nupdoc_ba(i,k,ibac)*
exbac_c(ng)*                    &
 
 1522                fv2=nupdoc_ba(i,k,ibac)*
exbac_c(ng)*                    &
 
 1524                fv3=nupdon_ba(i,k,ibac)*
exbac_n(ng)
 
 1542                nupdoc_ba(i,k,ibac)=nupdoc_ba(i,k,ibac)-                &
 
 1544                nupdon_ba(i,k,ibac)=nupdon_ba(i,k,ibac)-                &
 
 1551                bac_g(1)=nupdoc_ba(i,k,ibac)*
bac_ceff(ng)
 
 1552                bac_g(2)=(nupdon_ba(i,k,ibac)+                          &
 
 1553     &                    nupnh4_ba(i,k,ibac))*                         &
 
 1555                bac_g(3)=(nupdop_ba(i,k,ibac)+                          &
 
 1556     &                    nuppo4_ba(i,k,ibac))*                         &
 
 1558                bac_g(4)=nupfe_ba(i,k,ibac)*
c2febac(ng)
 
 1562                IF ((bac_g(1).le.bac_g(2)).and.                         &
 
 1563     &              (bac_g(1).le.bac_g(3)).and.                         &
 
 1564     &              (bac_g(1).le.bac_g(4))) 
THEN 
 1569                  bio_new(i,k,
ibacn(ibac))=bio_new(i,k,
ibacn(ibac))+    &
 
 1571                  bio_new(i,k,
ibacp(ibac))=bio_new(i,k,
ibacp(ibac))+    &
 
 1573                  bio_new(i,k,
ibacf(ibac))=bio_new(i,k,
ibacf(ibac))+    &
 
 1580                  nupnh4_ba(i,k,ibac)=fv1-nupdon_ba(i,k,ibac)
 
 1581                  nuppo4_ba(i,k,ibac)=fv2-nupdop_ba(i,k,ibac)
 
 1586                  relfe=nupfe_ba(i,k,ibac)-fv3
 
 1587                  nupfe_ba(i,k,ibac)=fv3
 
 1592                ELSE IF ((bac_g(2).le.bac_g(3)).and.                    &
 
 1593     &                   (bac_g(2).le.bac_g(4))) 
THEN 
 1597                  bio_new(i,k,
ibacn(ibac))=bio_new(i,k,
ibacn(ibac))+    &
 
 1598     &                                     (nupdon_ba(i,k,ibac)+        &
 
 1599     &                                      nupnh4_ba(i,k,ibac))
 
 1600                  bio_new(i,k,
ibacp(ibac))=bio_new(i,k,
ibacp(ibac))+    &
 
 1602                  bio_new(i,k,
ibacf(ibac))=bio_new(i,k,
ibacf(ibac))+    &
 
 1609                  nupdoc_ba(i,k,ibac)=nupdoc_ba(i,k,ibac)-fv1
 
 1618                  fv5=fv2-(nupdop_ba(i,k,ibac)+                         &
 
 1619                           nuppo4_ba(i,k,ibac)-fv4)
 
 1623                  IF (fv5.lt.0.0_r8) 
THEN 
 1625                    nuppo4_ba(i,k,ibac)=nuppo4_ba(i,k,ibac)+fv5
 
 1629                  nupdop_ba(i,k,ibac)=nupdop_ba(i,k,ibac)-reldop1
 
 1633                  relfe=nupfe_ba(i,k,ibac)-fv3
 
 1634                  nupfe_ba(i,k,ibac)=fv3
 
 1639                ELSE IF (bac_g(3).le.bac_g(4)) 
THEN 
 1643                  bio_new(i,k,
ibacn(ibac))=bio_new(i,k,
ibacn(ibac))+    &
 
 1645                  bio_new(i,k,
ibacp(ibac))=bio_new(i,k,
ibacp(ibac))+    &
 
 1646     &                                     (nupdop_ba(i,k,ibac)+        &
 
 1647     &                                      nuppo4_ba(i,k,ibac))
 
 1648                  bio_new(i,k,
ibacf(ibac))=bio_new(i,k,
ibacf(ibac))+    &
 
 1655                  nupdoc_ba(i,k,ibac)=nupdoc_ba(i,k,ibac)-fv1
 
 1664                  fv5=fv2-(nupdon_ba(i,k,ibac)+                         &
 
 1665     &                     nupnh4_ba(i,k,ibac)-fv4)
 
 1669                  IF (fv5.lt.0.0_r8) 
THEN 
 1671                    nupnh4_ba(i,k,ibac)=nupnh4_ba(i,k,ibac)+fv5
 
 1675                  nupdon_ba(i,k,ibac)=nupdon_ba(i,k,ibac)-reldon1
 
 1679                  relfe=nupfe_ba(i,k,ibac)-fv3
 
 1680                  nupfe_ba(i,k,ibac)=fv3
 
 1689                  bio_new(i,k,
ibacn(ibac))=bio_new(i,k,
ibacn(ibac))+    &
 
 1691                  bio_new(i,k,
ibacp(ibac))=bio_new(i,k,
ibacp(ibac))+    &
 
 1693                  bio_new(i,k,
ibacf(ibac))=bio_new(i,k,
ibacf(ibac))+    &
 
 1694     &                                     nupfe_ba(i,k,ibac)
 
 1700                  nupdoc_ba(i,k,ibac)=nupdoc_ba(i,k,ibac)-fv1
 
 1709                  fv5=fv2-(nupdon_ba(i,k,ibac)+                         &
 
 1710     &                     nupnh4_ba(i,k,ibac)-fv4)
 
 1714                  IF (fv5.lt.0.0_r8) 
THEN 
 1716                    nupnh4_ba(i,k,ibac)=nupnh4_ba(i,k,ibac)+fv5
 
 1720                  nupdon_ba(i,k,ibac)=nupdon_ba(i,k,ibac)-reldon1
 
 1728                  fv5=fv2-(nupdop_ba(i,k,ibac)+                         &
 
 1729     &                     nuppo4_ba(i,k,ibac)-fv4)
 
 1733                  IF (fv5.lt.0.0_r8) 
THEN 
 1735                    nuppo4_ba(i,k,ibac)=nuppo4_ba(i,k,ibac)+fv5
 
 1739                  nupdop_ba(i,k,ibac)=nupdop_ba(i,k,ibac)-reldop1
 
 1744                bio_new(i,k,
ibacc(ibac))=bio_new(i,k,
ibacc(ibac))+      &
 
 1746                fv1=nupdoc_ba(i,k,ibac)-het_bac
 
 1755     &                                   (totdoc_d(i,k)-reldoc1)
 
 1765     &                                   nupdon_ba(i,k,ibac)
 
 1767     &                                   nupdop_ba(i,k,ibac)
 
 1769     &                             nupnh4_ba(i,k,ibac)
 
 1771     &                             nuppo4_ba(i,k,ibac)
 
 1773     &                             nupfe_ba(i,k,ibac)
 
 1788                IF ((c2nalg(i,k,iphy).ge.                               &
 
 1790     &              (c2palg(i,k,iphy).ge.                               &
 
 1794                  bio_new(i,k,
iphyc(iphy))=bio_new(i,k,
iphyc(iphy))-    &
 
 1801                ELSE IF ((c2nalg(i,k,iphy).ge.                          &
 
 1803     &                   (c2palg(i,k,iphy).ge.                          &
 
 1805     &                   (c2salg(i,k,iphy).ge.                          &
 
 1808                  bio_new(i,k,
iphyc(iphy))=bio_new(i,k,
iphyc(iphy))-    &
 
 1819                IF (bio(i,k,
iphyc(iphy)).gt.refuge(i,k,iphy)) 
THEN 
 1823                  fv1=graz_act(i,k,iphy)*bio(i,k,
iphyc(iphy))
 
 1824                  bio_new(i,k,
iphyc(iphy))=bio_new(i,k,
iphyc(iphy))-    &
 
 1841                  fv2=graz_act(i,k,iphy)*bio(i,k,
iphyn(iphy))
 
 1842                  bio_new(i,k,
iphyn(iphy))=bio_new(i,k,
iphyn(iphy))-    &
 
 1855                  IF (
iphys(iphy).gt.0) 
THEN 
 1856                    fv2=graz_act(i,k,iphy)*bio(i,k,
iphys(iphy))
 
 1857                    bio_new(i,k,
iphys(iphy))=bio_new(i,k,
iphys(iphy))-  &
 
 1866     &                                       (1.0_r8-
fecdoc(iphy,ng))*  &
 
 1872                  fv2=graz_act(i,k,iphy)*bio(i,k,
iphyp(iphy))
 
 1873                  bio_new(i,k,
iphyp(iphy))=bio_new(i,k,
iphyp(iphy))-    &
 
 1886                  fv2=graz_act(i,k,iphy)*bio(i,k,
iphyf(iphy))
 
 1887                  bio_new(i,k,
iphyf(iphy))=bio_new(i,k,
iphyf(iphy))-    &
 
 1905              IF (
ipigs(iphy,ipig).gt.0) 
THEN 
 1906                itrc=
ipigs(iphy,ipig)
 
 1909                    IF (bio(i,k,
iphyc(iphy)).gt.refuge(i,k,iphy)) 
THEN 
 1910                      fv1=graz_act(i,k,iphy)*bio(i,k,itrc)
 
 1911                      bio_new(i,k,itrc)=bio_new(i,k,itrc) - fv1
 
 1943                bio_new(i,k,
ibacc(ibac))=bio_new(i,k,
ibacc(ibac))-      &
 
 1945     &                                   bio_new(i,k,
ibacc(ibac))
 
 1948     &                                   bio_new(i,k,
ibacc(ibac))*      &
 
 1953     &                                   bio_new(i,k,
ibacc(ibac))*      &
 
 1958     &                                   bio_new(i,k,
ibacc(ibac))*      &
 
 1962     &                             bio_new(i,k,
ibacc(ibac))*            &
 
 1967                bio_new(i,k,
ibacn(ibac))=bio_new(i,k,
ibacn(ibac))-      &
 
 1969     &                                   bio_new(i,k,
ibacn(ibac))
 
 1972     &                                   bio_new(i,k,
ibacn(ibac))*      &
 
 1976     &                                   bio_new(i,k,
ibacn(ibac))*      &
 
 1980     &                             bio_new(i,k,
ibacn(ibac))*            &
 
 1985                bio_new(i,k,
ibacp(ibac))=bio_new(i,k,
ibacp(ibac))-      &
 
 1987     &                                   bio_new(i,k,
ibacp(ibac))
 
 1990     &                                   bio_new(i,k,
ibacp(ibac))*      &
 
 1994     &                                   bio_new(i,k,
ibacp(ibac))*      &
 
 1998     &                             bio_new(i,k,
ibacp(ibac))*            &
 
 2003                bio_new(i,k,
ibacf(ibac))=bio_new(i,k,
ibacf(ibac))-      &
 
 2005     &                                   bio_new(i,k,
ibacf(ibac))
 
 2008     &                                   bio_new(i,k,
ibacf(ibac))*      &
 
 2012     &                             bio_new(i,k,
ibacf(ibac))*            &
 
 2028                fv3=regen_c(i,k,ifec)*bio(i,k,
ifecc(ifec))
 
 2029                bio_new(i,k,
ifecc(ifec))=bio_new(i,k,
ifecc(ifec))-      &
 
 2036                fv2=regen_n(i,k,ifec)*bio(i,k,
ifecn(ifec))
 
 2037                bio_new(i,k,
ifecn(ifec))=bio_new(i,k,
ifecn(ifec))-      &
 
 2044                fv2=regen_s(i,k,ifec)*bio(i,k,
ifecs(ifec))
 
 2045                bio_new(i,k,
ifecs(ifec))=bio_new(i,k,
ifecs(ifec))-      &
 
 2052                fv2=regen_p(i,k,ifec)*bio(i,k,
ifecp(ifec))
 
 2053                bio_new(i,k,
ifecp(ifec))=bio_new(i,k,
ifecp(ifec))-      &
 
 2060                fv2=regen_f(i,k,ifec)*bio(i,k,
ifecf(ifec))
 
 2061                bio_new(i,k,
ifecf(ifec))=bio_new(i,k,
ifecf(ifec))-      &
 
 2079              IF (ed_nz(i,n(ng)).ge.0.01) 
THEN 
 2081                fv1=
rtuvr_dic(ng)*ed_nz(i,n(ng))/1500.0_r8
 
 2082                fv2=
rtuvr_doc(ng)*ed_nz(i,n(ng))/1500.0_r8
 
 2091                photo_decay=0.5_r8*hz(i,j,n(ng))*                       &
 
 2093                fv3=exp(-photo_decay)
 
 2094                photo_decay=2.0_r8*photo_decay
 
 2098                DO k=n(ng),keuphotic(i),-1
 
 2099                  IF (fv3.gt.0.01_r8) 
THEN 
 2101                    IF (fv6.gt.0.0_r8) 
THEN 
 2103                      photo_dic=fv3*fv1*fv6
 
 2104                      photo_doc=fv3*fv2*fv6
 
 2105                      total_photo=photo_dic+photo_doc
 
 2109                      fv4=(1.0_r8-fv7)*total_photo
 
 2127     &                 0.5_r8*hz(i,j,k)*(0.2_r8+(fv4+fv5)*
adoc300(
ilab))
 
 2137                    photo_decay=photo_decay+2.0_r8*fv7
 
 2149            IF (keuphotic(i).le.n(ng)) 
THEN 
 2157                DO k=n(ng),keuphotic(i),-1
 
 2159                    c2chl_w(k,iphy)=min((
b_c2cl(iphy,ng)+               &
 
 2160     &                                   
mxc2cl(iphy,ng)*e0_nz(i,k)),   &
 
 2162                  ELSE IF (c2nalg(i,k,iphy).gt.                         &
 
 2164                    c2chl_w(k,iphy)=
b_c2cn(iphy,ng)+                    &
 
 2166     &                              (c2nalg(i,k,iphy)-                  &
 
 2169                    c2chl_w(k,iphy)=min((
b_c2cl(iphy,ng)+               &
 
 2170     &                                   
mxc2cl(iphy,ng)*e0_nz(i,k)),   &
 
 2177                DO k=n(ng),keuphotic(i),-1
 
 2178                  pigs_w(k,iphy,
ichl)=1.0_r8/c2chl_w(k,iphy)
 
 2183                IF (
ipigs(iphy,2).gt.0) 
THEN 
 2184                  DO k=n(ng),keuphotic(i),-1
 
 2185                    pigs_w(k,iphy,2)=
b_chlb(iphy,ng)+                   &
 
 2187     &                               (c2chl_w(k,iphy)-                  &
 
 2189                    pigs_w(k,iphy,2)=pigs_w(k,iphy,2)*                  &
 
 2190     &                               pigs_w(k,iphy,
ichl)
 
 2196                IF (
ipigs(iphy,3).gt.0) 
THEN 
 2197                  DO k=n(ng),keuphotic(i),-1
 
 2198                    pigs_w(k,iphy,3)=
b_chlc(iphy,ng)+                   &
 
 2200     &                               (c2chl_w(k,iphy)-                  &
 
 2202                    pigs_w(k,iphy,3)=pigs_w(k,iphy,3)*                  &
 
 2203     &                               pigs_w(k,iphy,
ichl)
 
 2209                IF (
ipigs(iphy,4).gt.0) 
THEN 
 2210                  DO k=n(ng),keuphotic(i),-1
 
 2211                    pigs_w(k,iphy,4)=
b_psc(iphy,ng)+                    &
 
 2213     &                               (c2chl_w(k,iphy)-                  &
 
 2215                    pigs_w(k,iphy,4)=pigs_w(k,iphy,4)*                  &
 
 2216     &                               pigs_w(k,iphy,
ichl)
 
 2222                IF (
ipigs(iphy,5).gt.0) 
THEN 
 2223                  DO k=n(ng),keuphotic(i),-1
 
 2224                    pigs_w(k,iphy,5)=
b_ppc(iphy,ng)+                    &
 
 2226     &                               (c2chl_w(k,iphy)-                  &
 
 2228                    pigs_w(k,iphy,5)=pigs_w(k,iphy,5)*                  &
 
 2229     &                               pigs_w(k,iphy,
ichl)
 
 2235                IF (
ipigs(iphy,6).gt.0) 
THEN 
 2236                  DO k=n(ng),keuphotic(i),-1
 
 2237                    pigs_w(k,iphy,6)=
b_lpub(iphy,ng)+                   &
 
 2239     &                               (c2chl_w(k,iphy)-                  &
 
 2241                    pigs_w(k,iphy,6)=pigs_w(k,iphy,6)*                  &
 
 2242     &                               pigs_w(k,iphy,
ichl)
 
 2248                IF (
ipigs(iphy,7).gt.0) 
THEN 
 2249                  DO k=n(ng),keuphotic(i),-1
 
 2250                    pigs_w(k,iphy,7)=
b_hpub(iphy,ng)+                   &
 
 2252     &                               (c2chl_w(k,iphy)-                  &
 
 2254                    pigs_w(k,iphy,7)=pigs_w(k,iphy,7)*                  &
 
 2255     &                               pigs_w(k,iphy,
ichl)
 
 2265                  IF (
ipigs(iphy,ipig).gt.0) 
THEN 
 2266                    itrc=
ipigs(iphy,ipig)
 
 2267                    DO k=n(ng),keuphotic(i),-1
 
 2268                      IF ((bio(i,k,
iphyc(iphy)).gt.0.0_r8).and.         &
 
 2269     &                    (bio(i,k,itrc).gt.0.0_r8)) 
THEN 
 2270                        fv1=bio(i,k,
iphyc(iphy))*12.0_r8
 
 2271                        fv2=gtalg_r(i,k,iphy)
 
 2273     &                      (fv2/pigs_w(k,iphy,ipig)+                   &
 
 2274     &                       fv1*(1.0_r8-fv2)/                          &
 
 2276                        bio_new(i,k,itrc)=bio_new(i,k,itrc)+            &
 
 2277     &                                    (fv3-bio(i,k,itrc))
 
 2294          sink_loop: 
DO isink=1,nsink
 
 2304                qc(i,k)=bio(i,k,itrc)
 
 2310                fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
 
 2315                dltr=hz(i,j,k)*fc(i,k)
 
 2316                dltl=hz(i,j,k)*fc(i,k-1)
 
 2317                cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
 
 2324                IF ((dltr*dltl).le.0.0_r8) 
THEN 
 2327                ELSE IF (abs(dltr).gt.abs(cffl)) 
THEN 
 2329                ELSE IF (abs(dltl).gt.abs(cffr)) 
THEN 
 2342                cff=(dltr-dltl)*hz_inv3(i,k)
 
 2343                dltr=dltr-cff*hz(i,j,k+1)
 
 2344                dltl=dltl+cff*hz(i,j,k-1)
 
 2345                br(i,k)=qc(i,k)+dltr
 
 2346                bl(i,k)=qc(i,k)-dltl
 
 2347                wr(i,k)=(2.0_r8*dltr-dltl)**2
 
 2348                wl(i,k)=(dltr-2.0_r8*dltl)**2
 
 2354                dltl=max(cff,wl(i,k  ))
 
 2355                dltr=max(cff,wr(i,k+1))
 
 2356                br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
 
 2362#if defined LINEAR_CONTINUATION 
 2363              bl(i,n(ng))=br(i,n(ng)-1)
 
 2364              br(i,n(ng))=2.0_r8*qc(i,n(ng))-bl(i,n(ng))
 
 2365#elif defined NEUMANN 
 2366              bl(i,n(ng))=br(i,n(ng)-1)
 
 2367              br(i,n(ng))=1.5_r8*qc(i,n(ng))-0.5_r8*bl(i,n(ng))
 
 2369              br(i,n(ng))=qc(i,n(ng))       
 
 2370              bl(i,n(ng))=qc(i,n(ng))       
 
 2371              br(i,n(ng)-1)=qc(i,n(ng))
 
 2373#if defined LINEAR_CONTINUATION 
 2375              bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
 
 2376#elif defined NEUMANN 
 2378              bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
 
 2392                dltr=br(i,k)-qc(i,k)
 
 2393                dltl=qc(i,k)-bl(i,k)
 
 2396                IF ((dltr*dltl).lt.0.0_r8) 
THEN 
 2399                ELSE IF (abs(dltr).gt.abs(cffl)) 
THEN 
 2401                ELSE IF (abs(dltl).gt.abs(cffr)) 
THEN 
 2404                br(i,k)=qc(i,k)+dltr
 
 2405                bl(i,k)=qc(i,k)-dltl
 
 2423            cff=dtbio*abs(wbio(isink))
 
 2427                wl(i,k)=z_w(i,j,k-1)+cff
 
 2428                wr(i,k)=hz(i,j,k)*qc(i,k)
 
 2435                  IF (wl(i,k).gt.z_w(i,j,ks)) 
THEN 
 2437                    fc(i,k-1)=fc(i,k-1)+wr(i,ks)
 
 2448                cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
 
 2449                fc(i,k-1)=fc(i,k-1)+                                    &
 
 2452     &                     cu*(0.5_r8*(br(i,ks)-bl(i,ks))-              &
 
 2454     &                         (br(i,ks)+bl(i,ks)-                      &
 
 2455     &                          2.0_r8*qc(i,ks))))
 
 2460                bio(i,k,itrc)=qc(i,k)+(fc(i,k)-fc(i,k-1))*hz_inv(i,k)
 
 2474              IF (itrc.eq.
ifecn(ifec)) 
THEN 
 2476                  cff1=fc(i,0)*hz_inv(i,1)
 
 2479              ELSE IF (itrc.eq.
ifecc(ifec)) 
THEN 
 2481                  cff1=fc(i,0)*hz_inv(i,1)
 
 2484              ELSE IF (itrc.eq.
ifecp(ifec)) 
THEN 
 2486                  cff1=fc(i,0)*hz_inv(i,1)
 
 2489              ELSE IF (itrc.eq.
ifecs(ifec)) 
THEN 
 2491                  cff1=fc(i,0)*hz_inv(i,1)
 
 2494              ELSE IF (itrc.eq.
ifecf(ifec)) 
THEN 
 2496                  cff1=fc(i,0)*hz_inv(i,1)
 
 2502              IF (itrc.eq.
iphyn(iphy)) 
THEN 
 2504                  cff1=fc(i,0)*hz_inv(i,1)
 
 2507              ELSE IF (itrc.eq.
iphyc(iphy)) 
THEN 
 2509                  cff1=fc(i,0)*hz_inv(i,1)
 
 2512              ELSE IF (itrc.eq.
iphyp(iphy)) 
THEN 
 2514                  cff1=fc(i,0)*hz_inv(i,1)
 
 2517              ELSE IF (itrc.eq.
iphys(iphy)) 
THEN 
 2519                  cff1=fc(i,0)*hz_inv(i,1)
 
 2522              ELSE IF (itrc.eq.
iphyf(iphy)) 
THEN 
 2524                  cff1=fc(i,0)*hz_inv(i,1)
 
 2540                  bio(i,k,itrc)=bio(i,k,itrc)+dtbio*bio_new(i,k,itrc)
 
 2566               cff=bio(i,k,itrc)-bio_old(i,k,itrc)
 
 2567               t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)+cff*hz(i,j,k)