Index: /NEMO/trunk/src/ICE/icedyn_adv_umx.F90
===================================================================
--- /NEMO/trunk/src/ICE/icedyn_adv_umx.F90	(revision 10944)
+++ /NEMO/trunk/src/ICE/icedyn_adv_umx.F90	(revision 10945)
@@ -32,28 +32,21 @@
 
    PUBLIC   ice_dyn_adv_umx   ! called by icedyn_adv.F90
-      
-   REAL(wp) ::   z1_6   = 1._wp /   6._wp   ! =1/6
-   REAL(wp) ::   z1_120 = 1._wp / 120._wp   ! =1/120
-   
-   ! advection for S and T:    dVS/dt = -div( uA * uHS / u ) => ll_advS = F
-   !                        or dVS/dt = -div( uV * uS  / u ) => ll_advS = T
-   LOGICAL ::   ll_advS = .FALSE.
    !
-   ! alternate directions for upstream
-   LOGICAL ::   ll_upsxy = .TRUE.
+   INTEGER, PARAMETER ::   np_advS = 1         ! advection for S and T:    dVS/dt = -div(      uVS     ) => np_advS = 1
+   !                                                                    or dVS/dt = -div( uA * uHS / u ) => np_advS = 2
+   !                                                                    or dVS/dt = -div( uV * uS  / u ) => np_advS = 3
+   INTEGER, PARAMETER ::   np_limiter = 1      ! limiter: 1 = nonosc
+   !                                                      2 = superbee
+   !                                                      3 = h3
+   LOGICAL            ::   ll_upsxy  = .TRUE.   ! alternate directions for upstream
+   LOGICAL            ::   ll_hoxy   = .TRUE.   ! alternate directions for high order
+   LOGICAL            ::   ll_neg    = .TRUE.   ! if T interpolated at u/v points is negative or v_i < 1.e-6
+   !                                                 then interpolate T at u/v points using the upstream scheme
+   LOGICAL            ::   ll_prelim = .FALSE.  ! prelimiter from: Zalesak(1979) eq. 14 => not well defined in 2D
    !
-   ! alternate directions for high order
-   LOGICAL ::   ll_hoxy = .TRUE.
+   REAL(wp)           ::   z1_6   = 1._wp /   6._wp   ! =1/6
+   REAL(wp)           ::   z1_120 = 1._wp / 120._wp   ! =1/120
    !
-   ! if T interpolated at u/v points is negative or v_i < 1.e-6
-   !    then interpolate T at u/v points using the upstream scheme
-   LOGICAL ::   ll_neg = .TRUE.
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zmski_small, zmskj_small
-   !
-   ! limiter: 1=nonosc_ice, 2=superbee, 3=h3(rachid)
-   INTEGER ::   kn_limiter = 1
-   !
-   ! prelimiter: use it to avoid overshoot in H
-   LOGICAL ::   ll_prelimiter_zalesak = .FALSE.  ! from: Zalesak(1979) eq. 14 => better for 1D. Not well defined in 2D
+   INTEGER, ALLOCATABLE, DIMENSION(:,:,:) ::   imsk_small, jmsk_small
    !
    !! * Substitutions
@@ -100,11 +93,12 @@
       REAL(wp) ::   zdt, zvi_cen
       REAL(wp), DIMENSION(1)           ::   zcflprv, zcflnow   ! for global communication
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zudy, zvdx, zcu_box, zcv_box 
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zudy, zvdx, zcu_box, zcv_box
       REAL(wp), DIMENSION(jpi,jpj)     ::   zati1, zati2
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zu_cat, zv_cat
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zua_ho, zva_ho, zua_ups, zva_ups
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zuv_ho, zvv_ho, zuv_ups, zvv_ups 
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_ai, z1_aip, z1_vi, z1_vs 
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zhvar
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_ai , z1_aip, zhvar
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zhi_max, zhs_max, zhip_max
+      !
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs 
       !!----------------------------------------------------------------------
       !
@@ -150,5 +144,11 @@
       zudy(:,:) = pu_ice(:,:) * e2u(:,:)
       zvdx(:,:) = pv_ice(:,:) * e1v(:,:)
