      SUBROUTINE BIOMASSE_old(TRBIO,FTRBIO,NTRBIO,NFLBIO,H,NZ,NZI,NPT,T,KBIO,
     *                XZM,XZE,XPAR,
     >                XPARZE,XPARZM,FLXBIO,XLE,XLM,XLT,XLNUT,flor)
CCC---------------------------------------------------------------------
CCC
CCC                       ROUTINE BIOMASSE
CCC                     *******************
CCC
CCC  PURPOSE :
CCC  ---------
CCC     COMPUTE THE FORCING TREND DUE TO BIOGEOCHEMICAL PROCESSES 
CCC     AND ADD IT TO THE GENERAL FORCING TREND OF PASSIVE TRACERS EQUATIONS.
CCC     WARNING : FORCING ARE FLUXES, ie FOR THE EQUATION OF C*H 
CCC     WHERE C IS THE TRACER AND H THE VERTICAL LAYER SCALE.
CCC
CCC
CC   METHOD :
CC   -------
CC      EACH NOW BIOLOGICAL FLUX IS CALCULATED  IN FUNCTION OF NOW 
CC      CONCENTRATIONS OF TRACERS.
CC      DEPENDING ON THE TRACER, THESE FLUXES ARE SOURCES OR SINKS.
CC      THE TOTAL OF THE SOURCES AND SINKS FOR EACH TRACER
CC      IS ADDED TO THE GENERAL TREND.
CC
CC		TRA = TRA + ZF...TRA - ZFTRA...
CC                             |         |
CC                             |         |
CC                          SOURCE      SINK
CC
CC
CC
CC
CC   INPUT :
CC   -----
CC      ARGUMENT
CC              TRBIO           : biological tracers array (in mmoleN/m3)
CC              FTRBIO          : biological tracer forcing (in mmoleN/m2/s)
CC              NTRBIO          : number of biological tracers
CC              NFLBIO          : number of biological fluxes
CC              KBIO            : number of biologicaly active vertical layers
CC              H               : vertical scale factors at T points (in m)
CC              NZ              : number of maximum vertical layers
CC              NZI             : number of vertical layers
CC              NPT             : number of horizontal grid points
CC              T               : temperature (in Celsius)
CC              XZM             : MIXED LAYER DEPTH (in m)
CC              XZE             : EUPHOTIC LAYER DEPTH (in m)
CC              XPAR            : PAR AT EACH T LEVEL
CC              XPARZE          : MEAN PAR IN THE EUPHOTIC LAYER
CC              XPARZM          : MEAN PAR IN THE MIXED LAYER
CC
CC      COMMON
CC            /COMOPT/          : OPTICAL PARAMETERS
CC            /COMBIO/          : BIOLOGICAL PARAMETERS
CC            /COMZEZM/         : DENSITY CRITERIA
CC
CC   OUTPUT :
CC   ------
CC      ARGUMENT
CC             FTRBIO   : GENERAL TRACER FORCING INCREASED BY THE
CC                                BIOLOGICAL FORCING 
CC         for plotting :
CC             FLXBIO   : INDIVIDUAL BIOLOGICAL FORCINGS (in mmoleN/m3/s)
CC             XLE      : LIGHT LIMITATION TERM
CC             XLT      : TEMPERATURE LIMITATION TERM
CC             XLM      : MIXING LIMITATION TERM
CC             XLNUT    : NUTRIENT LIMITATION TERM
CC
CC   WORKSPACE :
CC   ---------
CC      LOCAL
CC               GDEPT,ZDEPT                      : DEPTH AT LEVEL K
CC               ZDET,ZZOO,ZPHY,ZNUT              : NOW CONCENTRATIONS
CC               ZLM  ,ZLT,ZLNUT,ZLPE,ZLE         : TERMS FOR PHYTO GROWTH
CC               ZFLXNP,ZFLXPN,ZFLXPZ,ZFLXDZ      : FLUXES BETWEEN BIO BOXES
CC               ZFLXPD,ZFLXZD,ZFLXDN
CC               ZPHYA,ZZOOA,ZNUTA,ZDETA          : AFTER BIO TRENDS
CC               ZTESTP, ZTESTN, ZTESTZ, ZTESTD   : INTERMEDIATE VALUES FOR
CC                                                   POSSITIVITY TEST
CC               ZPHIMP, ZMP, ZPHIMZ, ZMZ         : MORTALITY TERMS
CC               ZPPZ, ZPDZ, ZPPPZ, ZPPDZ, ZFOOD  : PREFERENCES TERMS
CC               ZFILPZ, ZFILPD                   : FILTRATION TERMS
CC               INDICP, INDICZ, INDICN, INDICD   : INDICE FOR POSSITIVITY
CC
CC   EXTERNAL :                   NO
CC   --------
CC
CC   REFERENCES :                 NO
CC   ----------
CC   Levy, M., Memery, L. and Madec, G., 1997, The onset of a bloom in the
CC      NW Mediterranean sea : a 3D mesoscale process study with a primitive
CC      equation model, Journal of Marine System, accepted jan 97
CC
CC   MODIFICATIONS:
CC   --------------
CC       ORIGINAL : 95-05 (M. LEVY) for lodyc opa code
CC       MODIFICATION : 97-03 (M. LEVY) for lamont code

      include 'biology.h'
      include 'comm_para.h'

      REAL TRBIO(NPT,NZ,NTRBIO),T(NPT,NZ), H(NPT,NZ)
      REAL FTRBIO(NPT,NZ,NTRBIO), FLXBIO(NPT,NZ+1,NFLBIO)
      INTEGER NZI(NPT)
