Index: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 6140)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 6403)
@@ -11,4 +11,5 @@
    !!            3.2  !  2009-04  (G. Madec & NEMO team) 
    !!            3.6  !  2012-05  (C. Rousset) store attenuation coef for use in ice model 
+   !!            3.6  !  2015-12  (O. Aumont, J. Jouanno, C. Ethe) use vertical profile of chlorophyll
    !!            3.7  !  2015-11  (G. Madec, A. Coward)  remove optimisation for fix volume 
    !!----------------------------------------------------------------------
@@ -100,4 +101,5 @@
       !! Reference  : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
       !!              Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516.
+      !!              Morel, A. et Berthon, JF, 1989, Limnol Oceanogr 34(8), 1545-1562
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt     ! ocean time-step
@@ -109,7 +111,9 @@
       REAL(wp) ::   zzc0, zzc1, zzc2, zzc3   !    -         -
       REAL(wp) ::   zz0 , zz1                !    -         -
+      REAL(wp) ::   zCb, zCmax, zze, zpsi, zpsimax, zdelpsi, zCtot, zCze
+      REAL(wp) ::   zlogc, zlogc2, zlogc3 
       REAL(wp), POINTER, DIMENSION(:,:)   :: zekb, zekg, zekr
       REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zetot
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zetot, zchl3d
       !!----------------------------------------------------------------------
       !
@@ -158,28 +162,35 @@
          !
          CALL wrk_alloc( jpi,jpj,       zekb, zekg, zekr        ) 
-         CALL wrk_alloc( jpi,jpj,jpk,   ze0, ze1, ze2, ze3, zea ) 
+         CALL wrk_alloc( jpi,jpj,jpk,   ze0, ze1, ze2, ze3, zea, zchl3d ) 
          !
          IF( nqsr == np_RGBc ) THEN          !*  Variable Chlorophyll
             CALL fld_read( kt, 1, sf_chl )         ! Read Chl data and provides it at the current time step
-            DO jj = 2, jpjm1                       ! Separation in R-G-B depending of the surface Chl
-               DO ji = fs_2, fs_jpim1
-                  zchl = MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) )
-                  irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
-                  zekb(ji,jj) = rkrgb(1,irgb)
-                  zekg(ji,jj) = rkrgb(2,irgb)
-                  zekr(ji,jj) = rkrgb(3,irgb)
+            DO jk = 1, nksr + 1
+               DO jj = 2, jpjm1                       ! Separation in R-G-B depending of the surface Chl
+                  DO ji = fs_2, fs_jpim1
+                     zchl    = sf_chl(1)%fnow(ji,jj,1)
+                     zCtot   = 40.6  * zchl**0.459
+                     zze     = 568.2 * zCtot**(-0.746)
+                     IF( zze > 102. ) zze = 200.0 * zCtot**(-0.293)
+                     zpsi    = gdepw_n(ji,jj,jk) / zze
+                     !
+                     zlogc   = LOG( zchl )
+                     zlogc2  = zlogc * zlogc
+                     zlogc3  = zlogc * zlogc * zlogc
+                     zCb     = 0.768 + 0.087 * zlogc - 0.179 * zlogc2 - 0.025 * zlogc3
+                     zCmax   = 0.299 - 0.289 * zlogc + 0.579 * zlogc2
+                     zpsimax = 0.6   - 0.640 * zlogc + 0.021 * zlogc2 + 0.115 * zlogc3
+                     zdelpsi = 0.710 + 0.159 * zlogc + 0.021 * zlogc2
+                     zCze    = 1.12  * (zchl)**0.803 
+                     !
+                     zchl3d(ji,jj,jk) = zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) / zdelpsi )**2 ) )
+                  END DO
+                  !
                END DO
             END DO
          ELSE                                !* constant chrlorophyll
-            zchl = 0.05                            ! constant chlorophyll
-            !                                      ! Separation in R-G-B depending of the chlorophyll
-            irgb = NINT( 41 + 20.*LOG10( zchl ) + 1.e-15 )
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1
-                  zekb(ji,jj) = rkrgb(1,irgb)                      
-                  zekg(ji,jj) = rkrgb(2,irgb)
-                  zekr(ji,jj) = rkrgb(3,irgb)
-               END DO
-            END DO
+           DO jk = 1, nksr + 1
+              zchl3d(:,:,jk) = 0.05 
+            ENDDO
          ENDIF
          !
@@ -195,5 +206,15 @@
          END DO
          !
-         DO jk = 2, nksr+1                   !* interior equi-partition in R-G-B
+         DO jk = 2, nksr+1                   !* interior equi-partition in R-G-B depending of vertical profile of Chl
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  zchl = MIN( 10. , MAX( 0.03, zchl3d(ji,jj,jk) ) )
+                  irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
+                  zekb(ji,jj) = rkrgb(1,irgb)
+                  zekg(ji,jj) = rkrgb(2,irgb)
+                  zekr(ji,jj) = rkrgb(3,irgb)
+               END DO
+            END DO
+
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
@@ -220,5 +241,5 @@
          !
          CALL wrk_dealloc( jpi,jpj,        zekb, zekg, zekr        ) 
-         CALL wrk_dealloc( jpi,jpj,jpk,   ze0, ze1, ze2, ze3, zea ) 
+         CALL wrk_dealloc( jpi,jpj,jpk,   ze0, ze1, ze2, ze3, zea, zchl3d ) 
          !
       CASE( np_2BD  )            !==  2-bands fluxes  ==!