-
+      !
+      ! setup transport for each ice cat 
+      DO jl = 1, jpl
+         zu_cat(:,:,jl) = zudy(:,:)
+         zv_cat(:,:,jl) = zvdx(:,:)
+      END DO
+      !
       ! --- define velocity for advection: u*grad(H) --- !
       DO jj = 2, jpjm1
@@ -184,15 +184,15 @@
          ! setup a mask where advection will be upstream
          IF( ll_neg ) THEN
-            IF( .NOT. ALLOCATED(zmski_small) )   ALLOCATE( zmski_small(jpi,jpj,jpl) ) 
-            IF( .NOT. ALLOCATED(zmskj_small) )   ALLOCATE( zmskj_small(jpi,jpj,jpl) ) 
+            IF( .NOT. ALLOCATED(imsk_small) )   ALLOCATE( imsk_small(jpi,jpj,jpl) ) 
+            IF( .NOT. ALLOCATED(jmsk_small) )   ALLOCATE( jmsk_small(jpi,jpj,jpl) ) 
             DO jl = 1, jpl
                DO jj = 1, jpjm1
                   DO ji = 1, jpim1
                      zvi_cen = 0.5_wp * ( pv_i(ji+1,jj,jl) + pv_i(ji,jj,jl) )
-                     IF( zvi_cen < epsi06) THEN ; zmski_small(ji,jj,jl) = 0.
-                     ELSE                       ; zmski_small(ji,jj,jl) = 1. ; ENDIF
+                     IF( zvi_cen < epsi06) THEN   ;   imsk_small(ji,jj,jl) = 0
+                     ELSE                         ;   imsk_small(ji,jj,jl) = 1   ;   ENDIF
                      zvi_cen = 0.5_wp * ( pv_i(ji,jj+1,jl) + pv_i(ji,jj,jl) )
-                     IF( zvi_cen < epsi06) THEN ; zmskj_small(ji,jj,jl) = 0.
-                     ELSE                       ; zmskj_small(ji,jj,jl) = 1. ; ENDIF
+                     IF( zvi_cen < epsi06) THEN   ;   jmsk_small(ji,jj,jl) = 0
+                     ELSE                         ;   jmsk_small(ji,jj,jl) = 1   ;   ENDIF
                   END DO
                END DO
@@ -204,22 +204,17 @@
          ! ----------------------- !
          !
-         ! set u*A=u for advection of A only 
-         DO jl = 1, jpl
-            zua_ho(:,:,jl) = zudy(:,:)
-            zva_ho(:,:,jl) = zvdx(:,:)
-         END DO
-
          !== Ice area ==!
          zamsk = 1._wp
-         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
+         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zu_cat , zv_cat , zcu_box, zcv_box, &
             &                                      pa_i, pa_i, zua_ups, zva_ups, zua_ho , zva_ho )
-         zamsk = 0._wp
-
-         IF( .NOT. ll_advS ) THEN   !-- advection form: uA * uHS / u
+         !
+         !                             ! --------------------------------- !
+         IF( np_advS == 1 ) THEN       ! -- advection form: -div( uVS ) -- !
+            !                          ! --------------------------------- !
+            zamsk = 0._wp
             !== Ice volume ==!
             zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:)
             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
                &                                      zhvar, pv_i, zua_ups, zva_ups )
-            !
             !== Snw volume ==!         
             zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:)
@@ -227,9 +222,35 @@
                &                                      zhvar, pv_s, zua_ups, zva_ups )
             !
+            zamsk = 1._wp
+            !== Salt content ==!
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, &
+               &                                      psv_i, psv_i )
+            !== Ice heat content ==!
+            DO jk = 1, nlay_i
+               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, &
+                  &                                      pe_i(:,:,jk,:), pe_i(:,:,jk,:) )
+            END DO
+            !== Snw heat content ==!
+            DO jk = 1, nlay_s
+               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, &
+                  &                                      pe_s(:,:,jk,:), pe_s(:,:,jk,:) )
+            END DO
+            !
+            !                          ! ------------------------------------------ !
+         ELSEIF( np_advS == 2 ) THEN   ! -- advection form: -div( uA * uHS / u ) -- !
+            !                          ! ------------------------------------------ !
+            zamsk = 0._wp
+            !== Ice volume ==!
+            zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
+               &                                      zhvar, pv_i, zua_ups, zva_ups )
+            !== Snw volume ==!         
+            zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
+               &                                      zhvar, pv_s, zua_ups, zva_ups )
             !== Salt content ==!
             zhvar(:,:,:) = psv_i(:,:,:) * z1_ai(:,:,:)
             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, &
                &                                      zhvar, psv_i, zua_ups, zva_ups )