c
      REAL XZM(NPT), XZE(NPT), XPAR(NPT,NZ), XPARZE(NPT), XPARZM(NPT)

      REAL XLM(NPT,NZ),XLE(NPT,NZ),XLT(NPT,NZ),XLNUT(NPT,NZ)

      REAL GDEPT(MAXNZ)

      INTEGER KBIO(NPT)
C
C 0. INITIALISATION
C -----------------
C
              
      DO I=1,NPT
        ZDEPT=0.
        NZMAX= NZI(I)
        DO K=1, NZMAX
           ZDEPT=ZDEPT+H(I,K)
           GDEPT(K)=ZDEPT
        ENDDO
        DO K=1,NZMAX
          IF ( K.LE.JPBIO) THEN ! in the biologicaly active layer
C
C
C 1. INHIBITION BY DEEP MIXED-LAYER  - DETERMINATION OF PAR
C ---------------------------------------------------------
C
C TEMPERATURE EFFECT
C
           ZLT=1.
C
C DEEP MIXED LAYER REGIME 
C
           IF ( XZM(I) .GE. 2.*XZE(I) ) THEN
               IF(GDEPT(K).LE.XZE(I)) THEN
                 ZLM = 0.1
                 ZLE = ( 1. - EXP(-XPARZE(I)/AKI/ZLT) )
               ELSE
                 ZLM = 0.
                 ZLE = 0.
               ENDIF
           ELSE
C
C TRANSITION PHASE
C
              IF( ( ( 2.*XZE(I)-XZM(I) ).GE.0.).AND. 
     $            ( (    XZM(I)-XZE(I) ).GE.0.)        )  THEN 
                  IF(GDEPT(K).LE.XZE(I)) THEN
                    ZLM = 1. - 0.9*( (XZM(I)/XZE(I)) - 1. )
                    ZLE = ( 1. - EXP(-XPARZE(I)/AKI/ZLT) )
                  ELSE
                    ZLM = 0.
                    ZLE = 0.
                  ENDIF
              ELSE
C
C STRATIFICATION REGIME (ou si on supprime la limitation par le melange):
C
                ZLM = 1.
                  IF(GDEPT(K).LE.XZM(I)) THEN
                    ZLE = ( 1. - EXP(-XPARZM(I)/AKI/ZLT) )
                  ELSE
                    ZLE = ( 1. - EXP(-XPAR(I,K)/AKI/ZLT) )
                  ENDIF
