Index: branches/2011/dev_r2782_NOCS_Griffies/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
===================================================================
--- branches/2011/dev_r2782_NOCS_Griffies/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 2840)
+++ branches/2011/dev_r2782_NOCS_Griffies/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 2969)
@@ -402,8 +402,8 @@
       !!----------------------------------------------------------------------
       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
-      USE oce     , ONLY:   zalbet  => ua   ! use ua as workspace
+      USE oce     , ONLY:   zalbet  => ua       ! use ua as workspace
       USE wrk_nemo, ONLY:   z1_mlbw => wrk_2d_1
       !!
-      INTEGER, INTENT( in ) ::   kt         ! ocean time-step index
+      INTEGER, INTENT( in ) ::   kt             ! ocean time-step index
       !!
       INTEGER  ::   ji, jj, jk, jl, ip, jp, kp  ! dummy loop indices
@@ -425,12 +425,12 @@
       !--------------------------------!
       !
-      CALL eos_alpbet( tsb, zalbet, zbeta0 )         !==  before local thermal/haline expension ratio at T-points  ==!
-      !
-      DO jl = 0, 1                           !==  unmasked before density i- j-, k-gradients  ==!
+      CALL eos_alpbet( tsb, zalbet, zbeta0 )  !==  before local thermal/haline expension ratio at T-points  ==!
+      !
+      DO jl = 0, 1                            !==  unmasked before density i- j-, k-gradients  ==!
          !
-         ip = jl   ;   jp = jl        ! guaranteed nonzero gradients ( absolute value larger than repsln)
-         DO jk = 1, jpkm1                          ! done each pair of triad
-            DO jj = 1, jpjm1                       ! NB: not masked ==>  a minimum value is set
-               DO ji = 1, fs_jpim1   ! vector opt.
+         ip = jl   ;   jp = jl                ! guaranteed nonzero gradients ( absolute value larger than repsln)
+         DO jk = 1, jpkm1                     ! done each pair of triad
+            DO jj = 1, jpjm1                  ! NB: not masked ==>  a minimum value is set
+               DO ji = 1, fs_jpim1            ! vector opt.
                   zdit = ( tsb(ji+1,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) )    ! i-gradient of T & S at u-point
                   zdis = ( tsb(ji+1,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) )
@@ -445,8 +445,8 @@
          END DO
          !
-         IF( ln_zps.and.ln_grad_zps ) THEN                       ! partial steps: correction of i- & j-grad on bottom
+         IF( ln_zps.and.l_grad_zps ) THEN     ! partial steps: correction of i- & j-grad on bottom
 # if defined key_vectopt_loop
             DO jj = 1, 1
-               DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
+               DO ji = 1, jpij-jpi            ! vector opt. (forced unrolling)
 # else
             DO jj = 1, jpjm1
@@ -466,9 +466,9 @@
       END DO
 
-      DO kp = 0, 1                           !==  unmasked before density i- j-, k-gradients  ==!
-         DO jk = 1, jpkm1                          ! done each pair of triad
-            DO jj = 1, jpj                       ! NB: not masked ==>  a minimum value is set
-               DO ji = 1, jpi  ! vector opt.
-                  IF( jk+kp > 1 ) THEN                                       ! k-gradient of T & S a jk+kp
+      DO kp = 0, 1                            !==  unmasked before density i- j-, k-gradients  ==!
+         DO jk = 1, jpkm1                     ! done each pair of triad
+            DO jj = 1, jpj                    ! NB: not masked ==>  a minimum value is set
+               DO ji = 1, jpi                 ! vector opt.
+                  IF( jk+kp > 1 ) THEN        ! k-gradient of T & S a jk+kp
                      zdkt = ( tsb(ji,jj,jk+kp-1,jp_tem) - tsb(ji,jj,jk+kp,jp_tem) )
                      zdks = ( tsb(ji,jj,jk+kp-1,jp_sal) - tsb(ji,jj,jk+kp,jp_sal) )
