Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90	(revision 7646)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90	(revision 7698)
@@ -115,7 +115,16 @@
          
          !  Computation of ice albedo (free of snow)
-         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb(:,:,:) = ralb_im
-         ELSE WHERE                                              ;   zalb(:,:,:) = ralb_if
-         END  WHERE
+!$OMP PARALLEL DO schedule(static) private(jl,jj,ji)
+         DO jl = 1, ijpl
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  IF ( ph_snw(ji,jj,jl) == 0._wp .AND. pt_ice(ji,jj,jl) >= rt0_ice ) THEN
+                     zalb(ji,jj,jl) = ralb_im
+                  ELSE
+                     zalb(ji,jj,jl) = ralb_if
+                  END IF
+               END DO
+            END DO
+         END DO
       
          WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb
@@ -126,5 +135,6 @@
          ELSE WHERE                                       ;  zalb_it = 0.1    + 3.6    * ph_ice
          END WHERE
-     
+!$OMP PARALLEL
+!$OMP DO schedule(static) private(jl, jj, ji,zswitch,zalb_sf,zalb_sm,zalb_st)
          DO jl = 1, ijpl
             DO jj = 1, jpj
@@ -156,5 +166,13 @@
          END DO
 
-         pa_ice_os(:,:,:) = pa_ice_cs(:,:,:) + rcloud       ! Oberhuber correction for overcast sky
+!$OMP DO schedule(static) private(jl, jj, ji)     
+         DO jl = 1, ijpl
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  pa_ice_os(ji,jj,jl) = pa_ice_cs(ji,jj,jl) + rcloud       ! Oberhuber correction for overcast sky
+               END DO
+            END DO
+         END DO
+!$OMP END PARALLEL
 
       !------------------------------------------
@@ -193,4 +211,5 @@
          z1_c2 = 1. / 0.03
          !  Computation of the snow/ice albedo
+!$OMP PARALLEL DO schedule(static) private(jl, jj, ji,zswitch,zalb_sf,zalb_sm,zalb_st)     
          DO jl = 1, ijpl
             DO jj = 1, jpj
@@ -230,9 +249,15 @@
       !!
       REAL(wp) :: zcoef 
+      INTEGER  ::   ji, jj                                   ! dummy loop indices
       !!----------------------------------------------------------------------
       !
       zcoef = 0.05 / ( 1.1 * rmue**1.4 + 0.15 )   ! Parameterization of Briegled and Ramanathan, 1982
-      pa_oce_cs(:,:) = zcoef 
-      pa_oce_os(:,:) = 0.06                       ! Parameterization of Kondratyev, 1969 and Payne, 1972
+!$OMP PARALLEL DO schedule(static) private(jj, ji)     
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            pa_oce_cs(ji,jj) = zcoef 
+            pa_oce_os(ji,jj) = 0.06                       ! Parameterization of Kondratyev, 1969 and Payne, 1972
+         END DO
+      END DO
       !
    END SUBROUTINE albedo_oce
Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90	(revision 7646)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90	(revision 7698)
@@ -66,4 +66,5 @@
       !                                                             ! 'ij->e' = (i,j) components to east
       !                                                             ! 'ij->n' = (i,j) components to north
+      INTEGER  ::   ji, jj                                          ! dummy loop indices
       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   prot      
       !!----------------------------------------------------------------------
@@ -82,32 +83,128 @@
       CASE( 'en->i' )                  ! east-north to i-component
          SELECT CASE (cd_type)
-         CASE ('T')   ;   prot(:,:) = pxin(:,:) * gcost(:,:) + pyin(:,:) * gsint(:,:)
-         CASE ('U')   ;   prot(:,:) = pxin(:,:) * gcosu(:,:) + pyin(:,:) * gsinu(:,:)
-         CASE ('V')   ;   prot(:,:) = pxin(:,:) * gcosv(:,:) + pyin(:,:) * gsinv(:,:)
-         CASE ('F')   ;   prot(:,:) = pxin(:,:) * gcosf(:,:) + pyin(:,:) * gsinf(:,:)
+         CASE ('T')   
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  prot(ji,jj) = pxin(ji,jj) * gcost(ji,jj) + pyin(ji,jj) * gsint(ji,jj)
+               END DO
+            END DO
+         CASE ('U')
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  prot(ji,jj) = pxin(ji,jj) * gcosu(ji,jj) + pyin(ji,jj) * gsinu(ji,jj)
+               END DO
+            END DO
+         CASE ('V')
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  prot(ji,jj) = pxin(ji,jj) * gcosv(ji,jj) + pyin(ji,jj) * gsinv(ji,jj)
+               END DO
+            END DO
+         CASE ('F')
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  prot(ji,jj) = pxin(ji,jj) * gcosf(ji,jj) + pyin(ji,jj) * gsinf(ji,jj)
+               END DO
+            END DO
          CASE DEFAULT   ;   CALL ctl_stop( 'Only T, U, V and F grid points are coded' )
          END SELECT
       CASE ('en->j')                   ! east-north to j-component
          SELECT CASE (cd_type)
-         CASE ('T')   ;   prot(:,:) = pyin(:,:) * gcost(:,:) - pxin(:,:) * gsint(:,:)
-         CASE ('U')   ;   prot(:,:) = pyin(:,:) * gcosu(:,:) - pxin(:,:) * gsinu(:,:)
-         CASE ('V')   ;   prot(:,:) = pyin(:,:) * gcosv(:,:) - pxin(:,:) * gsinv(:,:)   
-         CASE ('F')   ;   prot(:,:) = pyin(:,:) * gcosf(:,:) - pxin(:,:) * gsinf(:,:)   
+         CASE ('T')
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  prot(ji,jj) = pyin(ji,jj) * gcost(ji,jj) - pxin(ji,jj) * gsint(ji,jj)
+               END DO
+            END DO
+         CASE ('U')
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  prot(ji,jj) = pyin(ji,jj) * gcosu(ji,jj) - pxin(ji,jj) * gsinu(ji,jj)
+               END DO
+            END DO
+         CASE ('V')   
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  prot(ji,jj) = pyin(ji,jj) * gcosv(ji,jj) - pxin(ji,jj) * gsinv(ji,jj)
+               END DO
+            END DO
+         CASE ('F')   
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  prot(ji,jj) = pyin(ji,jj) * gcosf(ji,jj) - pxin(ji,jj) * gsinf(ji,jj)
+               END DO
+            END DO
          CASE DEFAULT   ;   CALL ctl_stop( 'Only T, U, V and F grid points are coded' )
          END SELECT
       CASE ('ij->e')                   ! (i,j)-components to east
          SELECT CASE (cd_type)
-         CASE ('T')   ;   prot(:,:) = pxin(:,:) * gcost(:,:) - pyin(:,:) * gsint(:,:)
-         CASE ('U')   ;   prot(:,:) = pxin(:,:) * gcosu(:,:) - pyin(:,:) * gsinu(:,:)
-         CASE ('V')   ;   prot(:,:) = pxin(:,:) * gcosv(:,:) - pyin(:,:) * gsinv(:,:)
-         CASE ('F')   ;   prot(:,:) = pxin(:,:) * gcosf(:,:) - pyin(:,:) * gsinf(:,:)
+         CASE ('T')
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  prot(ji,jj) = pxin(ji,jj) * gcost(ji,jj) - pyin(ji,jj) * gsint(ji,jj)
+               END DO
+            END DO
+         CASE ('U')
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  prot(ji,jj) = pxin(ji,jj) * gcosu(ji,jj) - pyin(ji,jj) * gsinu(ji,jj)
+               END DO
+            END DO
+         CASE ('V')
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  prot(ji,jj) = pxin(ji,jj) * gcosv(ji,jj) - pyin(ji,jj) * gsinv(ji,jj)
+               END DO
+            END DO
+         CASE ('F')
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  prot(ji,jj) = pxin(ji,jj) * gcosf(ji,jj) - pyin(ji,jj) * gsinf(ji,jj)
+               END DO
+            END DO
          CASE DEFAULT   ;   CALL ctl_stop( 'Only T, U, V and F grid points are coded' )
          END SELECT
       CASE ('ij->n')                   ! (i,j)-components to north 
          SELECT CASE (cd_type)
-         CASE ('T')   ;   prot(:,:) = pyin(:,:) * gcost(:,:) + pxin(:,:) * gsint(:,:)
-         CASE ('U')   ;   prot(:,:) = pyin(:,:) * gcosu(:,:) + pxin(:,:) * gsinu(:,:)
-         CASE ('V')   ;   prot(:,:) = pyin(:,:) * gcosv(:,:) + pxin(:,:) * gsinv(:,:)
-         CASE ('F')   ;   prot(:,:) = pyin(:,:) * gcosf(:,:) + pxin(:,:) * gsinf(:,:)
+         CASE ('T')
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  prot(ji,jj) = pyin(ji,jj) * gcost(ji,jj) + pxin(ji,jj) * gsint(ji,jj)
+               END DO
+            END DO
+         CASE ('U')
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  prot(ji,jj) = pyin(ji,jj) * gcosu(ji,jj) + pxin(ji,jj) * gsinu(ji,jj)
+               END DO
+            END DO
+         CASE ('V')
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  prot(ji,jj) = pyin(ji,jj) * gcosv(ji,jj) + pxin(ji,jj) * gsinv(ji,jj)
+               END DO
+            END DO
+         CASE ('F')
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  prot(ji,jj) = pyin(ji,jj) * gcosf(ji,jj) + pxin(ji,jj) * gsinf(ji,jj)
+               END DO
+            END DO
          CASE DEFAULT   ;   CALL ctl_stop( 'Only T, U, V and F grid points are coded' )
          END SELECT
@@ -157,4 +254,7 @@
       ! (computation done on the north stereographic polar plane)
       !
+!$OMP PARALLEL
+!$OMP DO schedule(static) private(jj,ji,zlam,zphi,zxnpt,zynpt,znnpt,zxnpu,zynpu,znnpu,zxnpv,zynpv,znnpv,zxnpf) &
+!$OMP& private(zynpf,znnpf,zlan,zphh,zxvvt,zyvvt,znvvt,zxffu,zyffu,znffu,zxffv,zyffv,znffv,zxuuf,zyuuf,znuuf)
       DO jj = 2, jpjm1
          DO ji = fs_2, jpi   ! vector opt.
@@ -248,4 +348,5 @@
       ! =============== !
 
+!$OMP DO schedule(static) private(jj,ji)
       DO jj = 2, jpjm1
          DO ji = fs_2, jpi   ! vector opt.
@@ -268,4 +369,6 @@
          END DO
       END DO
+!$OMP END DO NOWAIT
+!$OMP END PARALLEL
 
       ! =========================== !
Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90	(revision 7646)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90	(revision 7698)
@@ -316,14 +316,30 @@
 #if defined key_cice
       IF( MOD( kt - 1, nn_fsbc ) == 0 )   THEN