C
              ENDIF
            ENDIF
C
C SAUVEGARDE
C
            XLE(I,K)=ZLE
            XLM(I,K)=ZLM
            XLT(I,K)=ZLT
C
C
C 2. TROPHIC VARIABLES( DET, ZOO, PHY, NUT)
C -----------------------------------------
C
C NEGATIVE TROPHIC VARIABLES DO NOT CONTRIBUTE TO THE FLUXES
C
            ZDET = MAX(0.,TRBIO(I,K,NDET))
            ZZOO = MAX(0.,TRBIO(I,K,NZOO))
            ZPHY = MAX(0.,TRBIO(I,K,NPHY))
            ZNUT = MAX(0.,TRBIO(I,K,NNUT))

C
C 3. NUTRIENT LIMITATION
C ----------------------
C
           ZLNUT = ZNUT / ( ZNUT + AKNUT )
C
C SAUVEGARDE
C
           XLNUT(I,K)=ZLNUT
C
C
C 4. SINKS AND SOURCES        
C --------------------
C
C
C 4.1 PHYTOPLANKTON PRODUCTION AND EXSUDATION
C
           ZFLXNP = TMUMAX * ZLE * ZLM * ZLNUT * ZLT * ZPHY
           ZFLXPN = RGAMMA * ZFLXNP
C
C 4.2 ZOOPLANKTON PRODUCTION
C
C
C PREFERENCES
C
            ZPPZ  = RPPZ
            ZPDZ  = 1. - RPPZ
            ZPPPZ = ( ZPPZ * ZPHY ) /
     $            ( ( ZPPZ * ZPHY + ZPDZ * ZDET ) + 1.E-13 )
            ZPPDZ = ( ZPDZ * ZDET ) /
     $            ( ( ZPPZ * ZPHY + ZPDZ * ZDET ) + 1.E-13 )
            ZFOOD = ZPPPZ * ZPHY + ZPPDZ * ZDET
C
C FILTRATION
C
            ZFILPZ = TAUS * ZPPPZ  / (AKS+ZFOOD)
            ZFILDZ = TAUS * ZPPDZ  / (AKS+ZFOOD)
C
C GRAZING  
C

            ZFLXPZ = ZFILPZ * ZPHY * ZZOO
            ZFLXDZ = ZFILDZ * ZDET * ZZOO
C
C 4.3 FECAL PELLETS PRODUCTION 
C
            ZFLXZL = RPNAZ * ZFLXPZ + RDNAZ * ZFLXDZ
C
C 4.4 ZOOPLANKTON LIQUIDE EXCRETION IF ZZOO GREATER THEN EGGZOO
C
            ZFLXZN = TAUZN * ZZOO * (1. + SIGN(1., ZZOO-EGGZOO))/2.
C
C 4.5 MORTALITY 
C
C PHYTOPLANKTON MORTALITY BY LACK OF NUTRIENTS 
C
            ZFLXPD = TMMINP * ZPHY
C 
C
C ZOOPLANKTON MORTALITY BY LACK OF FOOD IF ZZOO GREATER THEN EGGZOO
C
            ZFLXZD = TMMINZ * ZZOO * (1. + SIGN(1., ZZOO-EGGZOO))/2.
C
C
C 4.6 DETRITUS BREAKDOMN
C
            ZFLXDN = TAUDN * ZDET
 
C
C
C 5. DETERMINATION OF TRENDS
C --------------------------
C
C TOTAL TREND FOR EACH BIOLOGICAL TRACER
C
            ZPHYA = ZFLXNP - ZFLXPN - ZFLXPZ - ZFLXPD
            ZZOOA = ZFLXPZ + ZFLXDZ - ZFLXZD - ZFLXZL - ZFLXZN
            ZNUTA = ZFLXPN + ZFLXZN + ZFLXDN - ZFLXNP
            ZDETA = ZFLXPD + ZFLXZD + ZFLXZL - ZFLXDN - ZFLXDZ

