SUBROUTINE BIOMASSE(TRBIO,FTRBIO,NTRBIO,NFLBIO,H,NZ,NZI,NPT,T,KBIO,
     *     XPAR, FLXBIO,XLE,XLT,XLN,flor)
CCC---------------------------------------------------------------------
CCC
ccc npzd model following fasham :
c     - no temperature effect on grazing 
c     - no temperature effect on detritus remineralisation
c     - zooplankton mortality do not go to detritus, but are remineralized
c        under the euphotic layer immediately, following the martin profile
c     - detritus stop sinking at jpbio
c     - tminr is the remineralization depth parameter
c     - tmaxr is the time damping parameter under ze (for phy, zoo and det) 
ccc modif in optics pur is used instead of par
c-------------------------------------------------------------------------


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

      REAL TRBIO(NPT,NZ,NTRBIO),T(NPT,NZ), H(NPT,NZ)
cml
      REAL FTRBIO(NPT,NZ,NTRBIO), FLXBIO(NPT,NZ+1,NFLBIO)
cml
      INTEGER NZI(NPT)
c
      REAL XPAR(NPT,NZ)

      REAL XLE(NPT,NZ),XLT(NPT,NZ),XLN(NPT,NZ)

      REAL GDEPT(MAXNZ)

      INTEGER KBIO(NPT)
C
C 0. INITIALISATION
C -----------------

      xhr = tminr
      xkpuropt=tmumax/alph
      facT = SLOPET**TOPTP
C
c P/I CURVE :
c
      DO I=1,NPT
         NZMAX = NZI(I)
         DO K=1,NZMAX
           IF ( K.LE.JPBIO) THEN ! in the biologicaly active layer
              XLT(I,K) = SLOPET**( T(I,K) - TOPTP)
              XKPUR= XKPUROPT * XLT(I,K)
              XLE(I,K) = 1. - EXP (- XPAR(I,K)/XKPUR) 
           ENDIF
         ENDDO
      ENDDO

      DO I=1,NPT
         ZDEPT=0.
         NZMAX = NZI(I)
         DO K=1,NZMAX
            ZDEPT=ZDEPT+H(I,K)
            GDEPT(K)=ZDEPT
         ENDDO
         ZDTA =0.    ! depth integrated implicit detritus flux
         ZE = gdept (jpbio)

         DO K=1,NZMAX
            IF ( K.LE.JPBIO) THEN ! in the biologicaly active layer
C
C 1. 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 2. LIMITATION TERMS OF PRODUCTION
C ---------------------------------------------------------
C
C  NUTRIENT LIMITATION
C
           ZLN = ZNUT / ( ZNUT + AKNUT )
C
C SAUVEGARDE
C
           XLN(I,K)=ZLN
C
C
C 4. SINKS AND SOURCES        
C --------------------
C
C
C 4.1 PHYTOPLANKTON PRODUCTION AND EXSUDATION
C
           ZFLXNP = TMUMAX * facT * XLT(I,K) * XLE(I,K) * ZLN * 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 
C
            ZFLXPD = TMMINP * ZPHY
C 
C
C ZOOPLANKTON MORTALITY 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 + ZFLXZL - ZFLXDN - ZFLXDZ

c total detritus flux :
            ZDTA = ZDTA + (ZFLXZD) * H(I,K)

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*(ZFLXZD-ZFLXDZ)  !new2
            FLXBIO(I,K,5) = FLXBIO(I,K,5) + flor*ZFLXZL           !new2
            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
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 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
                 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)
C
C
C 6. REMINERALIZATION
C --------------------
C
          ELSE      ! under the biological active layer ( k gt jpbio)
            xnp = zdta + vsed* MAX(0., TRBIO(I, JPBIO, NDET))
            fze = xnp
            ff = fze * ((gdept(k-1)/ze)**xhr)
            if (k. lt. nzmax) ff1= fze * ((gdept(k)/ze)**xhr)
            if (k. eq. nzmax) ff1= 0.
 
            RESTOR=tmaxr
            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 + (ff-ff1) / h(i,k)

            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 END OF LOOP FOR THE VERTICAL LAYERS
C ===========
C
        END DO

C
C
C