-         qlw_ice(:,:,1)   = sf(jp_qlw )%fnow(:,:,1)
-         IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )
-         ELSE                ; qsr_ice(:,:,1) =          sf(jp_qsr)%fnow(:,:,1) 
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               qlw_ice(ji,jj,1)   = sf(jp_qlw)%fnow(ji,jj,1)
+            END DO
+         END DO
+         IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1)   = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )
+         ELSE                
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  qsr_ice(ji,jj,1)   = sf(jp_qsr)%fnow(ji,jj,1) 
+               END DO
+            END DO
          ENDIF 
-         tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1)
-         qatm_ice(:,:)    = sf(jp_humi)%fnow(:,:,1)
-         tprecip(:,:)     = sf(jp_prec)%fnow(:,:,1) * rn_pfac
-         sprecip(:,:)     = sf(jp_snow)%fnow(:,:,1) * rn_pfac
-         wndi_ice(:,:)    = sf(jp_wndi)%fnow(:,:,1)
-         wndj_ice(:,:)    = sf(jp_wndj)%fnow(:,:,1)
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               tatm_ice(ji,jj)    = sf(jp_tair)%fnow(ji,jj,1)
+               qatm_ice(ji,jj)    = sf(jp_humi)%fnow(ji,jj,1)
+               tprecip(ji,jj)     = sf(jp_prec)%fnow(ji,jj,1) * rn_pfac
+               sprecip(ji,jj)     = sf(jp_snow)%fnow(ji,jj,1) * rn_pfac
+               wndi_ice(ji,jj)    = sf(jp_wndi)%fnow(ji,jj,1)
+               wndj_ice(ji,jj)    = sf(jp_wndj)%fnow(ji,jj,1)
+            END DO
+         END DO
       ENDIF
 #endif
@@ -382,17 +398,23 @@
       !
 
-      ! local scalars ( place there for vector optimisation purposes)
-      zst(:,:) = pst(:,:) + rt0      ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K)
-
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+         ! local scalars ( place there for vector optimisation purposes)
+            zst(ji,jj) = pst(ji,jj) + rt0      ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K)
+
+            ! ... components ( U10m - U_oce ) at T-point (unmasked)
+!!gm    move zwnd_i (_j) set to zero  inside the key_cyclone ???
+            zwnd_i(ji,jj) = 0._wp
+            zwnd_j(ji,jj) = 0._wp
+         END DO
+      END DO
       ! ----------------------------------------------------------------------------- !
       !      0   Wind components and module at T-point relative to the moving ocean   !
       ! ----------------------------------------------------------------------------- !
 
-      ! ... components ( U10m - U_oce ) at T-point (unmasked)
-!!gm    move zwnd_i (_j) set to zero  inside the key_cyclone ???
-      zwnd_i(:,:) = 0._wp
-      zwnd_j(:,:) = 0._wp
 #if defined key_cyclone
       CALL wnd_cyc( kt, zwnd_i, zwnd_j )    ! add analytical tropical cyclone (Vincent et al. JGR 2012)
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vect. opt.
@@ -402,4 +424,5 @@
       END DO
 #endif
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vect. opt.
@@ -411,7 +434,12 @@
       CALL lbc_lnk( zwnd_j(:,:) , 'T', -1. )
       ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked)
-      wndm(:,:) = SQRT(  zwnd_i(:,:) * zwnd_i(:,:)   &
-         &             + zwnd_j(:,:) * zwnd_j(:,:)  ) * tmask(:,:,1)
-
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            wndm(ji,jj) = SQRT(  zwnd_i(ji,jj) * zwnd_i(ji,jj)   &
+               &             + zwnd_j(ji,jj) * zwnd_j(ji,jj)  ) * tmask(ji,jj,1)
+
+         END DO
+      END DO
       ! ----------------------------------------------------------------------------- !
       !      I   Radiative FLUXES                                                     !
@@ -421,8 +449,19 @@
       zztmp = 1. - albo
       IF( ln_dm2dc ) THEN   ;   qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1)
-      ELSE                  ;   qsr(:,:) = zztmp *          sf(jp_qsr)%fnow(:,:,1)   * tmask(:,:,1)
+      ELSE         
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               qsr(ji,jj) = zztmp *          sf(jp_qsr)%fnow(ji,jj,1)   * tmask(ji,jj,1)
+            END DO
+         END DO
       ENDIF
 
-      zqlw(:,:) = (  sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:)  ) * tmask(:,:,1)   ! Long  Wave
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zqlw(ji,jj) = (  sf(jp_qlw)%fnow(ji,jj,1) - Stef * zst(ji,jj)*zst(ji,jj)*zst(ji,jj)*zst(ji,jj)  ) * tmask(ji,jj,1)   ! Long  Wave
+         END DO
+      END DO
 
 
@@ -461,6 +500,13 @@
       END IF
 
-      Cd_oce(:,:) = Cd(:,:)  ! record value of pure ocean-atm. drag (clem)
-
+!$OMP PARALLEL
+!$OMP DO schedule(static) private(jj, ji)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            Cd_oce(ji,jj) = Cd(ji,jj)  ! record value of pure ocean-atm. drag (clem)
+         END DO
+      END DO
+
+!$OMP DO schedule(static) private(jj, ji)
       DO jj = 1, jpj             ! tau module, i and j component
          DO ji = 1, jpi
@@ -471,7 +517,15 @@
          END DO
       END DO
+!$OMP END PARALLEL
 
       !                          ! add the HF tau contribution to the wind stress module
-      IF( lhftau )   taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1)
+      IF( lhftau ) THEN
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               taum(ji,jj) = taum(ji,jj) + sf(jp_tdif)%fnow(ji,jj,1)
+            END DO
+         END DO
+      END IF
 
       CALL iom_put( "taum_oce", taum )   ! output wind stress module
@@ -480,4 +534,5 @@
       !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
       !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
       DO jj = 1, jpjm1
          DO ji = 1, fs_jpim1
@@ -496,14 +551,29 @@
 
       ! zqla used as temporary array, for rho*U (common term of bulk formulae):
-      zqla(:,:) = zrhoa(:,:) * zU_zu(:,:)
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zqla(ji,jj) = zrhoa(ji,jj) * zU_zu(ji,jj)
+         END DO
+      END DO
 
       IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN
          !! q_air and t_air are given at 10m (wind reference height)
-         zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) ) ! Evaporation, using bulk wind speed
-         zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch(:,:)*(zst(:,:) - ztpot(:,:)             )   ! Sensible Heat, using bulk wind speed
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zevap(ji,jj) = rn_efac*MAX( 0._wp,             zqla(ji,jj)*Ce(ji,jj)*(zsq(ji,jj) - sf(jp_humi)%fnow(ji,jj,1)) ) ! Evaporation, using bulk wind speed
+            END DO
+         END DO
+         zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch(:,:)*(zst(:,:) - ztpot(:,:) )   ! Sensible Heat, using bulk wind speed
       ELSE
          !! q_air and t_air are not given at 10m (wind reference height)
          ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!!
-         zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce(:,:)*(zsq(:,:) - zq_zu(:,:) ) ) ! Evaporation ! using bulk wind speed
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zevap(ji,jj) = rn_efac*MAX( 0._wp,             zqla(ji,jj)*Ce(ji,jj)*(zsq(ji,jj) - zq_zu(ji,jj) ) ) ! Evaporation ! using bulk wind speed
+            END DO
+         END DO
          zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch(:,:)*(zst(:,:) - zt_zu(:,:) )   ! Sensible Heat ! using bulk wind speed
       ENDIF
@@ -527,19 +597,24 @@
       ! ----------------------------------------------------------------------------- !
       !
-      emp (:,:) = (  zevap(:,:)                                          &   ! mass flux (evap. - precip.)
-         &         - sf(jp_prec)%fnow(:,:,1) * rn_pfac  ) * tmask(:,:,1)
-      !
-      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar
-         &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * lfus                         &   ! remove latent melting heat for solid precip
-         &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST
-         &     + ( sf(jp_prec)%fnow(:,:,1) - sf(jp_snow)%fnow(:,:,1) ) * rn_pfac  &   ! add liquid precip heat content at Tair
-         &     * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp                          &
-         &     + sf(jp_snow)%fnow(:,:,1) * rn_pfac                                &   ! add solid  precip heat content at min(Tair,Tsnow)
-         &     * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1)
-      !
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            emp (ji,jj) = (  zevap(ji,jj)                                          &   ! mass flux (evap. - precip.)
+               &         - sf(jp_prec)%fnow(ji,jj,1) * rn_pfac  ) * tmask(ji,jj,1)
+            !
+            qns(ji,jj) = zqlw(ji,jj) - zqsb(ji,jj) - zqla(ji,jj)                                &   ! Downward Non Solar
+               &     - sf(jp_snow)%fnow(ji,jj,1) * rn_pfac * lfus                         &   ! remove latent melting heat for solid precip
+               &     - zevap(ji,jj) * pst(ji,jj) * rcp                                      &   ! remove evap heat content at SST
+               &     + ( sf(jp_prec)%fnow(ji,jj,1) - sf(jp_snow)%fnow(ji,jj,1) ) * rn_pfac  &   ! add liquid precip heat content at Tair
+               &     * ( sf(jp_tair)%fnow(ji,jj,1) - rt0 ) * rcp                          &
+               &     + sf(jp_snow)%fnow(ji,jj,1) * rn_pfac                                &   ! add solid  precip heat content at min(Tair,Tsnow)
+               &     * ( MIN( sf(jp_tair)%fnow(ji,jj,1), rt0_snow ) - rt0 ) * cpic * tmask(ji,jj,1)
+            !
 #if defined key_lim3
-      qns_oce(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                ! non solar without emp (only needed by LIM3)
-      qsr_oce(:,:) = qsr(:,:)
+            qns_oce(ji,jj) = zqlw(ji,jj) - zqsb(ji,jj) - zqla(ji,jj)                                ! non solar without emp (only needed by LIM3)
+            qsr_oce(ji,jj) = qsr(ji,jj)
 #endif
+         END DO
+      END DO
       !
       IF ( nn_ice == 0 ) THEN
@@ -551,6 +626,11 @@
          CALL iom_put( "qsr_oce" ,   qsr  )                 ! output downward solar heat over the ocean
          CALL iom_put( "qt_oce"  ,   qns+qsr )              ! output total downward heat over the ocean
-         tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac   ! output total precipitation [kg/m2/s]
-         sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac   ! output solid precipitation [kg/m2/s]
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               tprecip(ji,jj) = sf(jp_prec)%fnow(ji,jj,1) * rn_pfac   ! output total precipitation [kg/m2/s]
+               sprecip(ji,jj) = sf(jp_snow)%fnow(ji,jj,1) * rn_pfac   ! output solid precipitation [kg/m2/s]
+            END DO
+         END DO
          CALL iom_put( 'snowpre', sprecip * 86400. )        ! Snow
          CALL iom_put( 'precip' , tprecip * 86400. )        ! Total precipitation
@@ -599,5 +679,10 @@
       CALL wrk_alloc( jpi,jpj, Cd )
 
-      Cd(:,:) = Cd_ice
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            Cd(ji,jj) = Cd_ice
+         END DO
+      END DO
 
       ! Make ice-atm. drag dependent on ice concentration (see Lupkes et al. 2012) (clem)
@@ -613,9 +698,14 @@
       zrhoa (:,:) =  rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1))
 
