Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90	(revision 7698)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90	(revision 7753)
@@ -316,30 +316,14 @@
 #if defined key_cice
       IF( MOD( kt - 1, nn_fsbc ) == 0 )   THEN
-!$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
+         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) 
          ENDIF 
-!$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
+         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)
       ENDIF
 #endif
@@ -398,23 +382,17 @@
       !
 
-!$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
+      ! 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)
+
       ! ----------------------------------------------------------------------------- !
       !      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.
@@ -424,5 +402,4 @@
       END DO
 #endif
-!$OMP PARALLEL DO schedule(static) private(jj, ji)
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vect. opt.
@@ -434,12 +411,7 @@
       CALL lbc_lnk( zwnd_j(:,:) , 'T', -1. )
       ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked)
-!$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
+      wndm(:,:) = SQRT(  zwnd_i(:,:) * zwnd_i(:,:)   &
+         &             + zwnd_j(:,:) * zwnd_j(:,:)  ) * tmask(:,:,1)
+
       ! ----------------------------------------------------------------------------- !
       !      I   Radiative FLUXES                                                     !
@@ -449,19 +421,8 @@
       zztmp = 1. - albo
       IF( ln_dm2dc ) THEN   ;   qsr(:,:) = zztmp * sbc_dcy( 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
-
-!$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
+      ELSE                  ;   qsr(:,:) = zztmp *          sf(jp_qsr)%fnow(:,:,1)   * tmask(:,:,1)
+      ENDIF
+
+      zqlw(:,:) = (  sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:)  ) * tmask(:,:,1)   ! Long  Wave
 
 
@@ -500,13 +461,6 @@
       END IF
 
-!$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)
+      Cd_oce(:,:) = Cd(:,:)  ! record value of pure ocean-atm. drag (clem)
+
       DO jj = 1, jpj             ! tau module, i and j component
          DO ji = 1, jpi
@@ -517,15 +471,7 @@
          END DO
       END DO
-!$OMP END PARALLEL
 
       !                          ! add the HF tau contribution to the wind stress module
-      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
+      IF( lhftau )   taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1)
 
       CALL iom_put( "taum_oce", taum )   ! output wind stress module
@@ -534,5 +480,4 @@
       !     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
@@ -551,29 +496,14 @@
 
       ! zqla used as temporary array, for rho*U (common term of bulk formulae):
-!$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
+      zqla(:,:) = zrhoa(:,:) * zU_zu(:,:)
 
       IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN
          !! q_air and t_air are given at 10m (wind reference height)
-!$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
+         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
       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!!!
-!$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
+         zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce(:,:)*(zsq(:,:) - zq_zu(:,:) ) ) ! Evaporation ! using bulk wind speed
          zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch(:,:)*(zst(:,:) - zt_zu(:,:) )   ! Sensible Heat ! using bulk wind speed
       ENDIF
@@ -597,24 +527,19 @@
       ! ----------------------------------------------------------------------------- !
       !
-!$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)
-            !
+      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)
+      !
 #if defined key_lim3
-            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)
+      qns_oce(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                ! non solar without emp (only needed by LIM3)
+      qsr_oce(:,:) = qsr(:,:)
 #endif
-         END DO
-      END DO
       !
       IF ( nn_ice == 0 ) THEN
@@ -626,11 +551,6 @@
          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
-!$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
+         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]
          CALL iom_put( 'snowpre', sprecip * 86400. )        ! Snow
          CALL iom_put( 'precip' , tprecip * 86400. )        ! Total precipitation
@@ -679,10 +599,5 @@
       CALL wrk_alloc( jpi,jpj, Cd )
 
-!$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
+      Cd(:,:) = Cd_ice
 
       ! Make ice-atm. drag dependent on ice concentration (see Lupkes et al. 2012) (clem)
@@ -698,14 +613,9 @@
       zrhoa (:,:) =  rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1))
 
-!$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
+      !!gm brutal....
+      utau_ice  (:,:) = 0._wp
+      vtau_ice  (:,:) = 0._wp
+      wndm_ice  (:,:) = 0._wp
+      !!gm end
 
       ! ----------------------------------------------------------------------------- !
@@ -715,5 +625,4 @@
       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
@@ -740,5 +649,4 @@
          !
       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.
@@ -748,5 +656,4 @@
             END DO
          END DO
-!$OMP PARALLEL DO schedule(static) private(jj,ji)
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vect. opt.
@@ -793,10 +700,8 @@
       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)
@@ -805,14 +710,9 @@
       IF( nn_timing == 1 )  CALL timing_start('blk_ice_flx')
       !
-      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,jpl,   z_qlw, z_qsb, z_dqlw, z_dqsb )
+      CALL wrk_alloc( jpi,jpj,       zrhoa)
       CALL wrk_alloc( jpi,jpj, Cd )
 
-!$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
+      Cd(:,:) = Cd_ice
 
       ! Make ice-atm. drag dependent on ice concentration (see Lupkes et al.  2012) (clem)
@@ -831,8 +731,7 @@
       !
       zztmp = 1. / ( 1. - albo )
-!$OMP PARALLEL
-!$OMP DO schedule(static) private(jl,jj,ji,zst2,zst3)            ! ========================== !
-      DO jl = 1, jpl                                             !  Loop over ice categories  !
-         !                                                       ! ========================== !
+      !                                     ! ========================== !
+      DO jl = 1, jpl                        !  Loop over ice categories  !
+         !                                  ! ========================== !
          DO jj = 1 , jpj
             DO ji = 1, jpi
@@ -882,12 +781,6 @@
       END DO
       !
-!$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
+      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]
       CALL iom_put( 'snowpre', sprecip * 86400. )                  ! Snow precipitation
       CALL iom_put( 'precip' , tprecip * 86400. )                  ! Total precipitation
@@ -898,98 +791,35 @@
       ! --- evaporation --- !
       z1_lsub = 1._wp / Lsub
-!$OMP PARALLEL
-!$OMP DO schedule(static) private(jl,jj,ji)
+      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
+      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 )
+
+      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) ---
       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
+         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
-      !
-!$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
-!$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
-         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 )
@@ -1001,11 +831,6 @@
       ! ( Maykut and Untersteiner, 1971 ; Ebert and Curry, 1993 )
       !
-!$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
+      fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
+      fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
       !
       !
@@ -1019,7 +844,7 @@
       ENDIF
 
-      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,jpl,   z_qlw, z_qsb, z_dqlw, z_dqsb )
       CALL wrk_dealloc( jpi,jpj,       zrhoa )
-      CALL wrk_dealloc( jpi,jpj, Cd, zevap_ice2d, zqns_ice2d, zqsr_ice2d)
+      CALL wrk_dealloc( jpi,jpj, Cd )
       !
       IF( nn_timing == 1 )  CALL timing_stop('blk_ice_flx')
@@ -1083,5 +908,4 @@
       !!----------------------------------------------------------------------------------
       !
-!$OMP PARALLEL DO schedule(static) private(jj,ji,ztmp,ze_sat)
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -1120,5 +944,4 @@
       !!----------------------------------------------------------------------------------
       !
-!$OMP PARALLEL DO schedule(static) private(jj,ji,zrv,ziRT)
       DO jj = 1, jpj
          DO ji = 1, jpi