C
            FLXBIO(I,K,1) = FLXBIO(I,K,1) + flor*ZFLXNP
            FLXBIO(I,K,2) = FLXBIO(I,K,2) + flor*ZFLXPZ
            FLXBIO(I,K,3) = FLXBIO(I,K,3) + flor*ZFLXPD
            FLXBIO(I,K,4) = FLXBIO(I,K,4) + flor*ZFLXDZ
            FLXBIO(I,K,5) = FLXBIO(I,K,5) + flor*(ZFLXZD+ZFLXZL)
            FLXBIO(I,K,6) = FLXBIO(I,K,6) + flor*ZFLXZN
            FLXBIO(I,K,7) = FLXBIO(I,K,7) + flor*ZFLXDN
            FLXBIO(I,K,9) = FLXBIO(I,K,9) + flor*(ZFLXPN + ZFLXZN + ZFLXDN)
C
C TRACER FORCING AT T-POINT ADDED TO THE GENERAL FORCING
C
            FTRBIO(I,K,NDET) =  FTRBIO(I,K,NDET) + ZDETA * H(I,K)
            FTRBIO(I,K,NZOO) =  FTRBIO(I,K,NZOO) + ZZOOA * H(I,K)
            FTRBIO(I,K,NPHY) =  FTRBIO(I,K,NPHY) + ZPHYA * H(I,K)
            FTRBIO(I,K,NNUT) =  FTRBIO(I,K,NNUT) + ZNUTA * H(I,K)
C
C

          ELSE      ! under the biological active layer ( k gt jpbio)
              RESTOR=TMINR + (TMAXR-TMINR)*
     $           (GDEPT(K)-GDEPT(JPBIO))/(GDEPT(NZI(I))-GDEPT(JPBIO))
              ZPHYA= - RESTOR * MAX(0.,TRBIO(I,K,NPHY)) 
              ZZOOA= - RESTOR * MAX(0.,TRBIO(I,K,NZOO)) 
              ZDETA= - RESTOR * MAX(0.,TRBIO(I,K,NDET)) 
              ZNUTA= -ZPHYA - ZZOOA - ZDETA
              FTRBIO(I,K,NDET) =  FTRBIO(I,K,NDET) + ZDETA * H(I,K)
              FTRBIO(I,K,NZOO) =  FTRBIO(I,K,NZOO) + ZZOOA * H(I,K)
              FTRBIO(I,K,NPHY) =  FTRBIO(I,K,NPHY) + ZPHYA * H(I,K)
              FTRBIO(I,K,NNUT) =  FTRBIO(I,K,NNUT) + ZNUTA * H(I,K)
              FLXBIO(I,K,9)=FLXBIO(I,K,9) + flor*ZNUTA
          ENDIF
C
C
C
C 7. DETRITUS SEDIMENTATION
C -------------------------
C
C WARNING : ZDETA DOES NOT HAVE THE SAME DIMENTION AS PREVIOUSLY, IT WAS
C IN mmole/m3/s AND NOW IT IS IN mmole/m2/s
C
         IF (K.EQ.1)  THEN
             ZDETA=-VSED*ZDET
         ELSE IF (K.EQ.NZI(I)) THEN
             ZDETA=VSED*MAX(0.,TRBIO(I,K-1,NDET))
         ELSE
             ZDETA=VSED*(MAX(0.,TRBIO(I,K-1,NDET))-ZDET)
         ENDIF
         FTRBIO(I,K,NDET) =  FTRBIO(I,K,NDET) + ZDETA
         FLXBIO(I,K,8)= FLXBIO(I,K,8) + flor* ZDETA / H(I,K)

      END DO     ! (end of k loop)



C END OF HORIZONTAL LOOP 
C ======================
C
      END DO
c
      RETURN
      END
CCC---------------------------------------------------------------------

      SUBROUTINE BIOINI_old(TRBIO,FTRBIO,NTRBIO,NFLBIO,NZ,NZI,NPT,KBIO,XZM,XZE,XPAR,
     >     XPARZE,XPARZM,FLXBIO,XLE,XLM,XLT,XLNUT)
