Index: NEMO/trunk/src/OCE/TRA/traqsr.F90
===================================================================
--- NEMO/trunk/src/OCE/TRA/traqsr.F90	(revision 13165)
+++ NEMO/trunk/src/OCE/TRA/traqsr.F90	(revision 13205)
@@ -110,10 +110,9 @@
       REAL(wp) ::   zc0 , zc1 , zc2 , zc3    !    -         -
       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), ALLOCATABLE, DIMENSION(:,:)   :: zekb, zekg, zekr
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zetot, zchl3d
+      REAL(wp) ::   zz0 , zz1 , ze3t, zlui   !    -         -
+      REAL(wp) ::   zCb, zCmax, zpsi, zpsimax, zrdpsi, zCze
+      REAL(wp) ::   zlogc, zlogze, zlogCtot, zlogCze
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   :: ze0, ze1, ze2, ze3
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdt, zetot, ztmp3d
       !!----------------------------------------------------------------------
       !
@@ -160,75 +159,95 @@
       CASE( np_RGB , np_RGBc )         !==  R-G-B fluxes  ==!
          !
-         ALLOCATE( zekb(jpi,jpj)     , zekg(jpi,jpj)     , zekr  (jpi,jpj)     , &
-            &      ze0 (jpi,jpj,jpk) , ze1 (jpi,jpj,jpk) , ze2   (jpi,jpj,jpk) , &
-            &      ze3 (jpi,jpj,jpk) , zea (jpi,jpj,jpk) , zchl3d(jpi,jpj,jpk)   ) 
+         ALLOCATE( ze0 (jpi,jpj)           , ze1 (jpi,jpj) ,   &
+            &      ze2 (jpi,jpj)           , ze3 (jpi,jpj) ,   &
+            &      ztmp3d(jpi,jpj,nksr + 1)                     )
          !
          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
+            !
+            ! Separation in R-G-B depending on the surface Chl
+            ! perform and store as many of the 2D calculations as possible
+            ! before the 3D loop (use the temporary 2D arrays to replace the
+            ! most expensive calculations)
+            !
+            DO_2D_00_00
+                       ! zlogc = log(zchl)
+               zlogc = LOG ( MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) ) )     
+                       ! zc1 : log(zCze)  = log (1.12  * zchl**0.803)
+               zc1   = 0.113328685307 + 0.803 * zlogc                         
+                       ! zc2 : log(zCtot) = log(40.6  * zchl**0.459)
+               zc2   = 3.703768066608 + 0.459 * zlogc                           
+                       ! zc3 : log(zze)   = log(568.2 * zCtot**(-0.746))
+               zc3   = 6.34247346942  - 0.746 * zc2                           
+                       ! IF( log(zze) > log(102.) ) log(zze) = log(200.0 * zCtot**(-0.293))
+               IF( zc3 > 4.62497281328 ) zc3 = 5.298317366548 - 0.293 * zc2        
+               !   
+               ze0(ji,jj) = zlogc                                                 ! ze0 = log(zchl)
+               ze1(ji,jj) = EXP( zc1 )                                            ! ze1 = zCze
+               ze2(ji,jj) = 1._wp / ( 0.710 + zlogc * ( 0.159 + zlogc * 0.021 ) ) ! ze2 = 1/zdelpsi
+               ze3(ji,jj) = EXP( - zc3 )                                          ! ze3 = 1/zze
+            END_2D
+            
+!
+            DO_3D_00_00 ( 1, nksr + 1 )
+               ! zchl    = ALOG( ze0(ji,jj) )
+               zlogc = ze0(ji,jj)
+               !
+               zCb       = 0.768 + zlogc * ( 0.087 - zlogc * ( 0.179 + zlogc * 0.025 ) )
+               zCmax     = 0.299 - zlogc * ( 0.289 - zlogc * 0.579 )
+               zpsimax   = 0.6   - zlogc * ( 0.640 - zlogc * ( 0.021 + zlogc * 0.115 ) )
+               ! zdelpsi = 0.710 + zlogc * ( 0.159 + zlogc * 0.021 )
+               !
+               zCze   = ze1(ji,jj)
+               zrdpsi = ze2(ji,jj)                                                 ! 1/zdelpsi
+               zpsi   = ze3(ji,jj) * gdepw(ji,jj,jk,Kmm)                           ! gdepw/zze
+               !
+               ! NB. make sure zchl value is such that: zchl = MIN( 10. , MAX( 0.03, zchl ) )
+               zchl = MIN( 10. , MAX( 0.03, zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) * zrdpsi )**2 ) ) ) )
+               ! Convert chlorophyll value to attenuation coefficient look-up table index
+               ztmp3d(ji,jj,jk) = 41 + 20.*LOG10(zchl) + 1.e-15
+            END_3D
+         ELSE                                !* constant chlorophyll
+            zchl = 0.05
+            ! NB. make sure constant value is such that: 
+            zchl = MIN( 10. , MAX( 0.03, zchl ) )
+            ! Convert chlorophyll value to attenuation coefficient look-up table index
+            zlui = 41 + 20.*LOG10(zchl) + 1.e-15
             DO jk = 1, nksr + 1
