Index: trunk/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
===================================================================
--- trunk/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 7698)
+++ trunk/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 7753)
@@ -114,31 +114,11 @@
       zviold = v_i
       zvsold = v_s
-!$OMP PARALLEL
-!$OMP DO schedule(static) private(jj,ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zsmvold(ji,jj) = 0._wp
-         END DO
-      END DO
+      zsmvold(:,:) = SUM( smv_i(:,:,:), dim=3 )
+      zeiold (:,:) = et_i
+      zesold (:,:) = et_s 
+
+      !--- Thickness correction init. --- !
+      zatold(:,:) = at_i
       DO jl = 1, jpl
-!$OMP DO schedule(static) private(jj,ji)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zsmvold(ji,jj) = zsmvold(ji,jj) + smv_i(ji,jj,jl)
-            END DO
-         END DO
-      END DO
-!$OMP DO schedule(static) private(jj,ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zeiold (ji,jj) = et_i(ji,jj)
-            zesold (ji,jj) = et_s(ji,jj)
-
-            !--- Thickness correction init. --- !
-            zatold (ji,jj) = at_i(ji,jj)
-         END DO
-      END DO
-      DO jl = 1, jpl
-!$OMP DO schedule(static) private(jj,ji,rswitch)
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -150,15 +130,6 @@
       END DO
       ! --- Record max of the surrounding ice thicknesses for correction in case advection creates ice too thick --- !
+      zhimax(:,:,:) = ht_i(:,:,:) + ht_s(:,:,:)
       DO jl = 1, jpl
-!$OMP DO schedule(static) private(jj,ji)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zhimax(ji,jj,jl) = ht_i(ji,jj,jl) + ht_s(ji,jj,jl)
-            END DO
-         END DO
-      END DO
-!$OMP END PARALLEL
-      DO jl = 1, jpl
-!$OMP PARALLEL DO schedule(static) private(jj,ji)
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
@@ -202,16 +173,9 @@
          zdt = rdt_ice / REAL(initad)
          
-!$OMP PARALLEL
          ! transport
-!$OMP DO schedule(static) private(jj,ji)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zudy(ji,jj) = u_ice(ji,jj) * e2u(ji,jj)
-               zvdx(ji,jj) = v_ice(ji,jj) * e1v(ji,jj)
-            END DO
-         END DO
+         zudy(:,:) = u_ice(:,:) * e2u(:,:)
+         zvdx(:,:) = v_ice(:,:) * e1v(:,:)
          
          ! define velocity for advection: u*grad(H)
-!$OMP DO schedule(static) private(jj,ji)
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1
@@ -227,5 +191,4 @@
             END DO
          END DO
-!$OMP END PARALLEL
          
          ! advection
@@ -245,20 +208,8 @@
          END DO
          !
-!$OMP PARALLEL
-!$OMP DO schedule(static) private(jj,ji)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               at_i(ji,jj) = a_i(ji,jj,1)      ! total ice fraction
-            END DO
-         END DO
+         at_i(:,:) = a_i(:,:,1)      ! total ice fraction
          DO jl = 2, jpl
-!$OMP DO schedule(static) private(jj,ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl)
-               END DO
-            END DO
-         END DO
-!$OMP END PARALLEL
+            at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
+         END DO
          !
          CALL wrk_dealloc( jpi,jpj, zudy, zvdx, zcu_box, zcv_box )
@@ -279,38 +230,21 @@
          ENDIF
          
-!$OMP PARALLEL
-!$OMP DO schedule(static) private(jj,ji)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zarea(ji,jj) = e1e2t(ji,jj)
-         
-               !-------------------------
-               ! transported fields                                        
-               !-------------------------
-               z0opw(ji,jj,1) = ato_i(ji,jj) * e1e2t(ji,jj)             ! Open water area 
-            END DO
-         END DO
+         zarea(:,:) = e1e2t(:,:)
+         
+         !-------------------------
+         ! transported fields                                        
+         !-------------------------
+         z0opw(:,:,1) = ato_i(:,:) * e1e2t(:,:)             ! Open water area 
          DO jl = 1, jpl
-!$OMP DO schedule(static) private(jj,ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  z0snw (ji,jj,jl)  = v_s  (ji,jj,  jl) * e1e2t(ji,jj)  ! Snow volume
-                  z0ice(ji,jj,jl)   = v_i  (ji,jj,  jl) * e1e2t(ji,jj)  ! Ice  volume
-                  z0ai  (ji,jj,jl)  = a_i  (ji,jj,  jl) * e1e2t(ji,jj)  ! Ice area
-                  z0smi (ji,jj,jl)  = smv_i(ji,jj,  jl) * e1e2t(ji,jj)  ! Salt content
-                  z0oi (ji,jj,jl)   = oa_i (ji,jj,  jl) * e1e2t(ji,jj)  ! Age content
-                  z0es (ji,jj,jl)   = e_s  (ji,jj,1,jl) * e1e2t(ji,jj)  ! Snow heat content
-               END DO
-            END DO
+            z0snw (:,:,jl)  = v_s  (:,:,  jl) * e1e2t(:,:)  ! Snow volume
+            z0ice(:,:,jl)   = v_i  (:,:,  jl) * e1e2t(:,:)  ! Ice  volume
+            z0ai  (:,:,jl)  = a_i  (:,:,  jl) * e1e2t(:,:)  ! Ice area
+            z0smi (:,:,jl)  = smv_i(:,:,  jl) * e1e2t(:,:)  ! Salt content
+            z0oi (:,:,jl)   = oa_i (:,:,  jl) * e1e2t(:,:)  ! Age content
+            z0es (:,:,jl)   = e_s  (:,:,1,jl) * e1e2t(:,:)  ! Snow heat content
             DO jk = 1, nlay_i
-!$OMP DO schedule(static) private(jj,ji)
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     z0ei  (ji,jj,jk,jl) = e_i  (ji,jj,jk,jl) * e1e2t(ji,jj) ! Ice  heat content
-                  END DO
-               END DO
-            END DO
-         END DO
-!$OMP END PARALLEL
+               z0ei  (:,:,jk,jl) = e_i  (:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content
+            END DO
+         END DO
 
 
@@ -402,48 +336,21 @@
          ! Recover the properties from their contents
          !-------------------------------------------
-!$OMP PARALLEL
-!$OMP DO schedule(static) private(jj,ji)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ato_i(ji,jj) = z0opw(ji,jj,1) * r1_e1e2t(ji,jj)
-            END DO
-         END DO
+         ato_i(:,:) = z0opw(:,:,1) * r1_e1e2t(:,:)
          DO jl = 1, jpl
-!$OMP DO schedule(static) private(jj,ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  v_i  (ji,jj,  jl) = z0ice(ji,jj,jl) * r1_e1e2t(ji,jj)
-                  v_s  (ji,jj,  jl) = z0snw(ji,jj,jl) * r1_e1e2t(ji,jj)
-                  smv_i(ji,jj,  jl) = z0smi(ji,jj,jl) * r1_e1e2t(ji,jj)
-                  oa_i (ji,jj,  jl) = z0oi (ji,jj,jl) * r1_e1e2t(ji,jj)
-                  a_i  (ji,jj,  jl) = z0ai (ji,jj,jl) * r1_e1e2t(ji,jj)
-                  e_s  (ji,jj,1,jl) = z0es (ji,jj,jl) * r1_e1e2t(ji,jj)
-               END DO
-            END DO
+            v_i  (:,:,  jl) = z0ice(:,:,jl) * r1_e1e2t(:,:)
+            v_s  (:,:,  jl) = z0snw(:,:,jl) * r1_e1e2t(:,:)
+            smv_i(:,:,  jl) = z0smi(:,:,jl) * r1_e1e2t(:,:)
+            oa_i (:,:,  jl) = z0oi (:,:,jl) * r1_e1e2t(:,:)
+            a_i  (:,:,  jl) = z0ai (:,:,jl) * r1_e1e2t(:,:)
+            e_s  (:,:,1,jl) = z0es (:,:,jl) * r1_e1e2t(:,:)
             DO jk = 1, nlay_i
-!$OMP DO schedule(static) private(jj,ji)
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     e_i(ji,jj,jk,jl) = z0ei(ji,jj,jk,jl) * r1_e1e2t(ji,jj)
-                  END DO
-               END DO
-            END DO
-         END DO
-
-!$OMP DO schedule(static) private(jj,ji)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               at_i(ji,jj) = a_i(ji,jj,1)      ! total ice fraction
-            END DO
-         END DO
+               e_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:)
+            END DO
+         END DO
+
+         at_i(:,:) = a_i(:,:,1)      ! total ice fraction
          DO jl = 2, jpl
-!$OMP DO schedule(static) private(jj,ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl)
-               END DO
-            END DO
-         END DO
-!$OMP END PARALLEL
+            at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
+         END DO
          
          CALL wrk_dealloc( jpi,jpj,            zarea )
@@ -462,7 +369,5 @@
          !     mask eddy diffusivity coefficient at ocean U- and V-points
          jm=1
-!$OMP PARALLEL
          DO jl = 1, jpl
-!$OMP DO schedule(static) private(jj,ji)
             DO jj = 1, jpjm1                 ! NB: has not to be defined on jpj line and jpi row
                DO ji = 1 , fs_jpim1
@@ -474,57 +379,15 @@
             END DO
 
-!$OMP DO schedule(static) private(jj,ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zhdfptab(ji,jj,jm)= a_i  (ji,jj,  jl)
-               END DO
-            END DO
-            jm = jm + 1
-!$OMP DO schedule(static) private(jj,ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zhdfptab(ji,jj,jm)= v_i  (ji,jj,  jl)
-               END DO
-            END DO
-            jm = jm + 1
-!$OMP DO schedule(static) private(jj,ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zhdfptab(ji,jj,jm)= v_s  (ji,jj,  jl)
-               END DO
-            END DO
-            jm = jm + 1
-!$OMP DO schedule(static) private(jj,ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zhdfptab(ji,jj,jm)= smv_i(ji,jj,  jl)
-               END DO
-            END DO
-            jm = jm + 1
-!$OMP DO schedule(static) private(jj,ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zhdfptab(ji,jj,jm)= oa_i (ji,jj,  jl)
-               END DO
-            END DO
-            jm = jm + 1
-!$OMP DO schedule(static) private(jj,ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zhdfptab(ji,jj,jm)= e_s  (ji,jj,1,jl)
-               END DO
-            END DO
-            jm = jm + 1
+            zhdfptab(:,:,jm)= a_i  (:,:,  jl); jm = jm + 1
+            zhdfptab(:,:,jm)= v_i  (:,:,  jl); jm = jm + 1
+            zhdfptab(:,:,jm)= v_s  (:,:,  jl); jm = jm + 1
+            zhdfptab(:,:,jm)= smv_i(:,:,  jl); jm = jm + 1
+            zhdfptab(:,:,jm)= oa_i (:,:,  jl); jm = jm + 1
+            zhdfptab(:,:,jm)= e_s  (:,:,1,jl); jm = jm + 1
             ! Sample of adding more variables to apply lim_hdf (ihdf_vars must be increased)
             !   zhdfptab(:,:,jm) = variable_1 (:,:,1,jl); jm = jm + 1  
             !   zhdfptab(:,:,jm) = variable_2 (:,:,1,jl); jm = jm + 1 
             DO jk = 1, nlay_i
-!$OMP DO schedule(static) private(jj,ji)
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zhdfptab(ji,jj,jm)=e_i(ji,jj,jk,jl)
-                  END DO
-               END DO
-               jm= jm+1
+              zhdfptab(:,:,jm)=e_i(:,:,jk,jl); jm= jm+1
             END DO
          END DO
@@ -532,5 +395,4 @@
          ! --- Prepare diffusion for open water area --- !
          !     mask eddy diffusivity coefficient at ocean U- and V-points
-!$OMP DO schedule(static) private(jj,ji)
          DO jj = 1, jpjm1                    ! NB: has not to be defined on jpj line and jpi row
             DO ji = 1 , fs_jpim1
@@ -542,11 +404,5 @@
          END DO
          !
-!$OMP DO schedule(static) private(jj,ji)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zhdfptab(ji,jj,jm)= ato_i  (ji,jj);
-            END DO
-         END DO
-!$OMP END PARALLEL
+         zhdfptab(:,:,jm)= ato_i  (:,:);
 
          ! --- Apply diffusion --- !
@@ -555,74 +411,23 @@
          ! --- Recover properties --- !
          jm=1
-!$OMP PARALLEL
          DO jl = 1, jpl
-!$OMP DO schedule(static) private(jj,ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  a_i  (ji,jj,  jl)=zhdfptab(ji,jj,jm)
-               END DO
-            END DO
-            jm = jm + 1
-!$OMP DO schedule(static) private(jj,ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  v_i  (ji,jj,  jl)=zhdfptab(ji,jj,jm)
-               END DO
-            END DO
-            jm = jm + 1
-!$OMP DO schedule(static) private(jj,ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  v_s  (ji,jj,  jl)=zhdfptab(ji,jj,jm)
-               END DO
-            END DO
-            jm = jm + 1
-!$OMP DO schedule(static) private(jj,ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  smv_i(ji,jj,  jl)=zhdfptab(ji,jj,jm)
-               END DO
-            END DO
-            jm = jm + 1
-!$OMP DO schedule(static) private(jj,ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  oa_i (ji,jj,  jl)=zhdfptab(ji,jj,jm)
-               END DO
-            END DO
-            jm = jm + 1
-!$OMP DO schedule(static) private(jj,ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  e_s  (ji,jj,1,jl)=zhdfptab(ji,jj,jm)
-               END DO
-            END DO
-            jm = jm + 1
-
+            a_i  (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1
+            v_i  (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1
+            v_s  (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1
+            smv_i(:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1
+            oa_i (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1
+            e_s  (:,:,1,jl) = zhdfptab(:,:,jm); jm = jm + 1
             ! Sample of adding more variables to apply lim_hdf
             !   variable_1  (:,:,1,jl) = zhdfptab(:,:, jm  ) ; jm + 1 
             !   variable_2  (:,:,1,jl) = zhdfptab(:,:, jm  ) ; jm + 1
             DO jk = 1, nlay_i
-!$OMP DO schedule(static) private(jj,ji)
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     e_i(ji,jj,jk,jl) = zhdfptab(ji,jj,jm)
-                  END DO
-               END DO
-               jm = jm + 1
-            END DO
-         END DO
-!$OMP DO schedule(static) private(jj,ji)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ato_i  (ji,jj) = zhdfptab(ji,jj,jm)
-            END DO
-         END DO
-!$OMP END PARALLEL
+               e_i(:,:,jk,jl) = zhdfptab(:,:,jm);jm= jm + 1
+            END DO
+         END DO
+         ato_i  (:,:) = zhdfptab(:,:,jm)
               
       ENDIF
 
       ! --- diags ---
-!$OMP PARALLEL DO schedule(static) private(jj,ji)
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -641,7 +446,5 @@
            
          !--- Thickness correction in case too high --- !
-!$OMP PARALLEL
          DO jl = 1, jpl
-!$OMP DO schedule(static) private(jj,ji,rswitch,zdv)
             DO jj = 1, jpj
                DO ji = 1, jpi
@@ -678,5 +481,4 @@
          
          ! Force the upper limit of ht_i to always be < hi_max (99 m).
-!$OMP DO schedule(static) private(jj,ji,rswitch)
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -686,5 +488,4 @@
             END DO
          END DO
-!$OMP END PARALLEL
 
       ENDIF
@@ -694,24 +495,7 @@
       !------------------------------------------------------------
       !
-!$OMP PARALLEL
-!$OMP DO schedule(static) private(jj,ji)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               at_i(ji,jj) = 0._wp
-            END DO
-         END DO
-         DO jl = 1, jpl
-!$OMP DO schedule(static) private(jj,ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl)
-               END DO
-            END DO
-         END DO
-!$OMP END PARALLEL
-
+      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
       IF ( nn_limdyn == 1 .OR. ( ( nn_monocat == 2 ) .AND. ( jpl == 1 ) ) ) THEN ! simple conservative piling, comparable with LIM2
          DO jl = 1, jpl
-!$OMP PARALLEL DO schedule(static) private(jj,ji,rswitch,zda)
             DO jj = 1, jpj
                DO ji = 1, jpi
@@ -726,32 +510,10 @@
       
       ! --- agglomerate variables -----------------
-!$OMP PARALLEL
-!$OMP DO schedule(static) private(jj,ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            vt_i(ji,jj) = 0._wp
-            vt_s(ji,jj) = 0._wp
-            at_i(ji,jj) = 0._wp
-         END DO
-      END DO
-      DO jl = 1, jpl
-!$OMP DO schedule(static) private(jj,ji)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl)
-               vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl)
-               at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl)
-            END DO
-         END DO
-      END DO
+      vt_i(:,:) = SUM( v_i(:,:,:), dim=3 )
+      vt_s(:,:) = SUM( v_s(:,:,:), dim=3 )
+      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
       
       ! --- open water = 1 if at_i=0 --------------------------------
-!$OMP DO schedule(static) private(jj,ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            IF( at_i(ji,jj) == 0._wp ) ato_i(ji,jj) = 1._wp 
-         END DO
-      END DO
-!$OMP END PARALLEL
+      WHERE( at_i == 0._wp ) ato_i = 1._wp 
       
       ! conservation test