@@ -484,5 +484,5 @@
       END DO
       !
-      DO jj = 1, jpj                         !==  Reciprocal depth of the w-point below ML base  ==!
+      DO jj = 1, jpj                          !==  Reciprocal depth of the w-point below ML base  ==!
          DO ji = 1, jpi
             jk = MIN( nmln(ji,jj), mbkt(ji,jj) ) + 1     ! MIN in case ML depth is the ocean depth
@@ -491,10 +491,10 @@
       END DO
       !
-      !                                      !==  intialisations to zero  ==!
-      !
-      wslp2  (:,:,:)     = 0._wp                                        ! wslp2 will be cumulated 3D field set to zero
+      !                                       !==  intialisations to zero  ==!
+      !
+      wslp2  (:,:,:)     = 0._wp              ! wslp2 will be cumulated 3D field set to zero
       triadi_g(:,:,1,:,:) = 0._wp   ;   triadi_g(:,:,jpk,:,:) = 0._wp   ! set surface and bottom slope to zero
       triadj_g(:,:,1,:,:) = 0._wp   ;   triadj_g(:,:,jpk,:,:) = 0._wp
-!!gm _iso set to zero missing
+      !!gm _iso set to zero missing
       triadi  (:,:,1,:,:) = 0._wp   ;   triadj  (:,:,jpk,:,:) = 0._wp   ! set surface and bottom slope to zero
       triadj  (:,:,1,:,:) = 0._wp   ;   triadj  (:,:,jpk,:,:) = 0._wp
@@ -504,6 +504,6 @@
       !-------------------------------------!
       !
-      DO jl = 0, 1               ! calculate slope of the 4 triads immediately ONE level below mixed-layer base
-         DO kp = 0, 1            ! with only the slope-max limit   and   MASKED
+      DO jl = 0, 1                            ! calculate slope of the 4 triads immediately ONE level below mixed-layer base
+         DO kp = 0, 1                         ! with only the slope-max limit   and   MASKED
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1
@@ -527,10 +527,10 @@
       !-------------------------------------!
       !
-      DO kp = 0, 1                        ! k-index of triads
+      DO kp = 0, 1                            ! k-index of triads
          DO jl = 0, 1
-            ip = jl   ;   jp = jl         ! i- and j-indices of triads (i-k and j-k planes)
+            ip = jl   ;   jp = jl             ! i- and j-indices of triads (i-k and j-k planes)
             DO jk = 1, jpkm1
                DO jj = 1, jpjm1
-                  DO ji = 1, fs_jpim1   ! vector opt.
+                  DO ji = 1, fs_jpim1         ! vector opt.
                      !
                      ! Calculate slope relative to geopotentials used for GM skew fluxes
@@ -575,5 +575,5 @@
                      zti_lim  = ( zti_g_lim + zti_coord ) * umask(ji,jj,jk+kp)    ! remove coordinate slope => relative to coordinate surfaces
                      ztj_lim  = ( ztj_g_lim + ztj_coord ) * vmask(ji,jj,jk+kp)
-                     zti_lim2 = zti_lim * zti_lim      ! square of limited slopes            ! masked <<==
+                     zti_lim2 = zti_lim * zti_lim      ! square of limited slopes                          ! masked <<==
                      ztj_lim2 = ztj_lim * ztj_lim
                      !
@@ -596,6 +596,6 @@
                      zbtj = e1t(ji,jj+jp) * e2t(ji,jj+jp) * fse3w(ji,jj+jp,jk+kp)
                      !