C END OF HORIZONTAL LOOP 
C ======================
C
      END DO

      RETURN
      END

      SUBROUTINE BIOINI(TRBIO,FTRBIO,NTRBIO,NFLBIO,NZ,NZI,NPT,KBIO,XPAR,
     >     FPAR,FLXBIO,XLE,XLT,XLN,XZE,XZM)
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             XLN    : NUTRIENT LIMITATION TERM
CC             XPAR     : PAR AT EACH T LEVEL
CC             FPAR     : PAR PENETRATION FUNCTION
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
CC       MODIFICATION : 97-07 (M. LEVY)  new optic 

      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), FPAR(NPT,NZ)
      REAL XLE(NPT,NZ),XLT(NPT,NZ),XLN(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)
        XZE(I)=0.
        XZM(I)=0.
        do k = jpbio+1,nzi(i)
           trbio(i,k,nphy) = 0.
           trbio(i,k,ndet) = 0.
           trbio(i,k,nzoo) = 0.
        enddo
        do k = nzi(i)+1,nz
           trbio(i,k,nphy) = -1000.
           trbio(i,k,ndet) = -1000.
           trbio(i,k,nzoo) = -1000.
           trbio(i,k,nnut) = -1000.
        enddo
      ENDDO

      DO I=1,NPT
        DO K=1,NZ
          XPAR(I,K)=0.
          FPAR(I,K)=0.
          XLE(I,K)=0.
          XLT(I,K)=0.
          XLN(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 PROPAG(TRBIO,NTRBIO,H,NZ,NZI,NPT,FPAR)
CCC---------------------------------------------------------------------
CCC
CCC                       ROUTINE PROPAG
CCC                     *******************
CCC
CCC  PURPOSE :
CCC  ---------
CCC     COMPUTE THE LIGHT PROPAGATION IN THE WATER COLUMN,
CCC     (PAR for Photosynthetic available radiation)
CCC     COMPUTE THE MEAN PROPAGATION IN T LEVELS
CCC
CC   CALL :
CC   -----
CC   Idealy every time step, but once a day is probably enough
CC
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      COMMON
CC            /COMOPT/          : OPTICAL PARAMETERS
CC            /COMBIO/          : BIOLOGICAL PARAMETERS
CC            /COMZEZM/         : DENSITY CRITERIA
CC
CC   OUTPUT :
CC   ------
CC      ARGUMENT
CC              FPAR            : mean fraction of incident par AT EACH T LEVEL
cc                                last modif : it is now in pur
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
CC       MODIFICATION : 97-12 (M. LEVY) add of xpar2pur and xkd_a
C
CC----------------------------------------------------------------------
C
      include 'biology.h'
      include 'comm_para.h'

      REAL TRBIO(NPT,NZ,NTRBIO),H(NPT,NZ)
      INTEGER NZI(NPT)
      REAL ZPARR(MAXNZ),ZPARG(MAXNZ)
c
      REAL FPAR(NPT,NZ)
C
      DO I=1,NPT

        KBIO = MIN( NZI(I), JPBIO)
C
C 1. DETERMINATION OF SURFACE IRRADIANCE AND PAR 
C ----------------------------------------------
C
        ZPARR(1) = 0.5
        ZPARG(1) = 0.5
C
C 2. DETERMINATION OF XPAR
C ------------------------
C
C DETERMINATION OF LOCAL PAR IN W LEVELS
C
        DO K=2,KBIO
          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,KBIO
          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) ) )  
          FPAR(I,K) = ZPARR(K)*XKD_A*XPAR2PUR_R+ZPARG(K)*XKD_A*XPAR2PUR_G
        END DO
      END DO
          
      RETURN
      END

      subroutine bio_limit(nstep,trbio,ntrbio,nz,nzi,npt,clm_coef,trclm,
     *                         inittr,h,flor,dnt,ftr_sp)
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              TRCLM           : tracer climatology
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),TRCLM(NPT,NZ,NTRBIO),h(NPT,NZ),
cml
     *          ftr_sp(npt,nz+1,ntrbio)
cml
      INTEGER NZI(NPT),INITTR(NTRBIO)
C
C SPONGE AT DEPTH
C -------------
C
cml
      c = flor/dnt
cml

      do itrac = 1, ntrbio
         if (inittr(itrac).eq.3) then
            do i = 1, npt
               do k = jpbio+1, nzi(i)
                  dtr = - clm_coef*(trbio(i,k,itrac)-trclm(i,k,itrac))
                  trbio(i,k,itrac) = trbio(i,k,itrac) + dtr
cml               c = flor*h(i,k)/dnt
                  ftr_sp(i,k,itrac) = ftr_sp(i,k,itrac) + c*dtr
               enddo
            enddo
         endif
      enddo

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**********************
         if (nzbot.gt.jpbio) then
            TRBIO(I,NZBOT,NNUT) = TRBIO(I,NZBOT,NNUT) + ZTOT
         endif

         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

      SUBROUTINE DAYLENGTH (NPT, tmonth, XLAT, DAYL)
CCC---------------------------------------------------------------------
CCC
CCC                       ROUTINE DAYLENGTH
CCC                     *******************
CCC
CCC  PURPOSE :
CCC  ---------
ccc  Computes the length of the day in hour as a function of latitude and
ccc  day of year (between 0 and 365)
ccc  Also computes declinaison and excentricity