-      !!gm brutal....
-      utau_ice  (:,:) = 0._wp
-      vtau_ice  (:,:) = 0._wp
-      wndm_ice  (:,:) = 0._wp
-      !!gm end
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            !!gm brutal....
+            utau_ice  (ji,jj) = 0._wp
+            vtau_ice  (ji,jj) = 0._wp
+            wndm_ice  (ji,jj) = 0._wp
+            !!gm end
+         END DO
+      END DO
 
       ! ----------------------------------------------------------------------------- !
@@ -625,4 +715,5 @@
       CASE( 'I' )                  ! B-grid ice dynamics :   I-point (i.e. F-point with sea-ice indexation)
          !                           and scalar wind at T-point ( = | U10m - U_ice | ) (masked)
+!$OMP PARALLEL DO schedule(static) private(jj,ji,zwndi_f,zwndj_f,zwnorm_f,zwndi_t,zwndj_t)
          DO jj = 2, jpjm1
             DO ji = 2, jpim1   ! B grid : NO vector opt
@@ -649,4 +740,5 @@
          !
       CASE( 'C' )                  ! C-grid ice dynamics :   U & V-points (same as ocean)
+!$OMP PARALLEL DO schedule(static) private(jj,ji,zwndi_t,zwndj_t)
          DO jj = 2, jpj
             DO ji = fs_2, jpi   ! vect. opt.
@@ -656,4 +748,5 @@
             END DO
          END DO
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vect. opt.
@@ -700,8 +793,10 @@
       REAL(wp) ::   zztmp, z1_lsub           !   -      -
       REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qlw         ! long wave heat flux over ice
+      REAL(wp), DIMENSION(:,:,:), POINTER ::   zevap_ice3d, zqns_ice3d, zqsr_ice3d 
       REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qsb         ! sensible  heat flux over ice
       REAL(wp), DIMENSION(:,:,:), POINTER ::   z_dqlw        ! long wave heat sensitivity over ice
       REAL(wp), DIMENSION(:,:,:), POINTER ::   z_dqsb        ! sensible  heat sensitivity over ice
       REAL(wp), DIMENSION(:,:)  , POINTER ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (LIM3)
+      REAL(wp), DIMENSION(:,:)  , POINTER ::   zevap_ice2d, zqns_ice2d, zqsr_ice2d
       REAL(wp), DIMENSION(:,:)  , POINTER ::   zrhoa
       REAL(wp), DIMENSION(:,:)  , POINTER ::   Cd            ! transfer coefficient for momentum      (tau)
@@ -710,9 +805,14 @@
       IF( nn_timing == 1 )  CALL timing_start('blk_ice_flx')
       !
-      CALL wrk_alloc( jpi,jpj,jpl,   z_qlw, z_qsb, z_dqlw, z_dqsb )
-      CALL wrk_alloc( jpi,jpj,       zrhoa)
+      CALL wrk_alloc( jpi,jpj,jpl,   z_qlw, z_qsb, z_dqlw, z_dqsb, zevap_ice3d, zqns_ice3d, zqsr_ice3d )
+      CALL wrk_alloc( jpi,jpj,       zrhoa, zevap_ice2d, zqns_ice2d, zqsr_ice2d)
       CALL wrk_alloc( jpi,jpj, Cd )
 
-      Cd(:,:) = Cd_ice
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            Cd(ji,jj) = Cd_ice
+         END DO
+      END DO
 
       ! Make ice-atm. drag dependent on ice concentration (see Lupkes et al.  2012) (clem)
@@ -731,7 +831,8 @@
       !
       zztmp = 1. / ( 1. - albo )
-      !                                     ! ========================== !
-      DO jl = 1, jpl                        !  Loop over ice categories  !
-         !                                  ! ========================== !
+!$OMP PARALLEL
+!$OMP DO schedule(static) private(jl,jj,ji,zst2,zst3)            ! ========================== !
+      DO jl = 1, jpl                                             !  Loop over ice categories  !
+         !                                                       ! ========================== !
          DO jj = 1 , jpj
             DO ji = 1, jpi
@@ -781,6 +882,12 @@
       END DO
       !
-      tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac      ! total precipitation [kg/m2/s]
-      sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac      ! solid precipitation [kg/m2/s]
+!$OMP DO schedule(static) private(jj, ji)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            tprecip(ji,jj) = sf(jp_prec)%fnow(ji,jj,1) * rn_pfac      ! total precipitation [kg/m2/s]
+            sprecip(ji,jj) = sf(jp_snow)%fnow(ji,jj,1) * rn_pfac      ! solid precipitation [kg/m2/s]
+         END DO
+      END DO
+!$OMP END PARALLEL
       CALL iom_put( 'snowpre', sprecip * 86400. )                  ! Snow precipitation
       CALL iom_put( 'precip' , tprecip * 86400. )                  ! Total precipitation
@@ -791,35 +898,98 @@
       ! --- evaporation --- !
       z1_lsub = 1._wp / Lsub
-      evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_lsub    ! sublimation
-      devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_lsub    ! d(sublimation)/dT
-      zevap    (:,:)   = rn_efac * ( emp(:,:) + tprecip(:,:) )  ! evaporation over ocean
-
-      ! --- evaporation minus precipitation --- !
-      zsnw(:,:) = 0._wp
+!$OMP PARALLEL
+!$OMP DO schedule(static) private(jl,jj,ji)
+      DO jl = 1, jpl
+         DO jj = 1 , jpj
+            DO ji = 1, jpi
+               evap_ice (ji,jj,jl) = rn_efac * qla_ice (ji,jj,jl) * z1_lsub    ! sublimation
+               devap_ice(ji,jj,jl) = rn_efac * dqla_ice(ji,jj,jl) * z1_lsub    ! d(sublimation)/dT
+            END DO
+         END DO
+      END DO
+      !
+!$OMP DO schedule(static) private(jj, ji)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zevap    (ji,jj)   = rn_efac * ( emp(ji,jj) + tprecip(ji,jj) )  ! evaporation over ocean
+
+            ! --- evaporation minus precipitation --- !
+            zsnw(ji,jj) = 0._wp
+         END DO
+      END DO
+!$OMP END PARALLEL
       CALL lim_thd_snwblow( pfrld, zsnw )  ! snow distribution over ice after wind blowing
-      emp_oce(:,:) = pfrld(:,:) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw )
-      emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw
-      emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:)
-
-      ! --- heat flux associated with emp --- !
-      qemp_oce(:,:) = - pfrld(:,:) * zevap(:,:) * sst_m(:,:) * rcp                               & ! evap at sst
-         &          + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp  & ! liquid precip at Tair
-         &          +   sprecip(:,:) * ( 1._wp - zsnw ) *                                        & ! solid precip at min(Tair,Tsnow)
-         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
-      qemp_ice(:,:) =   sprecip(:,:) * zsnw *                                                    & ! solid precip (only)
-         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
-
-      ! --- total solar and non solar fluxes --- !
-      qns_tot(:,:) = pfrld(:,:) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 ) + qemp_ice(:,:) + qemp_oce(:,:)
-      qsr_tot(:,:) = pfrld(:,:) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 )
-
-      ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
-      qprec_ice(:,:) = rhosn * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
+!$OMP PARALLEL
+!$OMP DO schedule(static) private(jj,ji)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            emp_oce(ji,jj) = pfrld(ji,jj) * zevap(ji,jj) - ( tprecip(ji,jj) - sprecip(ji,jj) ) - sprecip(ji,jj) * (1._wp - zsnw(ji,jj))
+         END DO
+      END DO
+!$OMP END DO NOWAIT
+!$OMP DO schedule(static) private(jl,jj,ji)
+      DO jl = 1, jpl
+         DO jj = 1 , jpj
+            DO ji = 1, jpi
+               zevap_ice3d(ji,jj,jl) = a_i_b(ji,jj,jl) * evap_ice(ji,jj,jl)
+               zqns_ice3d(ji,jj,jl) = a_i_b(ji,jj,jl) * qns_ice(ji,jj,jl)
+               zqsr_ice3d(ji,jj,jl) = a_i_b(ji,jj,jl) * qsr_ice(ji,jj,jl)
+            END DO
+         END DO
+      END DO
+!$OMP END DO NOWAIT
+!$OMP DO schedule(static) private(jj,ji)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zevap_ice2d(ji,jj) = 0._wp 
+            zqns_ice2d(ji,jj) = 0._wp
+            zqsr_ice2d(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
+               zevap_ice2d(ji,jj) = zevap_ice2d(ji,jj) + zevap_ice3d(ji,jj,jl)
+               zqns_ice2d(ji,jj) = zqns_ice2d(ji,jj) + zqns_ice3d(ji,jj,jl)
+               zqsr_ice2d(ji,jj) = zqsr_ice2d(ji,jj) + zqsr_ice3d(ji,jj,jl)
+            END DO
+         END DO
+      END DO
+!$OMP DO schedule(static) private(jj,ji)
+      DO jj = 1 , jpj
+         DO ji = 1, jpi
+            emp_ice(ji,jj) = zevap_ice2d(ji,jj) - sprecip(ji,jj) * zsnw(ji,jj)
+            emp_tot(ji,jj) = emp_oce(ji,jj) + emp_ice(ji,jj)
+
+            ! --- heat flux associated with emp --- !
+            qemp_oce(ji,jj) = - pfrld(ji,jj) * zevap(ji,jj) * sst_m(ji,jj) * rcp                                & ! evap at sst
+               &          + ( tprecip(ji,jj) - sprecip(ji,jj) ) * ( sf(jp_tair)%fnow(ji,jj,1) - rt0 ) * rcp     & ! liquid precip at Tair
+               &          +   sprecip(ji,jj) * ( 1._wp - zsnw(ji,jj) ) *                                        & ! solid precip at min(Tair,Tsnow)
+               &              ( ( MIN( sf(jp_tair)%fnow(ji,jj,1), rt0_snow ) - rt0 ) * cpic * tmask(ji,jj,1) - lfus )
+            qemp_ice(ji,jj) =   sprecip(ji,jj) * zsnw(ji,jj) *                                                  & ! solid precip (only)
+               &              ( ( MIN( sf(jp_tair)%fnow(ji,jj,1), rt0_snow ) - rt0 ) * cpic * tmask(ji,jj,1) - lfus )
+
+            ! --- total solar and non solar fluxes --- !
+            qns_tot(ji,jj) = pfrld(ji,jj) * qns_oce(ji,jj) + zqns_ice2d(ji,jj) + qemp_ice(ji,jj) + qemp_oce(ji,jj)
+            qsr_tot(ji,jj) = pfrld(ji,jj) * qsr_oce(ji,jj) + zqsr_ice2d(ji,jj)
+
+            ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
+            qprec_ice(ji,jj) = rhosn * ( ( MIN( sf(jp_tair)%fnow(ji,jj,1), rt0_snow ) - rt0 ) * cpic * tmask(ji,jj,1) - lfus )
+         END DO
+      END DO
+!$OMP END DO NOWAIT
 
       ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) ---