-!!gm this may inhibit vectorization on Vect Computers, and even on scalar computers....  ==> to be checked
-! agn may need to change to using ztj_g_lim**2, as wslp2 just used for eddy growth rate, needs *real* slope
+                     !!gm this may inhibit vectorization on Vect Computers, and even on scalar computers....  ==> to be checked
+                     ! agn may need to change to using ztj_g_lim**2, as wslp2 just used for eddy growth rate, needs *real* slope
                      wslp2 (ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_lim2      ! masked
                      wslp2 (ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ztj_lim2
@@ -637,6 +637,6 @@
       REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   p_dzr          ! z-gradient of density      (T-point)
       !!
-      INTEGER  ::   ji , jj , jk         ! dummy loop indices
-      INTEGER  ::   iku, ikv, ik, ikm1   ! local integers
+      INTEGER  ::   ji , jj , jk                   ! dummy loop indices
+      INTEGER  ::   iku, ikv, ik, ikm1             ! local integers
       REAL(wp) ::   zeps, zm1_g, zm1_2g            ! local scalars
       REAL(wp) ::   zci, zfi, zau, zbu, zai, zbi   !   -      -
@@ -654,9 +654,9 @@
       wslpjml(1,:) = 0._wp      ;      wslpjml(jpi,:) = 0._wp
       !
-      !                          !==   surface mixed layer mask   !
-      DO jk = 1, jpk                      ! =1 inside the mixed layer, =0 otherwise
+      !                                            !==   surface mixed layer mask   !
+      DO jk = 1, jpk                               ! =1 inside the mixed layer, =0 otherwise
 # if defined key_vectopt_loop
          DO jj = 1, 1
-            DO ji = 1, jpij   ! vector opt. (forced unrolling)
+            DO ji = 1, jpij                        ! vector opt. (forced unrolling)
 # else
          DO jj = 1, jpj
@@ -689,29 +689,29 @@
          DO ji = 2, jpim1
 # endif
-            !                    !==   Slope at u- & v-points just below the Mixed Layer   ==!
+            !                        !==   Slope at u- & v-points just below the Mixed Layer   ==!
             !
-            !                          !- vertical density gradient for u- and v-slopes (from dzr at T-point)
+            !                        !- vertical density gradient for u- and v-slopes (from dzr at T-point)
             iku = MIN(  MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1  )   ! ML (MAX of T-pts, bound by jpkm1)
             ikv = MIN(  MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1  )   !
             zbu = 0.5_wp * ( p_dzr(ji,jj,iku) + p_dzr(ji+1,jj  ,iku) )
             zbv = 0.5_wp * ( p_dzr(ji,jj,ikv) + p_dzr(ji  ,jj+1,ikv) )
-            !                          !- horizontal density gradient at u- & v-points
+            !                        !- horizontal density gradient at u- & v-points
             zau = p_gru(ji,jj,iku) / e1u(ji,jj)
             zav = p_grv(ji,jj,ikv) / e2v(ji,jj)
-            !                          !- bound the slopes: abs(zw.)<= 1/100 and zb..<0
-            !                                kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
+            !                        !- bound the slopes: abs(zw.)<= 1/100 and zb..<0
+            !                           kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
             zbu = MIN(  zbu , -100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,iku)* ABS( zau )  )
             zbv = MIN(  zbv , -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,ikv)* ABS( zav )  )
-            !                          !- Slope at u- & v-points (uslpml, vslpml)
+            !                        !- Slope at u- & v-points (uslpml, vslpml)
             uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,iku)
             vslpml(ji,jj) = zav / ( zbv - zeps ) * vmask(ji,jj,ikv)
             !
-            !                    !==   i- & j-slopes at w-points just below the Mixed Layer   ==!
+            !                        !==   i- & j-slopes at w-points just below the Mixed Layer   ==!
             !
             ik   = MIN( nmln(ji,jj) + 1, jpk )
             ikm1 = MAX( 1, ik-1 )
-            !                          !- vertical density gradient for w-slope (from N^2)
+            !                        !- vertical density gradient for w-slope (from N^2)
             zbw = zm1_2g * pn2 (ji,jj,ik) * ( prd (ji,jj,ik) + prd (ji,jj,ikm1) + 2. )