-               DO jj = 2, jpjm1                       ! Separation in R-G-B depending of the surface Chl
-                  DO ji = 2, jpim1
-                     zchl    = MIN( 10. , MAX( 0.03, 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(ji,jj,jk,Kmm) / 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
+               ztmp3d(:,:,jk) = zlui 
             END DO
-         ELSE                                !* constant chrlorophyll
-           DO jk = 1, nksr + 1
-              zchl3d(:,:,jk) = 0.05 
-            ENDDO
          ENDIF
          !
          zcoef  = ( 1. - rn_abs ) / 3._wp    !* surface equi-partition in R-G-B
          DO_2D_00_00
-            ze0(ji,jj,1) = rn_abs * qsr(ji,jj)
-            ze1(ji,jj,1) = zcoef  * qsr(ji,jj)
-            ze2(ji,jj,1) = zcoef  * qsr(ji,jj)
-            ze3(ji,jj,1) = zcoef  * qsr(ji,jj)
-            zea(ji,jj,1) =          qsr(ji,jj)
+            ze0(ji,jj) = rn_abs * qsr(ji,jj)
+            ze1(ji,jj) = zcoef  * qsr(ji,jj)
+            ze2(ji,jj) = zcoef  * qsr(ji,jj)
+            ze3(ji,jj) = zcoef  * qsr(ji,jj)
+            ! store the surface SW radiation; re-use the surface ztmp3d array 
+            ! since the surface attenuation coefficient is not used
+            ztmp3d(ji,jj,1) =       qsr(ji,jj)
          END_2D
          !
-         DO jk = 2, nksr+1                   !* interior equi-partition in R-G-B depending of vertical profile of Chl
-            DO_2D_00_00
-               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_2D
-
-            DO_2D_00_00
-               zc0 = ze0(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * xsi0r       )
-               zc1 = ze1(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * zekb(ji,jj) )
-               zc2 = ze2(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * zekg(ji,jj) )
-               zc3 = ze3(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * zekr(ji,jj) )
-               ze0(ji,jj,jk) = zc0
-               ze1(ji,jj,jk) = zc1
-               ze2(ji,jj,jk) = zc2
-               ze3(ji,jj,jk) = zc3
-               zea(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * wmask(ji,jj,jk)
-            END_2D
-         END DO
+         !* interior equi-partition in R-G-B depending on vertical profile of Chl
+         DO_3D_00_00 ( 2, nksr + 1 )
+            ze3t = e3t(ji,jj,jk-1,Kmm)
+            irgb = NINT( ztmp3d(ji,jj,jk) )
+            zc0 = ze0(ji,jj) * EXP( - ze3t * xsi0r )
+            zc1 = ze1(ji,jj) * EXP( - ze3t * rkrgb(1,irgb) )
+            zc2 = ze2(ji,jj) * EXP( - ze3t * rkrgb(2,irgb) )
+            zc3 = ze3(ji,jj) * EXP( - ze3t * rkrgb(3,irgb) )
+            ze0(ji,jj) = zc0
+            ze1(ji,jj) = zc1
+            ze2(ji,jj) = zc2
+            ze3(ji,jj) = zc3
+            ztmp3d(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * wmask(ji,jj,jk)
+         END_3D
          !
          DO_3D_00_00( 1, nksr )
-            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zea(ji,jj,jk) - zea(ji,jj,jk+1) )
+            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( ztmp3d(ji,jj,jk) - ztmp3d(ji,jj,jk+1) )
          END_3D
          !
-         DEALLOCATE( zekb , zekg , zekr , ze0 , ze1 , ze2 , ze3 , zea , zchl3d ) 
+         DEALLOCATE( ze0 , ze1 , ze2 , ze3 , ztmp3d ) 
          !
       CASE( np_2BD  )            !==  2-bands fluxes  ==!