CCC---------------------------------------------------------------------
CCC
CCC                       ROUTINE BIOINI
CCC                     *******************
CCC
CCC  PURPOSE :
CCC  ---------
CCC     INITIALIZE ALL ARRAYS USED IN THE BIOLOGICAL MODEL
CCC       modified by NN - tracers initialized in biotr_init
CCC
CCC
CC   METHOD :
CC   -------
CC
CC   INPUT :
CC   -----
CC      ARGUMENT
CC              NZ              : number of maximum vertical layers
CC              NZI             : number of vertical layers
CC              NPT             : number of horizontal grid points
CC              NTRBIO          : number of biological tracers 
CC              NFLBIO          : number of biological fluxes 
CC
CC      COMMON
CC            /COMOPT/          : OPTICAL PARAMETERS
CC            /COMBIO/          : BIOLOGICAL PARAMETERS
CC            /COMZEZM/         : DENSITY CRITERIA
CC
CC   OUTPUT :
CC   ------
CC      ARGUMENT
CC             TRBIO    : biological tracers array (in mmoleN/m3)
CC             FTRBIO   : biological tracer forcing (in mmoleN/m2/s)
CC             KBIO     : number of biologicaly active vertical layers
CC             FLXBIO   : INDIVIDUAL BIOLOGICAL FORCINGS (in mmoleN/m3/s)
CC             XLE      : LIGHT LIMITATION TERM
CC             XLT      : TEMPERATURE LIMITATION TERM
CC             XLM      : MIXING LIMITATION TERM
CC             XLNUT    : NUTRIENT LIMITATION TERM
CC             XZM      : MIXED LAYER DEPTH (in m)
CC             XZE      : EUPHOTIC LAYER DEPTH (in m)
CC             XPAR     : PAR AT EACH T LEVEL
CC             XPARZE   : MEAN PAR IN THE EUPHOTIC LAYER
CC             XPARZM   : MEAN PAR IN THE MIXED LAYER
CC
CC   WORKSPACE :
CC   ---------
CC      LOCAL
CC
CC   EXTERNAL :                   NO
CC   --------
CC
CC   REFERENCES :                 NO
CC   ----------
CC   Levy, M., Memery, L. and Madec, G., 1997, The onset of a bloom in the
CC      NW Mediterranean sea : a 3D mesoscale process study with a primitive
CC      equation model, Journal of Marine System, accepted jan 97
CC
CC   MODIFICATIONS:
CC   --------------
CC       ORIGINAL : 95-05 (M. LEVY) for lodyc opa code
CC       MODIFICATION : 97-03 (M. LEVY) for lamont code

      include 'biology.h'

      INTEGER NZI(NPT)
cml
      REAL TRBIO(NPT,NZ,NTRBIO), FTRBIO(NPT,NZ,NTRBIO), FLXBIO(NPT,NZ+1,NFLBIO)
cml
      REAL XZM(NPT), XZE(NPT), XPAR(NPT,NZ), XPARZE(NPT), XPARZM(NPT)
      REAL XLM(NPT,NZ),XLE(NPT,NZ),XLT(NPT,NZ),XLNUT(NPT,NZ)
      INTEGER KBIO(NPT)
C
C 0. INITIALISATION
C -----------------
C
      if (ibio_key .ne. 0) 
     *     open (unit = IBIO_OUT, file = f_bio(1:n_bio))


      DO I=1,NPT
        KBIO(I)=MIN(NZI(I),JPBIO)
        XZM(I)=0.
        XZE(I)=0.
        XPARZM(I)=0.
        XPARZE(I)=0.
      ENDDO


      DO I=1,NPT
        DO K=1,NZ
          XPAR(I,K)=0.
          XLM(I,K)=0.
          XLE(I,K)=0.
          XLT(I,K)=0.
          XLNUT(I,K)=0.
        ENDDO
      ENDDO
          
      DO I=1,NPT
        DO K=1,NZ
          DO NT=1,NTRBIO
            FTRBIO(I,K,NT)=0.
          ENDDO
        ENDDO
      ENDDO
      DO I=1,NPT
