Index: trunk/NEMO/OPA_SRC/TRA/traqsr.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 296)
+++ trunk/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 457)
@@ -4,5 +4,11 @@
    !! Ocean physics: solar radiation penetration in the top ocean levels
    !!======================================================================
-
+   !! History :
+   !!   6.0  !  90-10  (B. Blanke)  Original code
+   !!   7.0  !  91-11  (G. Madec)
+   !!        !  96-01  (G. Madec)  s-coordinates
+   !!   8.5  !  02-06  (G. Madec)  F90: Free form and module
+   !!   9.0  !  04-08  (C. Talandier) New trends organization
+   !!   9.0  !  05-11  (G. Madec) zco, zps, sco coordinate
    !!----------------------------------------------------------------------
    !!   tra_qsr      : trend due to the solar radiation penetration
@@ -81,14 +87,8 @@
       !!              - save the trend in ttrd ('key_trdtra')
       !!
-      !! History :
-      !!   6.0  !  90-10  (B. Blanke)  Original code
-      !!   7.0  !  91-11  (G. Madec)
-      !!        !  96-01  (G. Madec)  s-coordinates
-      !!   8.5  !  02-06  (G. Madec)  F90: Free form and module
-      !!   9.0  !  04-08  (C. Talandier) New trends organization
       !!----------------------------------------------------------------------
       !! * Modules used     
-      USE oce, ONLY :    ztdta => ua,   & ! use ua as 3D workspace   
-                         ztdsa => va      ! use va as 3D workspace   
+      USE oce, ONLY :    ztrdt => ua,   & ! use ua as 3D workspace   
+                         ztrds => va      ! use va as 3D workspace   
 
       !! * Arguments
@@ -97,7 +97,5 @@
       !! * Local declarations
       INTEGER ::    ji, jj, jk            ! dummy loop indexes
-      REAL(wp) ::   zc0, zta              ! temporary scalars
-      REAL(wp) ::   zc1 , zc2 ,        &  ! temporary scalars
-                    zdp1, zdp2            !
+      REAL(wp) ::   zc0 , zta             ! temporary scalars
       !!----------------------------------------------------------------------
 
@@ -106,22 +104,24 @@
          IF ( lwp )  WRITE(numout,*) 'tra_qsr : penetration of the surface solar radiation'
          IF ( lwp )  WRITE(numout,*) '~~~~~~~'
+         CALL tra_qsr_init
       ENDIF
 
       ! Save ta and sa trends
       IF( l_trdtra )   THEN
-         ztdta(:,:,:) = ta(:,:,:) 
-         ztdsa(:,:,:) = 0.e0
-      ENDIF
-
-      IF( lk_qsr_sms .AND. ln_qsr_sms ) THEN    !  Biological fluxes  !
-         !                                      ! =================== !
-         !
-         !                                                ! ===============
-         DO jk = 1, jpkm1                                 ! Horizontal slab
-            !                                             ! ===============
+         ztrdt(:,:,:) = ta(:,:,:) 
+         ztrds(:,:,:) = 0.e0
+      ENDIF
+
+      ! ---------------------------------------------- !
+      !  Biological fluxes  : all vertical coordinate  !
+      ! ---------------------------------------------- !
+      IF( lk_qsr_sms .AND. ln_qsr_sms ) THEN      
+         !                                                   ! ===============
+         DO jk = 1, jpkm1                                    ! Horizontal slab
+            !                                                ! ===============
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   
-                  zc0 = ro0cpr  / fse3t(ji,jj,jk)      ! compute the qsr trend
+                  zc0 = ro0cpr  / fse3t(ji,jj,jk)         ! compute the qsr trend
                   zta = zc0 * ( etot3(ji,jj,jk  ) * tmask(ji,jj,jk)     &
                      &        - etot3(ji,jj,jk+1) * tmask(ji,jj,jk+1) )
@@ -131,110 +131,60 @@
                END DO
             END DO