-            !
             !== Ice heat content ==!
             DO jk = 1, nlay_i
@@ -238,5 +259,4 @@
                   &                                      zhvar, pe_i(:,:,jk,:), zua_ups, zva_ups )
             END DO
-            !
             !== Snw heat content ==!
             DO jk = 1, nlay_s
@@ -246,5 +266,12 @@
             END DO
             !
-         ELSE                       !-- advection form: uV * uS / u
+            !                          ! ----------------------------------------- !
+         ELSEIF( np_advS == 3 ) THEN   ! -- advection form: -div( uV * uS / u ) -- !
+            !                          ! ----------------------------------------- !
+            zamsk = 0._wp
+            !
+            ALLOCATE( zuv_ho (jpi,jpj,jpl), zvv_ho (jpi,jpj,jpl),  &
+               &      zuv_ups(jpi,jpj,jpl), zvv_ups(jpi,jpj,jpl), z1_vi(jpi,jpj,jpl), z1_vs(jpi,jpj,jpl) )
+            !
             ! inverse of Vi
             WHERE( pv_i(:,:,:) >= epsi20 )   ;   z1_vi(:,:,:) = 1._wp / pv_i(:,:,:)
@@ -264,10 +291,8 @@
             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
                &                                      zhvar, pv_i, zuv_ups, zvv_ups, zuv_ho , zvv_ho )
-            !
             !== Salt content ==!
             zhvar(:,:,:) = psv_i(:,:,:) * z1_vi(:,:,:)
             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zuv_ho , zvv_ho , zcu_box, zcv_box, &
                &                                      zhvar, psv_i, zuv_ups, zvv_ups )
-            !
             !== Ice heat content ==!
             DO jk = 1, nlay_i
@@ -276,5 +301,4 @@
                   &                                      zhvar, pe_i(:,:,jk,:), zuv_ups, zvv_ups )
             END DO
-            !
             !== Snow volume ==!         
             zuv_ups = zua_ups
@@ -283,5 +307,4 @@
             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
                &                                      zhvar, pv_s, zuv_ups, zvv_ups, zuv_ho , zvv_ho )
-            !
             !== Snw heat content ==!
             DO jk = 1, nlay_s
@@ -291,37 +314,23 @@
             END DO
             !
+            DEALLOCATE( zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs )
+            !
          ENDIF
          !
          !== Ice age ==!
          IF( iom_use('iceage') .OR. iom_use('iceage_cat') ) THEN
-            ! clem: in theory we should use the formulation below to advect the ice age, but the code is unable to deal with
-            !       fields that do not depend on volume (here oa_i depends on concentration). It creates abnormal ages that
-            !       spread into the domain. Therefore we cheat and consider that ice age should be advected as ice concentration
-            !!zhvar(:,:,:) = poa_i(:,:,:) * z1_ai(:,:,:)
-            !!CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, poa_i )
-            ! set u*A=u for advection of ice age
-            DO jl = 1, jpl
-               zua_ho(:,:,jl) = zudy(:,:)
-               zva_ho(:,:,jl) = zvdx(:,:)
-            END DO
             zamsk = 1._wp
-            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho, zva_ho, zcu_box, zcv_box, &
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, &
                &                                      poa_i, poa_i )
-            zamsk = 0._wp
          ENDIF
          !
          !== melt ponds ==!
          IF ( ln_pnd_H12 ) THEN
-            ! set u*A=u for advection of Ap only 
-            DO jl = 1, jpl
-               zua_ho(:,:,jl) = zudy(:,:)
-               zva_ho(:,:,jl) = zvdx(:,:)
-            END DO
             ! fraction
             zamsk = 1._wp
-            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, &
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat , zv_cat , zcu_box, zcv_box, &
                &                                      pa_ip, pa_ip, zua_ups, zva_ups, zua_ho , zva_ho )
+            ! volume
             zamsk = 0._wp
-            ! volume
             zhvar(:,:,:) = pv_ip(:,:,:) * z1_aip(:,:,:)
             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, &