+!$OMP DO schedule(static) private(jl,jj,ji)
       DO jl = 1, jpl
-         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * cpic * tmask(:,:,1) )
-                                   ! But we do not have Tice => consider it at 0degC => evap=0 
-      END DO
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               qevap_ice(ji,jj,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * cpic * tmask(:,:,1) )
+                                           ! But we do not have Tice => consider it at 0degC => evap=0 
+            END DO
+         END DO
+      END DO
+!$OMP END PARALLEL
 
       CALL wrk_dealloc( jpi,jpj,   zevap, zsnw )
@@ -831,6 +1001,11 @@
       ! ( Maykut and Untersteiner, 1971 ; Ebert and Curry, 1993 )
       !
-      fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
-      fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            fr1_i0(ji,jj) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
+            fr2_i0(ji,jj) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
+         END DO
+      END DO
       !
       !
@@ -844,7 +1019,7 @@
       ENDIF
 
-      CALL wrk_dealloc( jpi,jpj,jpl,   z_qlw, z_qsb, z_dqlw, z_dqsb )
+      CALL wrk_dealloc( jpi,jpj,jpl,   z_qlw, z_qsb, z_dqlw, z_dqsb, zevap_ice3d, zqns_ice3d, zqsr_ice3d )
       CALL wrk_dealloc( jpi,jpj,       zrhoa )
-      CALL wrk_dealloc( jpi,jpj, Cd )
+      CALL wrk_dealloc( jpi,jpj, Cd, zevap_ice2d, zqns_ice2d, zqsr_ice2d)
       !
       IF( nn_timing == 1 )  CALL timing_stop('blk_ice_flx')
@@ -908,4 +1083,5 @@
       !!----------------------------------------------------------------------------------
       !
+!$OMP PARALLEL DO schedule(static) private(jj,ji,ztmp,ze_sat)
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -944,4 +1120,5 @@
       !!----------------------------------------------------------------------------------
       !
+!$OMP PARALLEL DO schedule(static) private(jj,ji,zrv,ziRT)
       DO jj = 1, jpj
          DO ji = 1, jpi
Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_ncar.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_ncar.F90	(revision 7646)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_ncar.F90	(revision 7698)
@@ -114,4 +114,5 @@
       !
       INTEGER ::   j_itt
+      INTEGER ::   ji, jj             ! dummy loop indices
       LOGICAL ::   l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U
       INTEGER , PARAMETER ::   nb_itt = 4       ! number of itterations
@@ -141,6 +142,11 @@
       !! Neutral coefficients at 10m:
       IF( ln_cdgw ) THEN      ! wave drag case
-         cdn_wave(:,:) = cdn_wave(:,:) + rsmall * ( 1._wp - tmask(:,:,1) )
-         ztmp0   (:,:) = cdn_wave(:,:)
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               cdn_wave(ji,jj) = cdn_wave(ji,jj) + rsmall * ( 1._wp - tmask(ji,jj,1) )
+               ztmp0   (ji,jj) = cdn_wave(ji,jj)
+            END DO
+         END DO
       ELSE
          ztmp0 = cd_neutral_10m( U_blk )
@@ -245,4 +251,5 @@
       !!----------------------------------------------------------------------------------
       !
+!$OMP PARALLEL DO schedule(static) private(jj,ji,zw,zw6,zgt33)
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -284,4 +291,5 @@
       !!----------------------------------------------------------------------------------
       !
+!$OMP PARALLEL DO schedule(static) private(jj,ji,zx2,zx,zstab)
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -318,4 +326,5 @@
       !!----------------------------------------------------------------------------------
       !
+!$OMP PARALLEL DO schedule(static) private(jj,ji,zx2,zstab)
       DO jj = 1, jpj
          DO ji = 1, jpi
Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 7646)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 7698)
@@ -109,5 +109,5 @@
                                        !                    4 = Pure Coupled formulation)
       !!
-      INTEGER  ::   jl                 ! dummy loop index
+      INTEGER  ::   jl, jj, ji         ! dummy loop index
       REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_os, zalb_cs  ! ice albedo under overcast/clear sky
       REAL(wp), POINTER, DIMENSION(:,:  )   ::   zutau_ice, zvtau_ice 
@@ -133,15 +133,32 @@
 
          ! mean surface ocean current at ice velocity point (C-grid dynamics :  U- & V-points as the ocean)
-         u_oce(:,:) = ssu_m(:,:) * umask(:,:,1)
-         v_oce(:,:) = ssv_m(:,:) * vmask(:,:,1)
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               u_oce(ji,jj) = ssu_m(ji,jj) * umask(ji,jj,1)
+               v_oce(ji,jj) = ssv_m(ji,jj) * vmask(ji,jj,1)
+            END DO
+         END DO
 
          ! masked sea surface freezing temperature [Kelvin] (set to rt0 over land)
          CALL eos_fzp( sss_m(:,:) , t_bo(:,:) )
-         t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )
+!$OMP PARALLEL
+!$OMP DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               t_bo(ji,jj) = ( t_bo(ji,jj) + rt0 ) * tmask(ji,jj,1) + rt0 * ( 1._wp - tmask(ji,jj,1) )
+            END DO
+         END DO
 
          ! Mask sea ice surface temperature (set to rt0 over land)
          DO jl = 1, jpl
-            t_su(:,:,jl) = t_su(:,:,jl) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )
-         END DO
+!$OMP DO schedule(static) private(jj, ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  t_su(ji,jj,jl) = t_su(ji,jj,jl) * tmask(ji,jj,1) + rt0 * ( 1._wp - tmask(ji,jj,1) )
+               END DO
+            END DO
+         END DO
+!$OMP END PARALLEL
          !
          !------------------------------------------------!
@@ -161,6 +178,11 @@
             CALL wrk_alloc( jpi,jpj    , zutau_ice, zvtau_ice)
                                       CALL sbc_cpl_ice_tau( zutau_ice , zvtau_ice )
-            utau_ice(:,:) = utau_ice(:,:) * xcplmask(:,:,0) + zutau_ice(:,:) * ( 1. - xcplmask(:,:,0) )
-            vtau_ice(:,:) = vtau_ice(:,:) * xcplmask(:,:,0) + zvtau_ice(:,:) * ( 1. - xcplmask(:,:,0) )
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  utau_ice(ji,jj) = utau_ice(ji,jj) * xcplmask(ji,jj,0) + zutau_ice(ji,jj) * ( 1. - xcplmask(ji,jj,0) )
+                  vtau_ice(ji,jj) = vtau_ice(ji,jj) * xcplmask(ji,jj,0) + zvtau_ice(ji,jj) * ( 1. - xcplmask(ji,jj,0) )
+               END DO
+            END DO
             CALL wrk_dealloc( jpi,jpj  , zutau_ice, zvtau_ice)
          ENDIF
@@ -180,6 +202,11 @@
                                       CALL lim_dyn( kt )       !     rheology  
             ELSE
-               u_ice(:,:) = rn_uice * umask(:,:,1)             !     or prescribed velocity
-               v_ice(:,:) = rn_vice * vmask(:,:,1)
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     u_ice(ji,jj) = rn_uice * umask(ji,jj,1)             !     or prescribed velocity
+                     v_ice(ji,jj) = rn_vice * vmask(ji,jj,1)
+                  END DO
+               END DO
             ENDIF
                                       CALL lim_trp( kt )       ! -- Ice transport (Advection/diffusion)
@@ -200,6 +227,11 @@
                                       CALL lim_var_agg(1)      ! at_i for coupling (via pfrld) 
          !
-         pfrld(:,:)   = 1._wp - at_i(:,:)
-         phicif(:,:)  = vt_i(:,:)
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               pfrld(ji,jj)   = 1._wp - at_i(ji,jj)
+               phicif(ji,jj)  = vt_i(ji,jj)
+            END DO
+         END DO
 
          !------------------------------------------------------!
@@ -220,5 +252,12 @@
             CASE( jp_blk )                                          ! bulk formulation
                ! albedo depends on cloud fraction because of non-linear spectral effects
-               alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
+               DO jl = 1, jpl
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+                  DO jj = 1, jpj
+                     DO ji = 1, jpi
+                        alb_ice(ji,jj,jl) = ( 1. - cldf_ice ) * zalb_cs(ji,jj,jl) + cldf_ice * zalb_os(ji,jj,jl)
+                     END DO
+                  END DO
+               END DO
                                       CALL blk_ice_flx( t_su, alb_ice )
                IF( ln_mixcpl      )   CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su )
@@ -226,5 +265,12 @@
             CASE ( jp_purecpl )
                ! albedo depends on cloud fraction because of non-linear spectral effects
-               alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
+               DO jl = 1, jpl
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+                  DO jj = 1, jpj
+                     DO ji = 1, jpi
+                        alb_ice(ji,jj,jl) = ( 1. - cldf_ice ) * zalb_cs(ji,jj,jl) + cldf_ice * zalb_os(ji,jj,jl)
+                     END DO
+                  END DO
+               END DO
                                       CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su )
                IF( nn_limflx == 2 )   CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
@@ -285,5 +331,5 @@
       !! ** purpose :   Allocate all the dynamic arrays of the LIM-3 modules
       !!----------------------------------------------------------------------
-      INTEGER :: ji, jj, ierr
+      INTEGER :: jl, ji, jj, ierr
       !!----------------------------------------------------------------------
       IF(lwp) WRITE(numout,*)
@@ -334,7 +380,23 @@
       IF( ln_limdiahsb) CALL lim_diahsb_init  ! initialization for diags
       !
-      fr_i(:,:)     = at_i(:,:)         ! initialisation of sea-ice fraction
-      tn_ice(:,:,:) = t_su(:,:,:)       ! initialisation of surface temp for coupled simu
-      !
+!$OMP PARALLEL
+!$OMP DO schedule(static) private(jj, ji)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            fr_i(ji,jj)     = at_i(ji,jj)         ! initialisation of sea-ice fraction
+         END DO
+      END DO
+!$OMP END DO NOWAIT
+      DO jl = 1, jpl
+!$OMP DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               tn_ice(ji,jj,jl) = t_su(ji,jj,jl)       ! initialisation of surface temp for coupled simu
+            END DO
+         END DO
+!$OMP END DO NOWAIT
+      END DO
+      !
+!$OMP DO schedule(static) private(jj, ji)
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -344,4 +406,5 @@
          END DO
       END DO
+!$OMP END PARALLEL
       !
       nstart = numit  + nn_fsbc