CC  CALL :
CC  ------
CC  max once a day, but can be less
CC
cc  METHODE :
CC  ---------
CC    Compared to JMA routine, we have suppressed the following test :
CC    THIS TEST IS NEEDED FOR LATITUDE > 60N
CC
c         IF (abs(cosahsun).ge.1.) THEN
c           IF(dcln*zlat.gt.0.) THEN
c               dayl=24.
c           ELSE
c               dayl=0.
c           ENDIF
c         ELSE
c            ahsun=rtd*acos(cosahsun)
c            dayl =24.*ahsun*dtr/pi
c         ENDIF
c
CC   INPUT :
CC   -----
CC      ARGUMENT
CC      tmonth   : month of year between 0 and 12
CC      XLAT     : latitude

CC   OUTPUT :
CC   ------
CC      ARGUMENT
CC         DAYL  : lenght of day in hours (between 0 and 24)

CC
CC   REFERENCES :   from Jean-Michel Andre code (routine "astro")
CC   ----------
CC
CC 
CC   MODIFICATIONS:
CC   --------------
CC       ORIGINAL : 97-07 (M. LEVY) for lamont loam code
cc
cc---------------------------------------------------------------------------
        include 'biology.h'
        include 'comm_para.h'
cc---------------------------------------------------------------------------
        REAL XLAT(NPT)
        REAL DAYL(NPT)
cc---------------------------------------------------------------------------
C
C 0. INITIALISATION
C -----------------
C
        pi=3.14159265
        dtr=pi/180.
        rtd=1./dtr
        day = tmonth * 365. / 12.

C 1. EXCENTRICITY:
c ----------------
        xjr=(day-2)*2*pi/365.
        exc=1.000110+0.034221*cos(xjr)+0.00128*sin(xjr)
     >         +     0.000719*cos(2*xjr)+0.000077*sin(2*xjr)

c 2. DECLINAISON (dg):
C--------------------
        xjr=(day)*2*pi/366.

        dcln=0.32281-22.984*cos(xjr)-0.3499*cos(2*xjr)
     >                                -0.1398*cos(3*xjr)
     >                +3.7878*sin(xjr)+0.03205*sin(2*xjr)
     >                                +0.07187*sin(3*xjr)

c 3. Solar angle at end of day (dg) and day length (h) (test suppressed)
C ---------------------------------------------------
        DO I=1,NPT
          zlat=XLAT(I)
          cosahsun=-tan(zlat*dtr)*tan(dcln*dtr)
          ahsun=rtd*acos(cosahsun)
          IF (abs(cosahsun).ge.1.) THEN
            IF(dcln*zlat.gt.0.) THEN
                dayl(i)=24.
            ELSE
                dayl(i)=0.
            ENDIF
          ELSE
             ahsun=rtd*acos(cosahsun)
             dayl(i) =24.*ahsun*dtr/pi
          ENDIF

          dayl(i) = max(.001,dayl(i))

        ENDDO

        RETURN
        END


      SUBROUTINE PAR (NPT,NZ,KBIO,tmonth,DAYL,SW,FPAR,dtbio,XPAR)
CCC----------------------------------------------------------------------------
CCC
CCC                       ROUTINE PAR
CCC                     **************
CCC
CCC  PURPOSE :
CCC  ---------
ccc  Computes the diurnal cycle of the incoming PAR at the sea surface
ccc  given the daily averaged short wave radiation and the day length
ccc  Also computes the par in the water column, given the propagation function
CCC
CCC  METHODE :
CCC  ---------
CCC  Assumes a sin function for PAR0 
CCC  
CC
CC   CALL :
CC   ------
CC   to be called every bilogical time step ( has to resolve the diurnal cycle)
CC
CC      XHOUR    : hour of day (in hour)

CC   INPUT :
CC   -----
CC     ARGUMENT
CC      tmonth   : month of year (0 to 12)
CC      DAYL     : lenght of day in hours (between 0 and 24)
CC      SW       : daily averaged short wave radiation (in W/m2)
CC      FPAR     : par propagation fct in thw water column (between 0 and 1)
ccc                last modif : it is now pur
CC      dtbio    : time step in between two calls of this routine (in hour)

CC   OUTPUT :
CC   ------
CC      ARGUMENT
CC       XPAR  : par at each depth
cc               it is now pur

CC
CC   REFERENCES :   no
CC   ----------
CC
CC
CC   MODIFICATIONS:
CC   --------------
CC       ORIGINAL : 97-07 (M. LEVY) for lamont loam code
cc
cc---------------------------------------------------------------------------

        include 'biology.h'
        include 'comm_para.h'
cc---------------------------------------------------------------------------
        REAL XPAR(NPT,NZ)
        REAL FPAR(NPT,NZ)
        REAL SW(NPT)
        REAL DAYL(NPT)
        INTEGER KBIO(NPT)