@@ -341,5 +350,4 @@
          !
          ! --- Ensure non-negative fields and in-bound thicknesses --- !
-         !
          ! Remove negative values (conservation is ensured)
          !    (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20)
@@ -474,7 +482,7 @@
          CALL lbc_lnk( 'icedyn_adv_umx', zt_ups, 'T',  1. )
          !
-         IF    ( kn_limiter == 1 ) THEN
+         IF    ( np_limiter == 1 ) THEN
             CALL nonosc_ice( 1._wp, pdt, pu, pv, ptc, zt_ups, zfu_ups, zfv_ups, zfu_ho, zfv_ho )
-         ELSEIF( kn_limiter == 2 .OR. kn_limiter == 3 ) THEN
+         ELSEIF( np_limiter == 2 .OR. np_limiter == 3 ) THEN
             CALL limiter_x( pdt, pu, ptc, zfu_ups, zfu_ho )
             CALL limiter_y( pdt, pv, ptc, zfv_ups, zfv_ho )
@@ -659,7 +667,7 @@
          END DO
          !
-         IF    ( kn_limiter == 1 ) THEN
+         IF    ( np_limiter == 1 ) THEN
             CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
-         ELSEIF( kn_limiter == 2 .OR. kn_limiter == 3 ) THEN
+         ELSEIF( np_limiter == 2 .OR. np_limiter == 3 ) THEN
             CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
             CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
@@ -677,5 +685,5 @@
                END DO
             END DO
-            IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
+            IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
 
             DO jl = 1, jpl              !-- first guess of tracer from u-flux
@@ -698,5 +706,5 @@
                END DO
             END DO
-            IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
+            IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
 
          ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==!
@@ -709,5 +717,5 @@
                END DO
             END DO
-            IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
+            IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
             !
             DO jl = 1, jpl              !-- first guess of tracer from v-flux
@@ -730,8 +738,8 @@
                END DO
             END DO
-            IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
+            IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
 
          ENDIF
-         IF( kn_limiter == 1 )   CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
+         IF( np_limiter == 1 )   CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
          
       ENDIF
@@ -771,5 +779,5 @@
          CALL ultimate_x( pamsk, kn_umx, pdt, pt, pu, zt_u, pfu_ho )
          !                                                        !--  limiter in x --!
-         IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
+         IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
          !                                                        !--  advective form update in zpt  --!
          DO jl = 1, jpl
@@ -792,5 +800,5 @@
          ENDIF
          !                                                        !--  limiter in y --!
-         IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
+         IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
          !         
          !
@@ -800,5 +808,5 @@
          CALL ultimate_y( pamsk, kn_umx, pdt, pt, pv, zt_v, pfv_ho )
          !                                                        !--  limiter in y --!
-         IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
+         IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
          !                                                        !--  advective form update in zpt  --!
          DO jl = 1, jpl
@@ -821,9 +829,9 @@
          ENDIF
          !                                                        !--  limiter in x --!
-         IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
+         IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
          !
       ENDIF
 
-      IF( kn_limiter == 1 )   CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
+      IF( np_limiter == 1 )   CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
       !
    END SUBROUTINE macho