@@ -527,5 +590,5 @@
       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdevap_ice ! sublimation sensitivity
       !
-      INTEGER  ::   jl      ! dummy loop index
+      INTEGER  ::   jl, jj, ji      ! dummy loop index
       !
       REAL(wp), POINTER, DIMENSION(:,:) :: zalb_m    ! Mean albedo over all categories
@@ -550,14 +613,28 @@
          z_evap_m (:,:) = fice_ice_ave ( pevap_ice (:,:,:) )
          z_devap_m(:,:) = fice_ice_ave ( pdevap_ice (:,:,:) )
+
+!$OMP PARALLEL
          DO jl = 1, jpl
-            pdqn_ice  (:,:,jl) = z_dqn_m(:,:)
-            pdevap_ice(:,:,jl) = z_devap_m(:,:)
+!$OMP DO schedule(static) private(jj, ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  pdqn_ice  (ji,jj,jl) = z_dqn_m(ji,jj)
+                  pdevap_ice(ji,jj,jl) = z_devap_m(ji,jj)
+               END DO
+            END DO
+!$OMP END DO NOWAIT
          END DO
          !
          DO jl = 1, jpl
-            pqns_ice (:,:,jl) = z_qns_m(:,:)
-            pqsr_ice (:,:,jl) = z_qsr_m(:,:)
-            pevap_ice(:,:,jl) = z_evap_m(:,:)
-         END DO
+!$OMP DO schedule(static) private(jj, ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  pqns_ice (ji,jj,jl) = z_qns_m(ji,jj)
+                  pqsr_ice (ji,jj,jl) = z_qsr_m(ji,jj)
+                  pevap_ice(ji,jj,jl) = z_evap_m(ji,jj)
+               END DO
+            END DO
+         END DO
+!$OMP END PARALLEL
          !
          CALL wrk_dealloc( jpi,jpj, z_qsr_m, z_qns_m, z_evap_m, z_dqn_m, z_devap_m)
@@ -571,7 +648,12 @@
          ztem_m(:,:) = fice_ice_ave ( ptn_ice  (:,:,:) )
          DO jl = 1, jpl
-            pqns_ice (:,:,jl) = pqns_ice (:,:,jl) + pdqn_ice  (:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) )
-            pevap_ice(:,:,jl) = pevap_ice(:,:,jl) + pdevap_ice(:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) )
-            pqsr_ice (:,:,jl) = pqsr_ice (:,:,jl) * ( 1._wp - palb_ice(:,:,jl) ) / ( 1._wp - zalb_m(:,:) )
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  pqns_ice (ji,jj,jl) = pqns_ice (ji,jj,jl) + pdqn_ice  (ji,jj,jl) * ( ptn_ice(ji,jj,jl) - ztem_m(ji,jj) )
+                  pevap_ice(ji,jj,jl) = pevap_ice(ji,jj,jl) + pdevap_ice(ji,jj,jl) * ( ptn_ice(ji,jj,jl) - ztem_m(ji,jj) )
+                  pqsr_ice (ji,jj,jl) = pqsr_ice (ji,jj,jl) * ( 1._wp - palb_ice(ji,jj,jl) ) / ( 1._wp - zalb_m(ji,jj) )
+               END DO
+            END DO
          END DO
          !
@@ -590,15 +672,60 @@
       !! ** purpose :  store ice variables at "before" time step
       !!----------------------------------------------------------------------