cml
        DO K=1,NZ+1
cml
           DO NT=1,NFLBIO
            FLXBIO(I,K,NT)=0.
          ENDDO
        ENDDO
      ENDDO

      RETURN
      END

      SUBROUTINE OPTIC_old(TRBIO,NTRBIO,H,NZ,NZI,NPT,RAU,QSR,
     $                 XZM, XZE, XPAR, XPARZE, XPARZM,
     $                 ZPARR, ZPARG, ZPAR100)
CCC---------------------------------------------------------------------
CCC
CCC                       ROUTINE OPTIC
CCC                     *******************
CCC
CCC  PURPOSE :
CCC  ---------
CCC     COMPUTE THE LIGHT PROPAGATION IN THE WATER COLUMN,
CCC     (PAR for Photosynthetic available radiation)
CCC     COMPUTE THE MIXED LAYER AND EUPHOTIC LAYER DEPTHS,
CCC     AND MEAN PAR IN THE EUPHOTIC AND MIXED LAYER
CCC
CCC
CC   METHOD :
CC   -------
CC
CC      LOCAL PAR IS COMPUTED IN W LAYERS USING LIGHT PROPAGATION
CC      MEAN PAR IN T LAYERS ARE COMPUTED BY INTEGRATION 
CC
CC
CC   INPUT :
CC   -----
CC      ARGUMENT
CC              TRBIO           : biological tracers array
CC              NTRBIO          : number of biological tracers
CC              H               : vertical scale factors at T points 
CC              NZ              : number of maximum vertical layers
CC              NZI             : number of vertical layers
CC              NPT             : number of horizontal grid points
CC              RAU             : density array (which one ? see with Naomi)
CC              QSR             : solar radiation flux
CC      COMMON
CC            /COMOPT/          : OPTICAL PARAMETERS
CC            /COMBIO/          : BIOLOGICAL PARAMETERS
CC            /COMZEZM/         : DENSITY CRITERIA
CC
CC   OUTPUT :
CC   ------
CC      ARGUMENT
CC              XZM             : MIXED LAYER DEPTH
CC              XZE             : EUPHOTIC LAYER DEPTH
CC              XPAR            : PAR AT EACH T LEVEL
CC              XPARZE          : MEAN PAR IN THE EUPHOTIC LAYER
CC              XPARZM          : MEAN PAR IN THE MIXED LAYER
CC
CC   WORKSPACE :
CC   ---------
CC      LOCAL     ZPARR         : RED COMPOUND OF PAR
CC                ZPARG         : GREEN COMPOUND OF PAR
CC                ZPAR0M        : IRRADIANCE JUST BELOW THE SURFACE
CC                ZPAR100       : IRRADIANCE AT EUPHOTIC LAYER DEPTH
CC                ZKR           : TOTAL ABSORPTION COEFFICIENT IN RED
CC                ZKG           : TOTAL ABSORPTION COEFFICIENT IN GREEN
CC                ZPIG          : TOTAL PIGMENT
CC
CC
CC   REFERENCES :                
CC   ----------
CC   Levy, M., Memery, L. and Madec, G., 1997, The onset of a bloom in the
CC      NW Mediterranean sea : a 3D mesoscale process study with a primitive
CC      equation model, Journal of Marine System, accepted jan 97
CC
CC   MODIFICATIONS:
CC   --------------
CC       ORIGINAL : 95-05 (M. LEVY) for lodyc opa code
CC       MODIFICATION : 97-03 (M. LEVY) for lamont code
C
CC----------------------------------------------------------------------
C
      include 'biology.h'

      REAL TRBIO(NPT,NZ,NTRBIO),RAU(NPT,NZ), QSR(NPT), H(NPT,NZ)
      INTEGER NZI(NPT)
      REAL ZPARR(NZ),ZPARG(NZ)