@@ -967,5 +975,5 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1
-                  IF( pt_u(ji,jj,jl) < 0._wp .OR. ( zmski_small(ji,jj,jl) == 0. .AND. pamsk == 0. ) ) THEN
+                  IF( pt_u(ji,jj,jl) < 0._wp .OR. ( imsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN
                      pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   &
                         &                                         - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
@@ -1121,5 +1129,5 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1
-                  IF( pt_v(ji,jj,jl) < 0._wp .OR. ( zmskj_small(ji,jj,jl) == 0. .AND. pamsk == 0. ) ) THEN
+                  IF( pt_v(ji,jj,jl) < 0._wp .OR. ( jmsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN
                      pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  &
                         &                                         - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
@@ -1185,5 +1193,5 @@
       !                        |      |      |        |    *
       !            t_ups :       i-1     i       i+1       i+2   
-      IF( ll_prelimiter_zalesak ) THEN
+      IF( ll_prelim ) THEN
          
          DO jl = 1, jpl
@@ -1352,5 +1360,5 @@
                Rjp = zslpx(ji+1,jj,jl)
 
-               IF( kn_limiter == 3 ) THEN
+               IF( np_limiter == 3 ) THEN
 
                   IF( pu(ji,jj) > 0. ) THEN   ;   Rr = Rjm
@@ -1368,5 +1376,5 @@
                   pfu_ho(ji,jj,jl) = pfu_ups(ji,jj,jl) + zlimiter
 
-               ELSEIF( kn_limiter == 2 ) THEN
+               ELSEIF( np_limiter == 2 ) THEN
                   IF( Rj /= 0. ) THEN
                      IF( pu(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj
@@ -1447,5 +1455,5 @@
                Rjp = zslpy(ji,jj+1,jl)
 
-               IF( kn_limiter == 3 ) THEN
+               IF( np_limiter == 3 ) THEN
 
                   IF( pv(ji,jj) > 0. ) THEN   ;   Rr = Rjm
@@ -1463,5 +1471,5 @@
                   pfv_ho(ji,jj,jl) = pfv_ups(ji,jj,jl) + zlimiter
 
-               ELSEIF( kn_limiter == 2 ) THEN
+               ELSEIF( np_limiter == 2 ) THEN
 
                   IF( Rj /= 0. ) THEN
Index: /NEMO/trunk/src/ICE/icedyn_rdgrft.F90
===================================================================
--- /NEMO/trunk/src/ICE/icedyn_rdgrft.F90	(revision 10944)
+++ /NEMO/trunk/src/ICE/icedyn_rdgrft.F90	(revision 10945)
@@ -945,5 +945,19 @@
          CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' )
       ENDIF
-      !                              ! allocate tke arrays
+      !
+      IF( .NOT. ln_icethd ) THEN
+         rn_porordg = 0._wp
+         rn_fsnwrdg = 1._wp ; rn_fsnwrft = 1._wp
+         rn_fpndrdg = 1._wp ; rn_fpndrft = 1._wp
+         IF( lwp ) THEN
+            WRITE(numout,*) '      ==> only ice dynamics is activated, thus some parameters must be changed'
+            WRITE(numout,*) '            rn_porordg   = ', rn_porordg
+            WRITE(numout,*) '            rn_fsnwrdg   = ', rn_fsnwrdg 
+            WRITE(numout,*) '            rn_fpndrdg   = ', rn_fpndrdg 
+            WRITE(numout,*) '            rn_fsnwrft   = ', rn_fsnwrft 
+            WRITE(numout,*) '            rn_fpndrft   = ', rn_fpndrft 
+         ENDIF
+      ENDIF
+      !                              ! allocate arrays
       IF( ice_dyn_rdgrft_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'ice_dyn_rdgrft_init: unable to allocate arrays' )
       !
Index: /NEMO/trunk/src/ICE/icevar.F90
===================================================================
--- /NEMO/trunk/src/ICE/icevar.F90	(revision 10944)
+++ /NEMO/trunk/src/ICE/icevar.F90	(revision 10945)
@@ -563,4 +563,7 @@
       DO jl = 1, jpl       !==  loop over the categories  ==!
          !
+         ! make sure a_i=0 where v_i<=0
+         WHERE( pv_i(:,:,:) <= 0._wp )   pa_i(:,:,:) = 0._wp
+
          !----------------------------------------
          ! zap ice energy and send it to the ocean
@@ -569,5 +572,5 @@
             DO jj = 1 , jpj
                DO ji = 1 , jpi
-                  IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN
+                  IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
                      hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * z1_dt ! W.m-2 >0
                      pe_i(ji,jj,jk,jl) = 0._wp
@@ -580,5 +583,5 @@
             DO jj = 1 , jpj
                DO ji = 1 , jpi
-                  IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN
+                  IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
                      hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * z1_dt ! W.m-2 <0
                      pe_s(ji,jj,jk,jl) = 0._wp
@@ -593,13 +596,13 @@
          DO jj = 1 , jpj
             DO ji = 1 , jpi
-               IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN
+               IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
                   wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * z1_dt
                   pv_i   (ji,jj,jl) = 0._wp
                ENDIF
-               IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN
+               IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
                   wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * z1_dt
                   pv_s   (ji,jj,jl) = 0._wp
                ENDIF
-               IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN
+               IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
                   sfx_res(ji,jj)    = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * z1_dt
                   psv_i  (ji,jj,jl) = 0._wp