cc---------------------------------------------------------------------------
C
C
C 0. INITIALISATION
C -----------------

        xday  = tmonth * 365./12.    ! day between 0 and 365
        xhour = 24.* (xday - int(xday)) ! hour of the day between 0 and 24
        

c zpar0 is defined as the derivative of time-integrated light (in order to
c be conservative) between hours xhourp and xhourm (that are forced to stay
c between 0 and 24, as xidt is non-continuous from one day to another)

        xhourp= min( xhour+ dtbio/2., 24.)
        xhourm= max( xhour- dtbio/2., 0.)
c#        delta=(xhourm-xhourp)
         
C
C 1. DIURNAL CYCLE OF PAR0
C ------------------------

        DO I=1,NPT
 
          dlen = dayl(i)
          ZPARDAY=24.0*0.43*SW(I)
          ZPAR0=ZPARDAY * ( XIDT(xhourp, dlen) - XIDT(xhourm, dlen) ) /dtbio

C
C 2. PAR IN THE WATER COLUMN
C --------------------------
          DO K=1,KBIO(I)
            XPAR(I,K)=FPAR(I,K)*ZPAR0
          ENDDO
        ENDDO
        
        RETURN
        END
      
C -----------------------------------------------------------------------------
        function xidt (xh, dayl)
C ------------------------------
C       FUNCTION XIDT
C
c--------------------------------------------------------------------------
c     purpose : computes the time integrated light field (dimensionless) over 
c      a day of length dayl, at the time xh (in hour)

c     xidt is the integral of sin ( pi2 * (xh-12+m) / dayl) between
c        12-m and xh, and normalized between 0 and 1
c     xidt is bounded between 0 and 1 (during the night)
c
c     M. Levy, aug 97 for lamont code
c----------------------------------------------------------------------------



      xm=dayl*0.5
      pi2=3.14159*0.5
      u= 0.5 *( 1- cos (pi2 * (xh-12+xm)/xm) )
      u= u* (xh.ge.(12-xm))
      xidt= u* (xh.le.(12+xm)) + 1 * (xh.gt.(12+xm))

      return
      end
      
C --------------------------
        subroutine ylat_degrees(npt,nxp,iox,ym,y_deg)
        dimension iox(npt),ym(1),y_deg(npt)
        
        do i=1,npt
           j=((iox(i)-1)/nxp)+1
           y_deg(i) = ym(j)
        enddo
        
        RETURN
        END


      SUBROUTINE ZEZM (DENS,FPAR,H,NZ,NZI,NPT,XZM,XZE)
CCC---------------------------------------------------------------------
CCC
CCC                       ROUTINE ZEZM
CCC                     *******************
CCC
CCC  PURPOSE :
CCC  ---------
CCC     COMPUTE THE MIXED LAYER AND EUPHOTIC LAYER DEPTHS,
CCC
CCC
CC   METHOD :
CC   -------
CC
CC      ZM computed on density gradient criterium
CC      ZE =1% of incident light
CC
CC
CC   INPUT :
CC   -----
CC      ARGUMENT
CC              DENS            : density 
CC              FPAR            : light attenuation function
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
CC   OUTPUT :
CC   ------
CC      ARGUMENT
CC              XZM             : MIXED LAYER DEPTH
CC              XZE             : EUPHOTIC LAYER DEPTH
CC
CC
CC   MODIFICATIONS:
CC   --------------
CC       ORIGINAL : 95-05 (M. LEVY) for lodyc opa code
CC       MODIFICATION : 97-07 (M. LEVY) for lamont code
C
CC----------------------------------------------------------------------
C
      include 'biology.h'

      REAL DENS(NPT,NZ), FPAR(NPT,NZ), H(NPT,NZ)
      INTEGER NZI(NPT)
c
      REAL XZM(NPT), XZE(NPT)
c
      DO I=1,NPT
C
C 1. DETERMINATION OF  MIXED LAYER DEPTH 
C --------------------------------------
C
C THE FIRST LAYER IS ALWAYS WITHIN THE MIXED LAYER
C
        XZM(I)=H(I,1)
        DO K=2,NZI(I)
           IF ( (DENS(I,K)-DENS(I,1)).LE.DELTAR ) THEN
               XZM(I)=XZM(I)+H(I,K)
           ENDIF
        END DO
C
C
C 2. DETERMINATION OF EUPHOTIC LAYER DEPTH AND EL MEAN PAR
C --------------------------------------------------------
C
C THE FIRST LAYER IS ALWAYS WITHIN THE EUPHOTIC LAYER
C
        XZE(I)=H(I,1)
        DO K=2,NZI(I)
            IF (FPAR(I,K).GE.0.01) THEN
               XZE(I)=XZE(I)+H(I,K)
            ENDIF
        END DO
C
      ENDDO 

      RETURN
      END