-            !                                             ! ===============
-         END DO                                           !   End of slab
-         !                                                ! ===============
-         ! save the trends for diagnostic
-         ! qsr tracers trends
-         IF( l_trdtra )   THEN
-            ztdta(:,:,:) = ta(:,:,:) - ztdta(:,:,:)
-            CALL trd_mod(ztdta, ztdsa, jpttdqsr, 'TRA', kt)
-         ENDIF
-
+            !                                                ! ===============
+         END DO                                              !   End of slab
+         !                                                   ! ===============
+
+      ! ---------------------------------------------- !
+      !  Ocean alone :
+      ! ---------------------------------------------- !
       ELSE
          !                                                ! =================== !
-         IF( lk_sco ) THEN                                !    s-coordinate     !
+         IF( ln_sco ) THEN                                !    s-coordinate     !
+            !                                             ! =================== !
+            DO jk = 1, jpkm1
+               ta(:,:,jk) = ta(:,:,jk) + etot3(:,:,jk) * qsr(:,:)
+            END DO
+         ENDIF
          !                                                ! =================== !
-         !
-         !                                                   ! ===============
-            DO jk = 1, jpkm1                                 ! Horizontal slab
-            !                                                ! ===============
+         IF( ln_zps ) THEN                                !    partial steps    !
+            !                                             ! =================== !
+            DO jk = 1, nksr
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1   ! vector opt.
-
-                     zdp1 = -fsdepw(ji,jj,jk  )              ! compute the qsr trend
-                     zdp2 = -fsdepw(ji,jj,jk+1)
-                     zc0  = qsr(ji,jj) * ro0cpr / fse3t(ji,jj,jk)
-                     zc1  =   (  rabs * EXP(zdp1/xsi1) + (1.-rabs) * EXP(zdp1/xsi2)  )
-                     zc2  = - (  rabs * EXP(zdp2/xsi1) + (1.-rabs) * EXP(zdp2/xsi2)  )
-                     zta  = zc0 * (  zc1 * tmask(ji,jj,jk) + zc2 * tmask(ji,jj,jk+1)  )
-                     
-                     ta(ji,jj,jk) = ta(ji,jj,jk) + zta       ! add qsr trend to the temperature trend
-                     
+                     ! qsr trend from gdsr
+                     zc0 = qsr(ji,jj) / fse3t(ji,jj,jk)
+                     zta = zc0 * ( gdsr(jk) * tmask(ji,jj,jk) - gdsr(jk+1) * tmask(ji,jj,jk+1) )
+                     ! add qsr trend to the temperature trend
+                     ta(ji,jj,jk) = ta(ji,jj,jk) + zta
                   END DO
                END DO
-               !                                                ! ===============
-            END DO                                              !   End of slab
-            !                                                   ! ===============
-            ! save the trends for diagnostic
-            ! qsr tracers trends
-            IF( l_trdtra )   THEN
-               ztdta(:,:,:) = ta(:,:,:) - ztdta(:,:,:)
-               CALL trd_mod(ztdta, ztdsa, jpttdqsr, 'TRA', kt)
-            ENDIF
-            !
+            END DO
          ENDIF
          !                                                ! =================== !
-         IF( lk_zps ) THEN                                !    partial steps    !
+         IF( ln_zco ) THEN                                !     z-coordinate    !
             !                                             ! =================== !
-            !
-            !                                                ! ===============
-            DO jk = 1, nksr                                  ! Horizontal slab
-               !                                             ! ===============
+            DO jk = 1, nksr
+               zc0 = 1. / e3t_0(jk)
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1   ! vector opt.
-                     
-                     zc0 = qsr(ji,jj) / fse3t(ji,jj,jk)      ! compute the qsr trend
-                     zta = zc0 * ( gdsr(jk) * tmask(ji,jj,jk) - gdsr(jk+1) * tmask(ji,jj,jk+1) )
-                     
-                     ta(ji,jj,jk) = ta(ji,jj,jk) + zta       ! add qsr trend to the temperature trend
-                     
+                     ! qsr trend
+                     zta = qsr(ji,jj) * zc0 * ( gdsr(jk)*tmask(ji,jj,jk) - gdsr(jk+1)*tmask(ji,jj,jk+1) )
+                     ! add qsr trend to the temperature trend
+                     ta(ji,jj,jk) = ta(ji,jj,jk) + zta      
                   END DO
                END DO