-      a_i_b  (:,:,:)   = a_i  (:,:,:)     ! ice area
-      e_i_b  (:,:,:,:) = e_i  (:,:,:,:)   ! ice thermal energy
-      v_i_b  (:,:,:)   = v_i  (:,:,:)     ! ice volume
-      v_s_b  (:,:,:)   = v_s  (:,:,:)     ! snow volume
-      e_s_b  (:,:,:,:) = e_s  (:,:,:,:)   ! snow thermal energy
-      smv_i_b(:,:,:)   = smv_i(:,:,:)     ! salt content
-      oa_i_b (:,:,:)   = oa_i (:,:,:)     ! areal age content
-      u_ice_b(:,:)     = u_ice(:,:)
-      v_ice_b(:,:)     = v_ice(:,:)
-      !
-      at_i_b (:,:)     = SUM( a_i_b(:,:,:), dim=3 )
+      INTEGER  ::   jn, jl, jj, ji         ! dummy loop index
+
+!$OMP PARALLEL
+      DO jl = 1, jpl
+!$OMP DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               a_i_b  (ji,jj,jl)   = a_i  (ji,jj,jl)     ! ice area
+               v_i_b  (ji,jj,jl)   = v_i  (ji,jj,jl)     ! ice volume
+               v_s_b  (ji,jj,jl)   = v_s  (ji,jj,jl)     ! snow volume
+               smv_i_b(ji,jj,jl)   = smv_i(ji,jj,jl)     ! salt content
+               oa_i_b (ji,jj,jl)   = oa_i (ji,jj,jl)     ! areal age content
+            END DO
+         END DO
+!$OMP END DO NOWAIT
+      END DO
+      DO jl = 1, jpl
+         DO jn = 1, nlay_i
+!$OMP DO schedule(static) private(jj, ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  e_i_b  (ji,jj,jn,jl) = e_i  (ji,jj,jn,jl)   ! ice thermal energy
+               END DO
+            END DO
+!$OMP END DO NOWAIT
+         END DO
+      END DO
+      DO jl = 1, jpl
+         DO jn = 1, nlay_s
+!$OMP DO schedule(static) private(jj, ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  e_s_b  (ji,jj,jn,jl) = e_s  (ji,jj,jn,jl)   ! snow thermal energy
+               END DO
+            END DO
+!$OMP END DO NOWAIT
+         END DO
+      END DO
+!$OMP DO schedule(static) private(jj, ji)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            u_ice_b(ji,jj)     = u_ice(ji,jj)
+            v_ice_b(ji,jj)     = v_ice(ji,jj)
+            at_i_b (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_b (ji,jj)     = at_i_b (ji,jj) + a_i_b(ji,jj,jl)
+            END DO
+         END DO
+      END DO
+!$OMP END PARALLEL
       
    END SUBROUTINE sbc_lim_bef
@@ -612,35 +739,42 @@
       !!               of the time step
       !!----------------------------------------------------------------------
-      sfx    (:,:) = 0._wp   ;
-      sfx_bri(:,:) = 0._wp   ;   sfx_lam(:,:) = 0._wp
-      sfx_sni(:,:) = 0._wp   ;   sfx_opw(:,:) = 0._wp
-      sfx_bog(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp
-      sfx_bom(:,:) = 0._wp   ;   sfx_sum(:,:) = 0._wp
-      sfx_res(:,:) = 0._wp   ;   sfx_sub(:,:) = 0._wp
-      !
-      wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp
-      wfx_sni(:,:) = 0._wp   ;   wfx_opw(:,:) = 0._wp
-      wfx_bog(:,:) = 0._wp   ;   wfx_dyn(:,:) = 0._wp
-      wfx_bom(:,:) = 0._wp   ;   wfx_sum(:,:) = 0._wp
-      wfx_res(:,:) = 0._wp   ;   wfx_sub(:,:) = 0._wp
-      wfx_spr(:,:) = 0._wp   ;   wfx_lam(:,:) = 0._wp  
+      INTEGER  ::   jj, ji         ! dummy loop index
+
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            sfx    (ji,jj) = 0._wp   ;
+            sfx_bri(ji,jj) = 0._wp   ;   sfx_lam(ji,jj) = 0._wp
+            sfx_sni(ji,jj) = 0._wp   ;   sfx_opw(ji,jj) = 0._wp
+            sfx_bog(ji,jj) = 0._wp   ;   sfx_dyn(ji,jj) = 0._wp
+            sfx_bom(ji,jj) = 0._wp   ;   sfx_sum(ji,jj) = 0._wp
+            sfx_res(ji,jj) = 0._wp   ;   sfx_sub(ji,jj) = 0._wp
+            !
+            wfx_snw(ji,jj) = 0._wp   ;   wfx_ice(ji,jj) = 0._wp
+            wfx_sni(ji,jj) = 0._wp   ;   wfx_opw(ji,jj) = 0._wp
+            wfx_bog(ji,jj) = 0._wp   ;   wfx_dyn(ji,jj) = 0._wp
+            wfx_bom(ji,jj) = 0._wp   ;   wfx_sum(ji,jj) = 0._wp
+            wfx_res(ji,jj) = 0._wp   ;   wfx_sub(ji,jj) = 0._wp
+            wfx_spr(ji,jj) = 0._wp   ;   wfx_lam(ji,jj) = 0._wp  
       
-      hfx_thd(:,:) = 0._wp   ;
-      hfx_snw(:,:) = 0._wp   ;   hfx_opw(:,:) = 0._wp
-      hfx_bog(:,:) = 0._wp   ;   hfx_dyn(:,:) = 0._wp
-      hfx_bom(:,:) = 0._wp   ;   hfx_sum(:,:) = 0._wp
-      hfx_res(:,:) = 0._wp   ;   hfx_sub(:,:) = 0._wp
-      hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp
-      hfx_err(:,:) = 0._wp   ;   hfx_err_rem(:,:) = 0._wp
-      hfx_err_dif(:,:) = 0._wp
-      wfx_err_sub(:,:) = 0._wp
-      !
-      afx_tot(:,:) = 0._wp   ;
-      afx_dyn(:,:) = 0._wp   ;   afx_thd(:,:) = 0._wp
-      !
-      diag_heat(:,:) = 0._wp ;   diag_smvi(:,:) = 0._wp
-      diag_vice(:,:) = 0._wp ;   diag_vsnw(:,:) = 0._wp
-
-      tau_icebfr(:,:) = 0._wp; ! landfast ice param only (clem: important to keep the init here)
+            hfx_thd(ji,jj) = 0._wp   ;
+            hfx_snw(ji,jj) = 0._wp   ;   hfx_opw(ji,jj) = 0._wp
+            hfx_bog(ji,jj) = 0._wp   ;   hfx_dyn(ji,jj) = 0._wp
+            hfx_bom(ji,jj) = 0._wp   ;   hfx_sum(ji,jj) = 0._wp
+            hfx_res(ji,jj) = 0._wp   ;   hfx_sub(ji,jj) = 0._wp
+            hfx_spr(ji,jj) = 0._wp   ;   hfx_dif(ji,jj) = 0._wp
+            hfx_err(ji,jj) = 0._wp   ;   hfx_err_rem(ji,jj) = 0._wp
+            hfx_err_dif(ji,jj) = 0._wp
+            wfx_err_sub(ji,jj) = 0._wp
+            !
+            afx_tot(ji,jj) = 0._wp   ;
+            afx_dyn(ji,jj) = 0._wp   ;   afx_thd(ji,jj) = 0._wp
+            !
+            diag_heat(ji,jj) = 0._wp ;   diag_smvi(ji,jj) = 0._wp
+            diag_vice(ji,jj) = 0._wp ;   diag_vsnw(ji,jj) = 0._wp
+      
+            tau_icebfr(ji,jj) = 0._wp; ! landfast ice param only (clem: important to keep the init here)
+         END DO
+      END DO
       
    END SUBROUTINE sbc_lim_diag0
Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 7646)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 7698)
@@ -84,4 +84,5 @@
       !!              - nsbc: type of sbc
       !!----------------------------------------------------------------------
+      INTEGER ::   ji, jj, jn                        ! dummy loop indices
       INTEGER ::   ios, icpt                         ! local integer
       LOGICAL ::   ll_purecpl, ll_opa, ll_not_nemo   ! local logical
@@ -240,15 +241,41 @@
       IF( .NOT.ln_isf ) THEN        !* No ice-shelf in the domain : allocate and set to zero
          IF( sbc_isf_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_isf arrays' )
-         fwfisf  (:,:)   = 0._wp   ;   risf_tsc  (:,:,:) = 0._wp
-         fwfisf_b(:,:)   = 0._wp   ;   risf_tsc_b(:,:,:) = 0._wp
+!$OMP PARALLEL
+!$OMP DO schedule(static) private(jj,ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               fwfisf  (ji,jj)   = 0.0_wp ; fwfisf_b  (ji,jj)   = 0.0_wp
+            END DO
+         END DO
+!$OMP END DO NOWAIT
+         DO jn = 1, jpts
+!$OMP DO schedule(static) private(jj,ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  risf_tsc(ji,jj,jn) = 0.0_wp ; risf_tsc_b(ji,jj,jn) = 0.0_wp
+               END DO
+            END DO
+         END DO
+!$OMP END PARALLEL
       END IF
       IF( nn_ice == 0 ) THEN        !* No sea-ice in the domain : ice fraction is always zero
-         IF( nn_components /= jp_iam_opa )   fr_i(:,:) = 0._wp    ! except for OPA in SAS-OPA coupled case
-      ENDIF
-      !
-      sfx   (:,:) = 0._wp           !* salt flux due to freezing/melting
-      fmmflx(:,:) = 0._wp           !* freezing minus melting flux
-
-      taum(:,:) = 0._wp             !* wind stress module (needed in GLS in case of reduced restart)
+         IF( nn_components /= jp_iam_opa ) THEN
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  fr_i(ji,jj) = 0._wp    ! except for OPA in SAS-OPA coupled case
+               END DO
+            END DO
+         END IF
+      ENDIF
+      !
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            sfx   (ji,jj) = 0._wp           !* salt flux due to freezing/melting
+            fmmflx(ji,jj) = 0._wp           !* freezing minus melting flux
+            taum  (ji,jj) = 0._wp           !* wind stress module (needed in GLS in case of reduced restart)
+         END DO
+      END DO
 
       !                          ! Choice of the Surface Boudary Condition (set nsbc)
@@ -356,4 +383,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time step
+      INTEGER ::   ji, jj, jn       ! dummy loop indices
       !
       LOGICAL ::   ll_sas, ll_opa   ! local logical
@@ -365,12 +393,32 @@
       IF( kt /= nit000 ) THEN                      !          Swap of forcing fields          !
          !                                         ! ---------------------------------------- !
-         utau_b(:,:) = utau(:,:)                         ! Swap the ocean forcing fields
-         vtau_b(:,:) = vtau(:,:)                         ! (except at nit000 where before fields
-         qns_b (:,:) = qns (:,:)                         !  are set at the end of the routine)
-         emp_b (:,:) = emp (:,:)
-         sfx_b (:,:) = sfx (:,:)
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               utau_b(ji,jj) = utau(ji,jj)                         ! Swap the ocean forcing fields
+               vtau_b(ji,jj) = vtau(ji,jj)                         ! (except at nit000 where before fields
+               qns_b (ji,jj) = qns (ji,jj)                         !  are set at the end of the routine)
+               emp_b (ji,jj) = emp (ji,jj)
+               sfx_b (ji,jj) = sfx (ji,jj)
+            END DO
+         END DO
          IF ( ln_rnf ) THEN
-            rnf_b    (:,:  ) = rnf    (:,:  )
-            rnf_tsc_b(:,:,:) = rnf_tsc(:,:,:)
+!$OMP PARALLEL
+!$OMP DO schedule(static) private(jj,ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  rnf_b    (ji,jj  ) = rnf    (ji,jj  )
+               END DO
+            END DO
+!$OMP END DO NOWAIT
+            DO jn = 1, jpts
+!$OMP DO schedule(static) private(jj,ji)
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     rnf_tsc_b(ji,jj,jn) = rnf_tsc(ji,jj,jn)
+                  END DO
+               END DO
+            END DO
+!$OMP END PARALLEL
          ENDIF
       ENDIF
@@ -401,7 +449,12 @@
       END SELECT
       IF ( ln_wave .AND. ln_tauoc) THEN                                 ! Wave stress subctracted
-            utau(:,:) = utau(:,:)*tauoc_wave(:,:)
-            vtau(:,:) = vtau(:,:)*tauoc_wave(:,:)
-            taum(:,:) = taum(:,:)*tauoc_wave(:,:)
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               utau(ji,jj) = utau(ji,jj)*tauoc_wave(ji,jj)
+               vtau(ji,jj) = vtau(ji,jj)*tauoc_wave(ji,jj)
+               taum(ji,jj) = taum(ji,jj)*tauoc_wave(ji,jj)
+            END DO
+         END DO
       !
             SELECT CASE( nsbc )
@@ -457,13 +510,23 @@
                CALL iom_get( numror, jpdom_autoglo, 'sfx_b', sfx_b )  ! before salt flux (T-point)
             ELSE
-               sfx_b (:,:) = sfx(:,:)
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     sfx_b (ji,jj) = sfx(ji,jj)
+                  END DO
+               END DO
             ENDIF
          ELSE                                                   !* no restart: set from nit000 values
             IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields set to nit000'
-            utau_b(:,:) = utau(:,:)
-            vtau_b(:,:) = vtau(:,:)
-            qns_b (:,:) = qns (:,:)
-            emp_b (:,:) = emp (:,:)
-            sfx_b (:,:) = sfx (:,:)
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  utau_b(ji,jj) = utau(ji,jj)
+                  vtau_b(ji,jj) = vtau(ji,jj)
+                  qns_b (ji,jj) = qns (ji,jj)
+                  emp_b (ji,jj) = emp(ji,jj)
+                  sfx_b (ji,jj) = sfx(ji,jj)
+               END DO
+            END DO
          ENDIF
       ENDIF
Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 7646)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 7698)
@@ -103,6 +103,6 @@
       INTEGER, INTENT(in) ::   kt          ! ocean time step
       !
-      INTEGER  ::   ji, jj    ! dummy loop indices
-      INTEGER  ::   z_err = 0 ! dummy integer for error handling
+      INTEGER  ::   ji, jj, jn    ! dummy loop indices
+      INTEGER  ::   z_err = 0     ! dummy integer for error handling
       !!----------------------------------------------------------------------
       REAL(wp), DIMENSION(:,:), POINTER       ::   ztfrz   ! freezing point used for temperature correction
@@ -120,22 +120,51 @@
       IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
          !
-         IF( .NOT. l_rnfcpl )   rnf(:,:) = rn_rfact * ( sf_rnf(1)%fnow(:,:,1) )       ! updated runoff value at time step kt
+         IF( .NOT. l_rnfcpl ) THEN                             ! updated runoff value at time step kt
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  rnf(ji,jj) = rn_rfact * ( sf_rnf(1)%fnow(ji,jj,1) )
+               END DO
+            END DO
+         END IF
          !
          !                                                     ! set temperature & salinity content of runoffs
          IF( ln_rnf_tem ) THEN                                       ! use runoffs temperature data
-            rnf_tsc(:,:,jp_tem) = ( sf_t_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rau0
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  rnf_tsc(ji,jj,jp_tem) = ( sf_t_rnf(1)%fnow(ji,jj,1) ) * rnf(ji,jj) * r1_rau0
+               END DO
+            END DO
             CALL eos_fzp( sss_m(:,:), ztfrz(:,:) )
-            WHERE( sf_t_rnf(1)%fnow(:,:,1) == -999._wp )             ! if missing data value use SST as runoffs temperature
-               rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rau0
-            END WHERE
-            WHERE( sf_t_rnf(1)%fnow(:,:,1) == -222._wp )             ! where fwf comes from melting of ice shelves or iceberg
-               rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rau0 - rnf(:,:) * rlfusisf * r1_rau0_rcp
-            END WHERE
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  IF ( sf_t_rnf(1)%fnow(ji,jj,1) == -999._wp ) THEN            ! if missing data value use SST as runoffs temperature
+                     rnf_tsc(ji,jj,jp_tem) = sst_m(ji,jj) * rnf(ji,jj) * r1_rau0
+                  END IF
+                  IF ( sf_t_rnf(1)%fnow(ji,jj,1) == -222._wp ) THEN            ! where fwf comes from melting of ice shelves or iceberg
+                     rnf_tsc(ji,jj,jp_tem) = ztfrz(ji,jj) * rnf(ji,jj) * r1_rau0 - rnf(ji,jj) * rlfusisf * r1_rau0_rcp
+                  END IF
+               END DO
+            END DO
          ELSE                                                        ! use SST as runoffs temperature
-            rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rau0
-         ENDIF
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  rnf_tsc(ji,jj,jp_tem) = sst_m(ji,jj) * rnf(ji,jj) * r1_rau0
+               END DO
+            END DO
+         END IF
          !                                                           ! use runoffs salinity data
-         IF( ln_rnf_sal )   rnf_tsc(:,:,jp_sal) = ( sf_s_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rau0
-         !                                                           ! else use S=0 for runoffs (done one for all in the init)
+         IF( ln_rnf_sal ) THEN
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  rnf_tsc(ji,jj,jp_sal) = ( sf_s_rnf(1)%fnow(ji,jj,1) ) * rnf(ji,jj) * r1_rau0
+               END DO
+            END DO
+         END IF
+         !                                                        ! else use S=0 for runoffs (done one for all in the init)
          CALL iom_put( "runoffs", rnf )         ! output runoffs arrays
       ENDIF
@@ -152,6 +181,21 @@
          ELSE                                                   !* no restart: set from nit000 values
             IF(lwp) WRITE(numout,*) '          nit000-1 runoff forcing fields set to nit000'
-            rnf_b    (:,:  ) = rnf    (:,:  )
-            rnf_tsc_b(:,:,:) = rnf_tsc(:,:,:)
+!$OMP PARALLEL
+!$OMP DO schedule(static) private(jj,ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  rnf_b    (ji,jj  ) = rnf    (ji,jj  )
+               END DO
+            END DO
+!$OMP END DO NOWAIT
+            DO jn = 1, jpts
+!$OMP DO schedule(static) private(jj,ji)
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     rnf_tsc_b(ji,jj,jn) = rnf_tsc(ji,jj,jn)
+                  END DO
+               END DO
+            END DO
+!$OMP END PARALLEL
          ENDIF
       ENDIF
@@ -187,5 +231,5 @@
       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
       !!
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
       REAL(wp) ::   zfact     ! local scalar
       !!----------------------------------------------------------------------
@@ -195,4 +239,5 @@
       IF( ln_rnf_depth .OR. ln_rnf_depth_ini ) THEN      !==   runoff distributed over several levels   ==!
          IF( ln_linssh ) THEN    !* constant volume case : just apply the runoff input flow
+!$OMP PARALLEL DO schedule(static) private(jj,ji,jk)
             DO jj = 1, jpj
                DO ji = 1, jpi
@@ -203,4 +248,5 @@
             END DO
          ELSE                    !* variable volume case
+!$OMP PARALLEL DO schedule(static) private(jj,ji,jk)
             DO jj = 1, jpj                   ! update the depth over which runoffs are distributed
                DO ji = 1, jpi
@@ -217,6 +263,11 @@
          ENDIF
       ELSE                       !==   runoff put only at the surface   ==!
-         h_rnf (:,:)   = e3t_n (:,:,1)        ! update h_rnf to be depth of top box
-         phdivn(:,:,1) = phdivn(:,:,1) - ( rnf(:,:) + rnf_b(:,:) ) * zfact * r1_rau0 / e3t_n(:,:,1)
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               h_rnf (ji,jj)   = e3t_n (ji,jj,1)        ! update h_rnf to be depth of top box
+               phdivn(ji,jj,1) = phdivn(ji,jj,1) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact * r1_rau0 / e3t_n(ji,jj,1)
+            END DO
+         END DO
       ENDIF
       !
@@ -235,5 +286,5 @@
       !!----------------------------------------------------------------------
       CHARACTER(len=32) ::   rn_dep_file   ! runoff file name
-      INTEGER           ::   ji, jj, jk, jm    ! dummy loop indices
+      INTEGER           ::   ji, jj, jk, jm, jn    ! dummy loop indices
       INTEGER           ::   ierror, inum  ! temporary integer
       INTEGER           ::   ios           ! Local integer output status for namelist read
@@ -256,8 +307,19 @@
          ln_rnf_mouth  = .FALSE.                   ! default definition needed for example by sbc_ssr or by tra_adv_muscl
          nkrnf         = 0
-         rnf     (:,:) = 0.0_wp
-         rnf_b   (:,:) = 0.0_wp
-         rnfmsk  (:,:) = 0.0_wp
-         rnfmsk_z(:)   = 0.0_wp
+!$OMP PARALLEL
+!$OMP DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               rnf     (ji,jj) = 0.0_wp
+               rnf_b   (ji,jj) = 0.0_wp
+               rnfmsk  (ji,jj) = 0.0_wp
+            END DO
+         END DO
+!$OMP END DO NOWAIT
+!$OMP DO schedule(static) private(jk)
+         DO jk = 1, jpk
+            rnfmsk_z(jk)   = 0.0_wp
+         END DO
+!$OMP END PARALLEL
          RETURN
       ENDIF
@@ -338,5 +400,12 @@
          CALL iom_close( inum )                                        ! close file
          !
-         nk_rnf(:,:) = 0                               ! set the number of level over which river runoffs are applied
+!$OMP PARALLEL
+!$OMP DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               nk_rnf(ji,jj) = 0                               ! set the number of level over which river runoffs are applied
+            END DO
+         END DO
+!$OMP DO schedule(static) private(jj, ji, jk)
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -354,4 +423,5 @@
             END DO
          END DO
+!$OMP DO schedule(static) private(jj, ji, jk)
          DO jj = 1, jpj                                ! set the associated depth
             DO ji = 1, jpi
@@ -362,4 +432,5 @@
             END DO
          END DO
+!$OMP END PARALLEL
          !
       ELSE IF( ln_rnf_depth_ini ) THEN           ! runoffs applied at the surface
@@ -381,10 +452,17 @@
          DEALLOCATE( zrnfcl )
          !
-         h_rnf(:,:) = 1.
-         !
          zacoef = rn_dep_max / rn_rnf_max            ! coef of linear relation between runoff and its depth (150m for max of runoff)
          !
-         WHERE( zrnf(:,:) > 0._wp )  h_rnf(:,:) = zacoef * zrnf(:,:)   ! compute depth for all runoffs
-         !
+!$OMP PARALLEL
+         IF( zrnf(ji,jj) > 0._wp ) THEN
+!$OMP DO schedule(static) private(jj, ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  h_rnf(ji,jj) = zacoef * zrnf(ji,jj)   ! compute depth for all runoffs
+               END DO
+            END DO
+         END IF
+         !
+!$OMP DO schedule(static) private(jj, ji, jk)
          DO jj = 1, jpj                     ! take in account min depth of ocean rn_hmin
             DO ji = 1, jpi
@@ -396,5 +474,11 @@
          END DO
          !
-         nk_rnf(:,:) = 0                       ! number of levels on which runoffs are distributed
+!$OMP DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               nk_rnf(ji,jj) = 0                       ! number of levels on which runoffs are distributed
+            END DO
+         END DO
+!$OMP DO schedule(static) private(jj, ji, jk)
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -409,7 +493,9 @@
             END DO
          END DO
+!$OMP END PARALLEL
          !
          DEALLOCATE( zrnf )
          !
+!$OMP PARALLEL DO schedule(static) private(jj, ji, jk)
          DO jj = 1, jpj                                ! set the associated depth
             DO ji = 1, jpi
@@ -428,10 +514,30 @@
          ENDIF
       ELSE                                       ! runoffs applied at the surface
-         nk_rnf(:,:) = 1
-         h_rnf (:,:) = e3t_n(:,:,1)
-      ENDIF
-      !
-      rnf(:,:) =  0._wp                         ! runoff initialisation
-      rnf_tsc(:,:,:) = 0._wp                    ! runoffs temperature & salinty contents initilisation
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               nk_rnf(ji,jj) = 1
+               h_rnf (ji,jj) = e3t_n(ji,jj,1)
+            END DO
+         END DO
+      ENDIF
+      !
+!$OMP PARALLEL
+!$OMP DO schedule(static) private(jj, ji)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            rnf(ji,jj) =  0._wp                         ! runoff initialisation
+         END DO
+      END DO
+!$OMP END DO NOWAIT
+      DO jn = 1, jpts
+!$OMP DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               rnf_tsc(ji,jj,jn) = 0._wp                    ! runoffs temperature & salinty contents initilisation
+            END DO
+         END DO
+      END DO
+!$OMP END PARALLEL
       !
       !                                   ! ========================
@@ -466,6 +572,17 @@
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) '          No specific treatment at river mouths'
-         rnfmsk  (:,:) = 0._wp
-         rnfmsk_z(:)   = 0._wp
+!$OMP PARALLEL
+!$OMP DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               rnfmsk  (ji,jj) = 0._wp
+            END DO
+         END DO
+!$OMP END DO NOWAIT
+!$OMP DO schedule(static) private(jk)
+         DO jk = 1, jpk
+            rnfmsk_z(jk)   = 0._wp
+         END DO
+!$OMP END PARALLEL
          nkrnf = 0
       ENDIF
Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90	(revision 7646)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90	(revision 7698)
@@ -59,4 +59,5 @@
       !
       !                                        !* surface T-, U-, V- ocean level variables (T, S, depth, velocity)
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -68,18 +69,52 @@
       IF( nn_fsbc == 1 ) THEN                             !   Instantaneous surface fields        !
          !                                                ! ---------------------------------------- !
-         ssu_m(:,:) = ub(:,:,1)
-         ssv_m(:,:) = vb(:,:,1)
-         IF( l_useCT )  THEN    ;   sst_m(:,:) = eos_pt_from_ct( zts(:,:,jp_tem), zts(:,:,jp_sal) )
-         ELSE                    ;   sst_m(:,:) = zts(:,:,jp_tem)
-         ENDIF
-         sss_m(:,:) = zts(:,:,jp_sal)
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ssu_m(ji,jj) = ub(ji,jj,1)
+               ssv_m(ji,jj) = vb(ji,jj,1)
+            END DO
+         END DO
+         IF( l_useCT )  THEN
+           sst_m(:,:) = eos_pt_from_ct( zts(:,:,jp_tem), zts(:,:,jp_sal) )
+         ELSE                    
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  sst_m(ji,jj) = zts(ji,jj,jp_tem)
+               END DO
+            END DO
+         ENDIF
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               sss_m(ji,jj) = zts(ji,jj,jp_sal)
+            END DO
+         END DO
          !                          ! removed inverse barometer ssh when Patm forcing is used (for sea-ice dynamics)
-         IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )
-         ELSE                    ;   ssh_m(:,:) = sshn(:,:)
-         ENDIF
-         !
-         e3t_m(:,:) = e3t_n(:,:,1)
-         !
-         frq_m(:,:) = fraqsr_1lev(:,:)
+         IF( ln_apr_dyn ) THEN  
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ssh_m(ji,jj) = sshn(ji,jj) - 0.5 * ( ssh_ib(ji,jj) + ssh_ibb(ji,jj) )
+               END DO
+            END DO
+         ELSE                    
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ssh_m(ji,jj) = sshn(ji,jj)
+               END DO
+            END DO
+         ENDIF
+         !
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               e3t_m(ji,jj) = e3t_n(ji,jj,1)
+         !
+               frq_m(ji,jj) = fraqsr_1lev(ji,jj)
+            END DO
+         END DO
          !
       ELSE
@@ -91,46 +126,119 @@
             IF(lwp) WRITE(numout,*) '~~~~~~~   '
             zcoef = REAL( nn_fsbc - 1, wp )
-            ssu_m(:,:) = zcoef * ub(:,:,1)
-            ssv_m(:,:) = zcoef * vb(:,:,1)
-            IF( l_useCT )  THEN    ;   sst_m(:,:) = zcoef * eos_pt_from_ct( zts(:,:,jp_tem), zts(:,:,jp_sal) )
-            ELSE                    ;   sst_m(:,:) = zcoef * zts(:,:,jp_tem)
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ssu_m(ji,jj) = zcoef * ub(ji,jj,1)
+                  ssv_m(ji,jj) = zcoef * vb(ji,jj,1)
+               END DO
+            END DO
+            IF( l_useCT )  THEN
+              sst_m(:,:) = zcoef * eos_pt_from_ct( zts(:,:,jp_tem), zts(:,:,jp_sal) )
+            ELSE                    
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+              DO jj = 1, jpj
+                 DO ji = 1, jpi
+                    sst_m(ji,jj) = zcoef * zts(ji,jj,jp_tem)
+                 END DO
+              END DO
             ENDIF
-            sss_m(:,:) = zcoef * zts(:,:,jp_sal)
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  sss_m(ji,jj) = zcoef * zts(ji,jj,jp_sal)
+               END DO
+            END DO
             !                          ! removed inverse barometer ssh when Patm forcing is used (for sea-ice dynamics)
-            IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = zcoef * ( sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) )
-            ELSE                    ;   ssh_m(:,:) = zcoef * sshn(:,:)
+            IF( ln_apr_dyn ) THEN   
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ssh_m(ji,jj) = zcoef * ( sshn(ji,jj) - 0.5 * ( ssh_ib(ji,jj) + ssh_ibb(ji,jj) ) )
+                  END DO
+               END DO
+            ELSE                    
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ssh_m(ji,jj) = zcoef * sshn(ji,jj)
+                  END DO
+               END DO
             ENDIF
             !
-            e3t_m(:,:) = zcoef * e3t_n(:,:,1)
-            !
-            frq_m(:,:) = zcoef * fraqsr_1lev(:,:)
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  e3t_m(ji,jj) = zcoef * e3t_n(ji,jj,1)
+                  !
+                  frq_m(ji,jj) = zcoef * fraqsr_1lev(ji,jj)
+               END DO
+            END DO
             !                                             ! ---------------------------------------- !
          ELSEIF( MOD( kt - 2 , nn_fsbc ) == 0 ) THEN      !   Initialisation: New mean computation   !
             !                                             ! ---------------------------------------- !
-            ssu_m(:,:) = 0._wp     ! reset to zero ocean mean sbc fields
-            ssv_m(:,:) = 0._wp
-            sst_m(:,:) = 0._wp
-            sss_m(:,:) = 0._wp
-            ssh_m(:,:) = 0._wp
-            e3t_m(:,:) = 0._wp
-            frq_m(:,:) = 0._wp
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ssu_m(ji,jj) = 0._wp     ! reset to zero ocean mean sbc fields
+                  ssv_m(ji,jj) = 0._wp
+                  sst_m(ji,jj) = 0._wp
+                  sss_m(ji,jj) = 0._wp
+                  ssh_m(ji,jj) = 0._wp
+                  e3t_m(ji,jj) = 0._wp
+                  frq_m(ji,jj) = 0._wp
+               END DO
+            END DO
          ENDIF
          !                                                ! ---------------------------------------- !
          !                                                !        Cumulate at each time step        !
          !                                                ! ---------------------------------------- !
-         ssu_m(:,:) = ssu_m(:,:) + ub(:,:,1)
-         ssv_m(:,:) = ssv_m(:,:) + vb(:,:,1)
-         IF( l_useCT )  THEN    ;   sst_m(:,:) = sst_m(:,:) + eos_pt_from_ct( zts(:,:,jp_tem), zts(:,:,jp_sal) )
-         ELSE                    ;   sst_m(:,:) = sst_m(:,:) + zts(:,:,jp_tem)
-         ENDIF
-         sss_m(:,:) = sss_m(:,:) + zts(:,:,jp_sal)
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ssu_m(ji,jj) = ssu_m(ji,jj) + ub(ji,jj,1)
+               ssv_m(ji,jj) = ssv_m(ji,jj) + vb(ji,jj,1)
+            END DO
+         END DO
+         IF( l_useCT )  THEN   
+           sst_m(:,:) = sst_m(:,:) + eos_pt_from_ct( zts(:,:,jp_tem), zts(:,:,jp_sal) )
+         ELSE                   
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+           DO jj = 1, jpj
+              DO ji = 1, jpi
+                 sst_m(ji,jj) = sst_m(ji,jj) + zts(ji,jj,jp_tem)
+              END DO
+           END DO
+         ENDIF
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               sss_m(ji,jj) = sss_m(ji,jj) + zts(ji,jj,jp_sal)
+            END DO
+         END DO
          !                          ! removed inverse barometer ssh when Patm forcing is used (for sea-ice dynamics)
-         IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = ssh_m(:,:) + sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )
-         ELSE                    ;   ssh_m(:,:) = ssh_m(:,:) + sshn(:,:)
-         ENDIF
-         !
-         e3t_m(:,:) = e3t_m(:,:) + e3t_n(:,:,1)
-         !
-         frq_m(:,:) = frq_m(:,:) + fraqsr_1lev(:,:)
+         IF( ln_apr_dyn ) THEN   
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ssh_m(ji,jj) = ssh_m(ji,jj) + sshn(ji,jj) - 0.5 * ( ssh_ib(ji,jj) + ssh_ibb(ji,jj) )
+               END DO
+            END DO
+         ELSE                    
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+           DO jj = 1, jpj
+              DO ji = 1, jpi
+                 ssh_m(ji,jj) = ssh_m(ji,jj) + sshn(ji,jj)
+              END DO
+           END DO
+         ENDIF
+         !
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               e3t_m(ji,jj) = e3t_m(ji,jj) + e3t_n(ji,jj,1)
+               !
+               frq_m(ji,jj) = frq_m(ji,jj) + fraqsr_1lev(ji,jj)
+            END DO
+         END DO
 
          !                                                ! ---------------------------------------- !
@@ -138,11 +246,16 @@
             !                                             ! ---------------------------------------- !
             zcoef = 1. / REAL( nn_fsbc, wp )
-            sst_m(:,:) = sst_m(:,:) * zcoef     ! mean SST             [Celsius]
-            sss_m(:,:) = sss_m(:,:) * zcoef     ! mean SSS             [psu]
-            ssu_m(:,:) = ssu_m(:,:) * zcoef     ! mean suface current  [m/s]
-            ssv_m(:,:) = ssv_m(:,:) * zcoef     !
-            ssh_m(:,:) = ssh_m(:,:) * zcoef     ! mean SSH             [m]
-            e3t_m(:,:) = e3t_m(:,:) * zcoef     ! mean vertical scale factor [m]
-            frq_m(:,:) = frq_m(:,:) * zcoef     ! mean fraction of solar net radiation absorbed in the 1st T level [-]
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  sst_m(ji,jj) = sst_m(ji,jj) * zcoef     ! mean SST             [Celsius]
+                  sss_m(ji,jj) = sss_m(ji,jj) * zcoef     ! mean SSS             [psu]
+                  ssu_m(ji,jj) = ssu_m(ji,jj) * zcoef     ! mean suface current  [m/s]
+                  ssv_m(ji,jj) = ssv_m(ji,jj) * zcoef     !
+                  ssh_m(ji,jj) = ssh_m(ji,jj) * zcoef     ! mean SSH             [m]
+                  e3t_m(ji,jj) = e3t_m(ji,jj) * zcoef     ! mean vertical scale factor [m]
+                  frq_m(ji,jj) = frq_m(ji,jj) * zcoef     ! mean fraction of solar net radiation absorbed in the 1st T level [-]
+               END DO
+            END DO
             !
          ENDIF
@@ -190,4 +303,5 @@
       !!----------------------------------------------------------------------
       REAL(wp) ::   zcoef, zf_sbc   ! local scalar
+      INTEGER  ::   ji, jj          ! loop index
       !!----------------------------------------------------------------------
       !
@@ -217,5 +331,10 @@
                CALL iom_get( numror, jpdom_autoglo, 'frq_m'  , frq_m  )
             ELSE
-               frq_m(:,:) = 1._wp   ! default definition
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     frq_m(ji,jj) = 1._wp   ! default definition
+                  END DO
+               END DO
             ENDIF
             !
@@ -223,11 +342,16 @@
                IF(lwp) WRITE(numout,*) '   restart with a change in the frequency of mean from ', zf_sbc, ' to ', nn_fsbc 
                zcoef = REAL( nn_fsbc - 1, wp ) / zf_sbc 
-               ssu_m(:,:) = zcoef * ssu_m(:,:) 
-               ssv_m(:,:) = zcoef * ssv_m(:,:)
-               sst_m(:,:) = zcoef * sst_m(:,:)
-               sss_m(:,:) = zcoef * sss_m(:,:)
-               ssh_m(:,:) = zcoef * ssh_m(:,:)
-               e3t_m(:,:) = zcoef * e3t_m(:,:)
-               frq_m(:,:) = zcoef * frq_m(:,:)
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ssu_m(ji,jj) = zcoef * ssu_m(ji,jj) 
+                     ssv_m(ji,jj) = zcoef * ssv_m(ji,jj)
+                     sst_m(ji,jj) = zcoef * sst_m(ji,jj)
+                     sss_m(ji,jj) = zcoef * sss_m(ji,jj)
+                     ssh_m(ji,jj) = zcoef * ssh_m(ji,jj)
+                     e3t_m(ji,jj) = zcoef * e3t_m(ji,jj)
+                     frq_m(ji,jj) = zcoef * frq_m(ji,jj)
+                  END DO
+               END DO
             ELSE
                IF(lwp) WRITE(numout,*) '   mean fields read in the ocean restart file'
@@ -239,13 +363,23 @@
          !
          IF(lwp) WRITE(numout,*) '   default initialisation of ss._m arrays'
-         ssu_m(:,:) = ub(:,:,1)
-         ssv_m(:,:) = vb(:,:,1)
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ssu_m(ji,jj) = ub(ji,jj,1)
+                  ssv_m(ji,jj) = vb(ji,jj,1)
+               END DO
+            END DO
          IF( l_useCT )  THEN    ;   sst_m(:,:) = eos_pt_from_ct( tsn(:,:,1,jp_tem), tsn(:,:,1,jp_sal) )
          ELSE                   ;   sst_m(:,:) = tsn(:,:,1,jp_tem)
          ENDIF
-         sss_m(:,:) = tsn  (:,:,1,jp_sal)
-         ssh_m(:,:) = sshn (:,:)
-         e3t_m(:,:) = e3t_n(:,:,1)
-         frq_m(:,:) = 1._wp
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               sss_m(ji,jj) = tsn  (ji,jj,1,jp_sal)
+               ssh_m(ji,jj) = sshn (ji,jj)
+               e3t_m(ji,jj) = e3t_n(ji,jj,1)
+               frq_m(ji,jj) = 1._wp
+            END DO
+         END DO
          !
       ENDIF
Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90	(revision 7646)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90	(revision 7698)
@@ -93,4 +93,5 @@
             !
             IF( nn_sstr == 1 ) THEN                                   !* Temperature restoring term
+!$OMP PARALLEL DO schedule(static) private(jj,ji,zqrp)
                DO jj = 1, jpj
                   DO ji = 1, jpi
@@ -105,4 +106,5 @@
             IF( nn_sssr == 1 ) THEN                                   !* Salinity damping term (salt flux only (sfx))
                zsrp = rn_deds / rday                                  ! from [mm/day] to [kg/m2/s]
+!$OMP PARALLEL DO schedule(static) private(jj, ji, zerp)
                DO jj = 1, jpj
                   DO ji = 1, jpi
@@ -118,4 +120,5 @@
                zsrp = rn_deds / rday                                  ! from [mm/day] to [kg/m2/s]
                zerp_bnd = rn_sssr_bnd / rday                          !       -              -    
+!$OMP PARALLEL DO schedule(static) private(jj, ji, zerp)
                DO jj = 1, jpj
                   DO ji = 1, jpi                            