c
      REAL XZM(NPT), XZE(NPT), XPAR(NPT,NZ), XPARZE(NPT), XPARZM(NPT)
c
      REAL ZPAR0M, ZPAR100(NPT)
C
      DO I=1,NPT
C
C 1. DETERMINATION OF SURFACE IRRADIANCE AND PAR IN THE FIRST W LEVEL
C --------------------------------------------------------------------
C
        ZPAR0M = QSR(I)*0.43
        ZPAR100(I) = ZPAR0M*0.01
        ZPIG = MAX(0.,TRBIO(I,1,NPHY))*12*REDF/RCCHL/RPIG
        ZKR = XKR0+XKRP*ZPIG**XLR
        ZKG = XKG0+XKGP*ZPIG**XLG
        ZPARR(1) = 0.5* ZPAR0M 
        ZPARG(1) = 0.5* ZPAR0M
C
C 2. DETERMINATION OF XPAR
C ------------------------
C
C DETERMINATION OF LOCAL PAR IN W LEVELS
C
        DO K=2,NZI(I)
          ZPIG = MAX(0.,TRBIO(I,K,NPHY))*12*REDF/RCCHL/RPIG
          ZKR = XKR0+XKRP*ZPIG**XLR
          ZKG = XKG0+XKGP*ZPIG**XLG
          ZPARR(K) = ZPARR(K-1) *EXP( -ZKR*H(I,K-1) )
          ZPARG(K) = ZPARG(K-1) *EXP( -ZKG*H(I,K-1) )
        END DO

C
C MEAN PAR IN T LEVELS
        DO K=1,NZI(I)
          ZPIG = MAX(0.,TRBIO(I,K  ,NPHY))*12*REDF/RCCHL/RPIG
          ZKR = XKR0+XKRP*ZPIG**XLR
          ZKG = XKG0+XKGP*ZPIG**XLG
          ZPARR(K) =  ZPARR(K) / ZKR / H(I,K)
     >          * ( 1 - EXP( -ZKR*H(I,K) ) )  
          ZPARG(K) =  ZPARG(K) / ZKG / H(I,K)
     >          * ( 1 - EXP( -ZKG*H(I,K) ) )  
          XPAR(I,K) = MAX(ZPARR(K)+ZPARG(K),1.E-15)
        END DO
C
          
C
C 3. DETERMINATION OF  MIXED LAYER DEPTH AND ML MEAN PAR
C ------------------------------------------------------
C
C THE FIRST LAYER IS ALWAYS WITHIN THE MIXED LAYER
C
        XPARZM(I)=XPAR(I,1)*H(I,1)
        XZM(I)=H(I,1)
        DO K=2,NZI(I)
           IF ( (RAU(I,K)-RAU(I,1)).LE.DELTAR ) THEN
               XPARZM(I)=XPARZM(I) +XPAR(I,K)*H(I,K)
               XZM(I)=XZM(I)+H(I,K)
           ENDIF
        END DO
        XPARZM(I)=XPARZM(I)/XZM(I)
C
C
C 4. DETERMINATION OF EUPHOTIC LAYER DEPTH AND EL MEAN PAR
C --------------------------------------------------------
C
C THE FIRST LAYER IS ALWAYS WITHIN THE EUPHOTIC LAYER
C
        XPARZE(I)=XPAR(I,1)*H(I,1)
        XZE(I)=H(I,1)
        DO K=2,NZI(I)
            IF (XPAR(I,K).GE.ZPAR100(I)) THEN
               XPARZE(I)=XPARZE(I) +XPAR(I,K)*H(I,K)
               XZE(I)=XZE(I)+H(I,K)
            ENDIF
        END DO
        XPARZE(I)=XPARZE(I)/XZE(I)

C
      ENDDO 

      RETURN
      END

      subroutine bio_limit_old(nstep,trbio,ntrbio,nz,nzi,npt)
