MODULE trcopt
   !!======================================================================
   !!                         ***  MODULE trcopt  ***
   !! TOP :   LOBSTER Compute the light availability in the water column
   !!======================================================================
   !! History :    -   !  1995-05  (M. Levy) Original code
   !!              -   !  1999-09  (J.-M. Andre, M. Levy) 
   !!              -   !  1999-11  (C. Menkes, M.-A. Foujols) itabe initial
   !!              -   !  2000-02  (M.A. Foujols) change x**y par exp(y*log(x))
   !!             2.0  !  2007-12  (C. Deltel, G. Madec)  F90
   !!----------------------------------------------------------------------
#if defined key_lobster
   !!----------------------------------------------------------------------
   !!   'key_lobster'                                     LOBSTER bio-model
   !!----------------------------------------------------------------------
   !!   trc_opt        :   Compute the light availability in the water column
   !!----------------------------------------------------------------------
   USE oce_trc         !
   USE trc
   USE sms_lobster
   USE prtctl_trc      ! Print control for debbuging

   IMPLICIT NONE
   PRIVATE

   PUBLIC   trc_opt   ! called in trcprg.F90

   !!* Substitution
#  include "domzgr_substitute.h90"
   !!----------------------------------------------------------------------
   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
   !! $Id$ 
   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------

CONTAINS

   SUBROUTINE trc_opt( kt )
      !!---------------------------------------------------------------------
      !!                     ***  ROUTINE trc_opt  ***
      !!
      !! ** Purpose :   computes the light propagation in the water column
      !!              and the euphotic layer depth
      !!
      !! ** Method  :   local par is computed in w layers using light propagation
      !!              mean par in t layers are computed by integration
      !!---------------------------------------------------------------------
      INTEGER, INTENT( in ) ::   kt   ! index of the time stepping
      INTEGER  ::   ji, jj, jk
      INTEGER , DIMENSION(jpi,jpj)     ::   itabe           ! euphotic layer last k index
      INTEGER , DIMENSION(jpi,jpj,jpk) ::   imaske          ! euphotic layer mask
      REAL(wp) ::   zpig                                    ! total pigment
      REAL(wp) ::   zkr                                     ! total absorption coefficient in red
      REAL(wp) ::   zkg                                     ! total absorption coefficient in green
      REAL(wp), DIMENSION(jpi,jpj)     ::   zpar100         ! irradiance at euphotic layer depth
      REAL(wp), DIMENSION(jpi,jpj)     ::   zpar0m          ! irradiance just below the surface
      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zparr, zparg    ! red and green compound of par

      CHARACTER (len=25) :: charout
      !!---------------------------------------------------------------------

      IF( kt == nit000 ) THEN
         IF(lwp) WRITE(numout,*)
         IF(lwp) WRITE(numout,*) ' trc_opt: LOBSTER optic-model'
         IF(lwp) WRITE(numout,*) ' ~~~~~~~'
      ENDIF

      ! determination of surface irradiance
      ! -----------------------------------
      zpar0m (:,:)   = qsr   (:,:) * 0.43
      zpar100(:,:)   = zpar0m(:,:) * 0.01
      xpar   (:,:,1) = zpar0m(:,:)
      zparr  (:,:,1) = 0.5 * zpar0m(:,:)
      zparg  (:,:,1) = 0.5 * zpar0m(:,:)


      ! determination of xpar
      ! ---------------------

      DO jk = 2, jpk                     ! determination of local par in w levels
         DO jj = 1, jpj
            DO ji = 1, jpi
               zpig = MAX( TINY(0.), trn(ji,jj,jk-1,jpphy) ) * 12 * redf / rcchl / rpig
               zkr  = xkr0 + xkrp * EXP( xlr * LOG( zpig ) )
               zkg  = xkg0 + xkgp * EXP( xlg * LOG( zpig ) )
               zparr(ji,jj,jk) = zparr(ji,jj,jk-1) * EXP( -zkr * fse3t(ji,jj,jk-1) )
               zparg(ji,jj,jk) = zparg(ji,jj,jk-1) * EXP( -zkg * fse3t(ji,jj,jk-1) )
            END DO
        END DO
      END DO

      DO jk = 1, jpkm1                   ! mean par in t levels
         DO jj = 1, jpj
            DO ji = 1, jpi
               zpig = MAX( TINY(0.), trn(ji,jj,jk,jpphy) ) * 12 * redf / rcchl / rpig
               zkr  = xkr0 + xkrp * EXP( xlr * LOG( zpig ) )
               zkg  = xkg0 + xkgp * EXP( xlg * LOG( zpig ) )
               zparr(ji,jj,jk)    = zparr(ji,jj,jk) / zkr / fse3t(ji,jj,jk) * ( 1 - EXP( -zkr*fse3t(ji,jj,jk) ) )
               zparg(ji,jj,jk)    = zparg(ji,jj,jk) / zkg / fse3t(ji,jj,jk) * ( 1 - EXP( -zkg*fse3t(ji,jj,jk) ) )
               xpar (ji,jj,jk) = MAX( zparr(ji,jj,jk) + zparg(ji,jj,jk), 1.e-15 )
            END DO
         END DO
      END DO

      ! determination of euphotic layer depth (xze)
      ! -------------------------------------

      DO jk = 1, jpk                   ! imaske equal 1 in the euphotic layer, and 0 without
         DO jj = 1, jpj
            DO ji = 1,jpi
               IF( xpar(ji,jj,jk) >= zpar100(ji,jj) ) THEN
                  imaske(ji,jj,jk) = 1
               ELSE
                  imaske(ji,jj,jk) = 0
               ENDIF
            END DO
         END DO
      END DO
      !                                ! sum of imaske Cover the vertical with a minimim value of 1 
      itabe(:,:) = 1                   !   surface value setto 1 <=> set a ninimum value to 1
      DO jk = 2, jpk
         DO jj = 1, jpj
            DO ji = 1,jpi
               itabe(ji,jj) = itabe(ji,jj) + imaske(ji,jj,jk)
            END DO
         END DO
      END DO
      DO jj = 1, jpj                  ! converte the number of level into depth
         DO ji = 1,jpi
            xze(ji,jj) = fsdepw(ji,jj,itabe(ji,jj)+1)
         END DO
      END DO
      !

      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
         WRITE(charout, FMT="('opt')")
         CALL prt_ctl_trc_info(charout)
         CALL prt_ctl_trc( tab4d=trn, mask=tmask, clinfo=ctrcnm )
      ENDIF

   END SUBROUTINE trc_opt

#else
   !!======================================================================
   !!  Dummy module :                                   No PISCES bio-model
   !!======================================================================
CONTAINS
   SUBROUTINE trc_opt( kt )                   ! Empty routine
      INTEGER, INTENT( in ) ::   kt
      WRITE(*,*) 'trc_opt: You should not have seen this print! error?', kt
   END SUBROUTINE trc_opt
#endif 

   !!======================================================================
END MODULE  trcopt