-               !                                             ! ===============
-            END DO                                           !   End of slab
-            !                                                ! ===============
-            ! save the trends for diagnostic
-            ! qsr tracers trends
-            IF( l_trdtra )   THEN
-               ztdta(:,:,:) = ta(:,:,:) - ztdta(:,:,:)
-               CALL trd_mod(ztdta, ztdsa, jpttdqsr, 'TRA', kt)
-            ENDIF
-            !
-         ENDIF
-         !                                                ! =================== !
-         IF( lk_zco ) THEN                                !     z-coordinate    !
-            !                                             ! =================== !
-            !
-            !                                                ! ===============
-            DO jk = 1, nksr                                  ! Horizontal slab
-               !                                             ! ===============
-               zc0 = 1. / fse3t(1,1,jk)
-               DO jj = 2, jpjm1
-                  DO ji = fs_2, fs_jpim1   ! vector opt.
-                     !                                       ! compute qsr forcing trend
-                     zta = qsr(ji,jj) * zc0 * ( gdsr(jk)*tmask(ji,jj,jk) - gdsr(jk+1)*tmask(ji,jj,jk+1) )
-                     
-                     ta(ji,jj,jk) = ta(ji,jj,jk) + zta       ! add qsr trend to the temperature trend
-                     
-                  END DO
-               END DO
-               !                                             ! ===============
-            END DO                                           !   End of slab
-            !                                                ! ===============
-            ! save the trends for diagnostic
-            ! qsr tracers trends
-            IF( l_trdtra )   THEN
-               ztdta(:,:,:) = ta(:,:,:) - ztdta(:,:,:)
-               CALL trd_mod(ztdta, ztdsa, jpttdqsr, 'TRA', kt)
-            ENDIF
-            !
+            END DO
          ENDIF
          !
       ENDIF
 
-
-      IF(ln_ctl) THEN         ! print mean trends (used for debugging)
-         CALL prt_ctl(tab3d_1=ta, clinfo1=' qsr  - Ta: ', mask1=tmask, clinfo3='tra-ta')
-      ENDIF
+      ! qsr tracers trends saved the trends for diagnostics
+      IF( l_trdtra )   THEN
+         ztrdt(:,:,:) = ta(:,:,:) - ztrdt(:,:,:)
+         CALL trd_mod( ztrdt, ztrds, jpttdqsr, 'TRA', kt )
+      ENDIF
+
+      !                       ! print mean trends (used for debugging)
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' qsr  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
 
    END SUBROUTINE tra_qsr
@@ -258,12 +208,9 @@
       !! Reference :
       !!   Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
-      !!
-      !! History :
-      !!   8.5  !  02-06  (G. Madec) Original code
-      !!----------------------------------------------------------------------
-      !! * Local declarations
+      !!----------------------------------------------------------------------
       INTEGER ::    ji,jj,jk, &  ! dummy loop index
                     indic        ! temporary integer
-      REAL(wp) ::   zdp1         ! temporary scalar
+      REAL(wp) ::   zc0 , zc1 , zc2 ,   & ! temporary scalars
+         &          zcst, zdp1, zdp2      !    "         "
 
       NAMELIST/namqsr/ ln_traqsr, rabs, xsi1, xsi2, ln_qsr_sms
@@ -278,23 +225,22 @@
       ! ---------------------------
       IF( ln_traqsr  ) THEN