-            !                          !- horizontal density i- & j-gradient at w-points
+            !                        !- horizontal density i- & j-gradient at w-points
             zci = MAX(   umask(ji-1,jj,ik  ) + umask(ji,jj,ik  )           &
                &       + umask(ji-1,jj,ikm1) + umask(ji,jj,ikm1) , zeps  ) * e1t(ji,jj)
@@ -722,14 +722,14 @@
             zaj =    (   p_grv(ji,jj-1,ik  ) + p_grv(ji,jj,ik  )           &
                &       + p_grv(ji,jj-1,ikm1) + p_grv(ji,jj,ikm1)  ) / zcj  * tmask(ji,jj,ik)
-            !                          !- bound the slopes: abs(zw.)<= 1/100 and zb..<0.
-            !                             kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
+            !                        !- bound the slopes: abs(zw.)<= 1/100 and zb..<0.
+            !                           kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
             zbi = MIN(  zbw , -100._wp* ABS( zai ) , -7.e+3_wp/fse3w(ji,jj,ik)* ABS( zai )  )
             zbj = MIN(  zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,ik)* ABS( zaj )  )
-            !                          !- i- & j-slope at w-points (wslpiml, wslpjml)
+            !                        !- i- & j-slope at w-points (wslpiml, wslpjml)
             wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask (ji,jj,ik)
             wslpjml(ji,jj) = zaj / ( zbj - zeps ) * tmask (ji,jj,ik)
          END DO
       END DO
-!!gm this lbc_lnk should be useless....
+      !!gm this lbc_lnk should be useless....
       CALL lbc_lnk( uslpml , 'U', -1. )   ;   CALL lbc_lnk( vslpml , 'V', -1. )   ! lateral boundary cond. (sign change)
       CALL lbc_lnk( wslpiml, 'W', -1. )   ;   CALL lbc_lnk( wslpjml, 'W', -1. )   ! lateral boundary conditions
@@ -765,7 +765,4 @@
          IF( ln_dynldf_iso )   CALL ctl_stop( 'ldf_slp_init: Griffies operator on momentum not supported' )
          !
-         ! IF( ( ln_traldf_hor .OR. ln_dynldf_hor ) .AND. ln_sco )   &
-         !    CALL ctl_stop( 'ldf_slp_init: horizontal Griffies operator in s-coordinate not supported' )
-         !
       ELSE                             ! Madec operator : slopes at u-, v-, and w-points
          ALLOCATE( uslp(jpi,jpj,jpk) , vslp(jpi,jpj,jpk) , wslpi(jpi,jpj,jpk) , wslpj(jpi,jpj,jpk) ,                &
@@ -780,5 +777,5 @@
          wslpj(:,:,:) = 0._wp   ;   wslpjml(:,:) = 0._wp
 
-!!gm I no longer understand this.....
+         !!gm I no longer understand this.....
          IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. (lk_vvl .AND. ln_rstart) ) THEN
             IF(lwp)   WRITE(numout,*) '          Horizontal mixing in s-coordinate: slope = slope of s-surfaces'
@@ -813,5 +810,5 @@
    LOGICAL, PUBLIC, PARAMETER ::   lk_ldfslp = .FALSE.    !: slopes flag
 CONTAINS
-   SUBROUTINE ldf_slp( kt, prd, pn2 )        ! Dummy routine
+   SUBROUTINE ldf_slp( kt, prd, pn2 )   ! Dummy routine
       INTEGER, INTENT(in) :: kt
       REAL, DIMENSION(:,:,:), INTENT(in) :: prd, pn2
@@ -822,5 +819,5 @@
       WRITE(*,*) 'ldf_slp_grif: You should not have seen this print! error?', kt
    END SUBROUTINE ldf_slp_grif
-   SUBROUTINE ldf_slp_init       ! Dummy routine
+   SUBROUTINE ldf_slp_init              ! Dummy routine
    END SUBROUTINE ldf_slp_init
 #endif
