Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 2325)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 2326)
@@ -103,5 +103,5 @@
       REAL(wp), DIMENSION(jpi,jpj)     ::   zekb, zekg, zekr            ! 2D workspace
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze0, ze1 , ze2, ze3, zea    ! 3D workspace
-      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ztrdt,  ztrds
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ztrdt
       !!----------------------------------------------------------------------
 
@@ -115,5 +115,4 @@
       IF( l_trdtra ) THEN      ! Save ta and sa trends
          ALLOCATE( ztrdt(jpi,jpj,jpk) )   ;    ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
-         ALLOCATE( ztrds(jpi,jpj,jpk) )   ;    ztrds(:,:,:) = 0.
       ENDIF
 
@@ -143,7 +142,12 @@
          !                                        ! ============================================== !
          DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
+            qsr_hc(:,:,jk) = ro0cpr * ( etot3(:,:,jk) - etot3(:,:,jk+1) )
+         END DO
+         !                                        Add to the general trend
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1 
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  qsr_hc(ji,jj,jk) =  ro0cpr * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) ) / fse3t(ji,jj,jk) 
+                  z1_e3t = zfact / fse3t(ji,jj,jk)
+                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) * z1_e3t
                END DO
             END DO
@@ -234,5 +238,5 @@
                         zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r )
                         zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r )
-                        qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0*tmask(ji,jj,jk) - zc1*tmask(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
+                        qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0*tmask(ji,jj,jk) - zc1*tmask(ji,jj,jk+1) ) 
                      END DO
                   END DO
@@ -251,15 +255,15 @@
          ENDIF
          !
-      ENDIF
-      !                                        Add to the general trend
-      !                                        ************************
-      DO jk = 1, nksr
-         DO jj = 2, jpjm1 
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               z1_e3t = zfact / fse3t(ji,jj,jk)
-               tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) * z1_e3t
+         !                                        Add to the general trend
+         DO jk = 1, nksr
+            DO jj = 2, jpjm1 
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  z1_e3t = zfact / fse3t(ji,jj,jk)
+                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) * z1_e3t
+               END DO
             END DO
          END DO
-      END DO
+         !
+      ENDIF
       !
       IF( lrst_oce ) THEN   !                  Write in the ocean restart file
@@ -276,6 +280,5 @@
          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_qsr, ztrdt )
-         CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_qsr, ztrds )
-         DEALLOCATE( ztrdt )    ;        DEALLOCATE( ztrds )
+         DEALLOCATE( ztrdt )
       ENDIF
       !                       ! print mean trends (used for debugging)
@@ -443,5 +446,5 @@
                   !
                   DO jk = 1, nksr
-                     etot3(:,:,jk) = ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) ) / fse3t_0(:,:,jk)
+                     etot3(:,:,jk) = ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) ) 
                   END DO
                   etot3(:,:,nksr+1:jpk) = 0.e0                ! below 400m set to zero
@@ -471,5 +474,5 @@
                         zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r )
                         zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r )
-                        etot3(ji,jj,jk) = (  zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1)  ) / fse3t_0(ji,jj,jk)
+                        etot3(ji,jj,jk) = (  zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1)  ) 
                      END DO
                   END DO
Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/trc_oce.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/trc_oce.F90	(revision 2325)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/trc_oce.F90	(revision 2326)
@@ -168,14 +168,9 @@
       REAL(wp), DIMENSION(3,61), INTENT(out) ::   prgb   ! tabulated attenuation coefficient
       !!
-      INTEGER  ::   jchl, jband   ! dummy loop indices
+      INTEGER  ::   jc, jb ! dummy loop indice
+      INTEGER  ::   irgb   ! temporary integer
+      REAL(wp) ::   zchl   ! temporary scalar
       INTEGER  ::   numlight
-      REAL(wp) ::   zchl
-      !!----------------------------------------------------------------------
-      !
-      CALL ctl_opn( numlight, 'kRGB61.txt', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
-      DO jchl = 1, 61
-         READ(numlight,*) zchl, ( prgb(jband,jchl), jband=1,3 )
-      END DO
-      CLOSE( numlight )
+      !!----------------------------------------------------------------------
       !
       IF(lwp) THEN                         ! control print
@@ -183,5 +178,21 @@
          WRITE(numout,*) ' trc_oce_rgb_read : optical look-up table read in kRGB61.txt file'
          WRITE(numout,*) ' ~~~~~~~~~~~~~~~~'
+         WRITE(numout,*) 
       ENDIF
+      !
+      CALL ctl_opn( numlight, 'kRGB61.txt', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
+      DO jc = 1, 61
+         READ(numlight,*) zchl, ( prgb(jb,jc), jb = 1, 3 )
+         irgb = NINT( 41 + 20.* LOG10( zchl ) + 1.e-15 )   
+         IF(lwp) WRITE(numout,*) '    jc =', jc, '  Chl = ', zchl, '  irgb = ', irgb  
+         IF( irgb /= jc ) THEN  
+            IF(lwp) WRITE(numout,*) '    jc =', jc, '  Chl = ', zchl, '  Chl class = ', irgb
+            CALL ctl_stop( 'trc_oce_rgb_read : inconsistency in Chl tabulated attenuation coeff.' )
+         ENDIF
+      END DO
+      CLOSE( numlight )
+      !
+      r_si2 = 1.e0 / prgb(1, 1)      ! blue with the smallest chlorophyll concentration)
+      IF(lwp) WRITE(numout,*) '      RGB longest depth of extinction    r_si2 = ', r_si2
       !
    END SUBROUTINE trc_oce_rgb_read
Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zopt.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zopt.F90	(revision 2325)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zopt.F90	(revision 2326)
@@ -234,14 +234,13 @@
       !!----------------------------------------------------------------------
 
-!!      CALL trc_oce_rgb( xkrgb )                  ! tabulated attenuation coefficients
-      CALL trc_oce_rgb_read( xkrgb )               ! tabulated attenuation coefficients
+      CALL trc_oce_rgb( xkrgb )                  ! tabulated attenuation coefficients
+!!      CALL trc_oce_rgb_read( xkrgb )               ! tabulated attenuation coefficients
       nksrp = trc_oce_ext_lev( r_si2, 0.33e2 )     ! max level of light extinction (Blue Chl=0.01)
       IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksrp, ' ref depth = ', gdepw_0(nksrp+1), ' m'
       !
-      etot (:,:,:) = 0.e0
-      enano(:,:,:) = 0.e0
-      ediat(:,:,:) = 0.e0
+                         etot (:,:,:) = 0.e0
+                         enano(:,:,:) = 0.e0
+                         ediat(:,:,:) = 0.e0
       IF( ln_qsr_bio )   etot3(:,:,:) = 0.e0
-
       ! 
    END SUBROUTINE p4z_opt_init