-        IF ( lwp ) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'tra_qsr_init : penetration of the surface solar radiation'
-         WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*) '    Namelist namqsr : set the parameter of penetration'
-         WRITE(numout,*) '        fraction associated with xsi     rabs        = ',rabs
-         WRITE(numout,*) '        first depth of extinction        xsi1        = ',xsi1
-         WRITE(numout,*) '        second depth of extinction       xsi2        = ',xsi2
-         IF( lk_qsr_sms ) THEN
-            WRITE(numout,*) '     Biological fluxes for light(Y/N) ln_qsr_sms  = ',ln_qsr_sms
-         ENDIF
-         WRITE(numout,*) ' '
-        END IF
+         IF(lwp) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) 'tra_qsr_init : penetration of the surface solar radiation'
+            WRITE(numout,*) '~~~~~~~~~~~~'
+            WRITE(numout,*) '    Namelist namqsr : set the parameter of penetration'
+            WRITE(numout,*) '        fraction associated with xsi     rabs        = ',rabs
+            WRITE(numout,*) '        first depth of extinction        xsi1        = ',xsi1
+            WRITE(numout,*) '        second depth of extinction       xsi2        = ',xsi2
+            IF( lk_qsr_sms ) THEN
+               WRITE(numout,*) '     Biological fluxes for light(Y/N) ln_qsr_sms  = ',ln_qsr_sms
+            ENDIF
+         ENDIF
       ELSE
-        IF ( lwp ) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'tra_qsr_init : NO solar flux penetration'
-         WRITE(numout,*) '~~~~~~~~~~~~'
-        END IF
+         IF(lwp) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) 'tra_qsr_init : NO solar flux penetration'
+            WRITE(numout,*) '~~~~~~~~~~~~'
+         ENDIF
       ENDIF
 
@@ -306,11 +252,12 @@
 
 
-      ! Initialization
-      ! --------------
-      IF( .NOT. lk_sco ) THEN
-         ! z-coordinate with or without partial step : same before last ocean w-level everywhere
+      ! Initialization of gdsr
+      ! ----------------------
+      IF( ln_zco .OR. ln_zps ) THEN
+
+         ! z-coordinate with or without partial step : same w-level everywhere inside the ocean
          gdsr(:) = 0.e0
          DO jk = 1, jpk
-            zdp1 = -fsdepw(1,1,jk)
+            zdp1 = -gdepw_0(jk)
             gdsr(jk) = ro0cpr * (  rabs  * EXP( zdp1/xsi1 ) + (1.-rabs) * EXP( zdp1/xsi2 )  )
             IF ( gdsr(jk) <= 1.e-10 ) EXIT
@@ -321,5 +268,4 @@
                gdsr(jk) = 0.e0
                nksr = jk
-               !!bug Edmee chg res   nksr = jk - 1
                indic = 1
             ENDIF
@@ -337,13 +283,28 @@
          IF( lk_qsr_sms .AND. ln_qsr_sms ) THEN
             DO jk = 1, jpkm1
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     etot3(ji,jj,jk) = qsr(ji,jj) * gdsr(jk) * tmask(ji,jj,jk) / ro0cpr
-                  END DO
+               zcst = gdsr(jk) / ro0cpr
+               etot3(:,:,jk) = qsr(:,:) * zcst * tmask(:,:,jk) 
+            END DO
+         ENDIF
+
+      ENDIF
+
+      ! Initialisation of etot3 (s-coordinate)
+      ! -----------------------
+      IF( ln_sco ) THEN
+         etot3(:,:,jpk) = 0.e0
+         DO jk = 1, jpkm1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zdp1 = -fsdepw(ji,jj,jk  )
+                  zdp2 = -fsdepw(ji,jj,jk+1)
+                  zc0  = ro0cpr / fse3t(ji,jj,jk)
+                  zc1  =   (  rabs * EXP(zdp1/xsi1) + (1.-rabs) * EXP(zdp1/xsi2)  )
+                  zc2  = - (  rabs * EXP(zdp2/xsi1) + (1.-rabs) * EXP(zdp2/xsi2)  )
+                  etot3(ji,jj,jk)  = zc0 * (  zc1 * tmask(ji,jj,jk) + zc2 * tmask(ji,jj,jk+1)  )
                END DO
             END DO
-         ENDIF
-
-      ENDIF
+         END DO 
+       ENDIF
 
    END SUBROUTINE tra_qsr_init