CCC---------------------------------------------------------------------
CCC
CCC                       ROUTINE BIO_LIMIT
CCC                     *******************
CCC
CCC  PURPOSE :
CCC  ---------
CCC     CHECK THE POSITIVITY OF BIOLOGICAL TRACERS
CCC
CC   METHOD :
CC   -------
CC
CC
CC   INPUT :
CC   -----
CC      ARGUMENT
CC              TRBIO           : biological tracers array (in mmoleN/m3)
CC              NTRBIO          : number of biological tracers
CC              NZ              : number of maximum vertical layers
CC              NZI             : number of vertical layers
CC              NPT             : number of horizontal grid points
CC
CC      COMMON
CC            /COMOPT/          : OPTICAL PARAMETERS
CC            /COMBIO/          : BIOLOGICAL PARAMETERS
CC            /COMZEZM/         : DENSITY CRITERIA
CC
CC   OUTPUT :
CC   ------
CC      ARGUMENT
CC              TRBIO           : biological tracers array (in mmoleN/m3)
CC
CC   WORKSPACE :
CC   ---------
CC      LOCAL
CC
CC   EXTERNAL :                   NO
CC   --------
CC
CC   REFERENCES :                 NO
CC   ----------
CC
CC   MODIFICATIONS:
CC   --------------
CC       MODIFICATION : 97-03 (M. LEVY) for lamont code
CC       MODIFICATION : 97-04 (Naomi) toss all negatives into last layer NUT

      include 'biology.h'

      REAL TRBIO(NPT,NZ,NTRBIO)
      INTEGER NZI(NPT)
C
C POSITIVITY 
C -------------
C

      INDICP = 0
      INDICZ = 0
      INDICN = 0
      INDICD = 0
      DTOT=0.
      TOT=0.
      
      DO I=1,NPT
         ZTOT = 0.
         NZBOT = NZI(I)
         DO K=1,NZBOT
            INDICP = INDICP + INT( ( 1. - SIGN (1.,TRBIO(I,K,NPHY)) )/2.)
            INDICZ = INDICZ + INT( ( 1. - SIGN (1.,TRBIO(I,K,NZOO)) )/2.)
            INDICN = INDICN + INT( ( 1. - SIGN (1.,TRBIO(I,K,NNUT)) )/2.)
            INDICD = INDICD + INT( ( 1. - SIGN (1.,TRBIO(I,K,NDET)) )/2.)
            DO N=1,NTRBIO
               ZTOT= ZTOT + TRBIO(I,K,N)-MAX(0., TRBIO(I,K,N) )
               TOT = TOT  + TRBIO(I,K,N)
               TRBIO(I,K,N) = MAX(0., TRBIO(I,K,N))
            ENDDO
         ENDDO

c**********************
c********WARNING*******
c  all negative tracer quantities are put into the last layer of nutrient
c**********************
         TRBIO(I,NZBOT,NNUT) = TRBIO(I,NZBOT,NNUT) + ZTOT

         DTOT = DTOT + ZTOT
      ENDDO
      
C
      if (ibio_key .ne. 0 .and. mod(nstep,nstep_bio_print).eq.0 ) then
         IF( INDICP+INDICZ+INDICN+INDICD.NE.0) THEN
            WRITE (IBIO_OUT,*) ' ===>>>> : NON CONSERVATIVE BIOLOGICAL SCHEME ' 
            WRITE (IBIO_OUT,*) ' =======                      ======='
            WRITE (IBIO_OUT,*) ' TIME STEP KT =', nstep
            WRITE (IBIO_OUT,*) ' INDICP =', INDICP
            WRITE (IBIO_OUT,*) ' INDICZ =', INDICZ
            WRITE (IBIO_OUT,*) ' INDICN =', INDICN
            WRITE (IBIO_OUT,*) ' INDICD =', INDICD
            WRITE (IBIO_OUT,*) ' DERIVE EN % DU TOTAL =', DTOT/TOT
         ENDIF
         call flush(IBIO_OUT)
      endif


C
      RETURN
      END

