Index: /trunk/NEMO/LIM_SRC/iceini.F90
===================================================================
--- /trunk/NEMO/LIM_SRC/iceini.F90	(revision 75)
+++ /trunk/NEMO/LIM_SRC/iceini.F90	(revision 76)
@@ -29,5 +29,5 @@
    !! * Share Module variables
    LOGICAL , PUBLIC  ::   & !!! ** init namelist (namicerun) **
-      ldyn   = .TRUE.        !: flag for ice dynamics (T) or not (F)
+      ln_limdyn   = .TRUE.   !: flag for ice dynamics (T) or not (F)
    INTEGER , PUBLIC  ::   &  !:
       nstart ,            &  !: iteration number of the begining of the run 
@@ -36,9 +36,9 @@
       numit                  !: iteration number
    REAL(wp), PUBLIC  ::   &  !:
-      hsndif = 0.0 ,      &  !: computation of temp. in snow (0) or not (9999)
-      hicdif = 0.0 ,      &  !: computation of temp. in ice (0) or not (9999)
+      hsndif = 0.e0 ,     &  !: computation of temp. in snow (0) or not (9999)
+      hicdif = 0.e0 ,     &  !: computation of temp. in ice (0) or not (9999)
       tpstot                 !: time of the run in seconds
    REAL(wp), PUBLIC, DIMENSION(2)  ::  &  !:
-      acrit  = (/ 1.0e-06 , 1.0e-06 /)    !: minimum fraction for leads in 
+      acrit  = (/ 1.e-06 , 1.e-06 /)    !: minimum fraction for leads in 
       !                                   !  north and south hemisphere
    !!----------------------------------------------------------------------
@@ -111,5 +111,5 @@
       !!   2.0  !  03-08 (C. Ethe)  Original code
       !!-------------------------------------------------------------------
-      NAMELIST/namicerun/ ldyn, acrit, hsndif, hicdif
+      NAMELIST/namicerun/ ln_limdyn, acrit, hsndif, hicdif
       !!-------------------------------------------------------------------
 
@@ -121,6 +121,6 @@
          WRITE(numout,*) 'ice_run : ice share parameters for dynamics/advection/thermo of sea-ice'
          WRITE(numout,*) ' ~~~~~~'
-         WRITE(numout,*) '   switch for ice dynamics (1) or not (0)           ldyn   = ', ldyn
-         WRITE(numout,*) '   minimum fraction for leads in the NH (SH)  acrit(1/2)   = ', acrit
+         WRITE(numout,*) '   switch for ice dynamics (1) or not (0)      ln_limdyn   = ', ln_limdyn
+         WRITE(numout,*) '   minimum fraction for leads in the NH (SH)  acrit(1/2)   = ', acrit(:)
          WRITE(numout,*) '   computation of temp. in snow (=0) or not (=9999) hsndif = ', hsndif
          WRITE(numout,*) '   computation of temp. in ice  (=0) or not (=9999) hicdif = ', hicdif
Index: /trunk/NEMO/LIM_SRC/limdyn.F90
===================================================================
--- /trunk/NEMO/LIM_SRC/limdyn.F90	(revision 75)
+++ /trunk/NEMO/LIM_SRC/limdyn.F90	(revision 76)
@@ -20,6 +20,7 @@
    USE iceini
    USE limistate
-   USE limrhg     ! ice rheology
+   USE limrhg          ! ice rheology
    USE lbclnk
+   USE lib_mpp
 
    IMPLICIT NONE
@@ -30,5 +31,5 @@
 
    !! * Module variables
-   REAL(wp)  ::  rone    = 1.0   ! constant value
+   REAL(wp)  ::  rone    = 1.e0   ! constant value
 
    !!----------------------------------------------------------------------
@@ -55,20 +56,21 @@
       !!---------------------------------------------------------------------
       !! * Loal variables
-      INTEGER ::   ji, jj,        &  ! dummy loop indices
-         jhemis                      ! jhemis = 1 (NH) ; jhemis = -1 (SH)
+      INTEGER ::   ji, jj             ! dummy loop indices
+      INTEGER ::   i_j1, i_jpj        ! Starting/ending j-indices for rheology
       REAL(wp) ::   &
-         ztairx, ztairy,          &  ! tempory scalars
-         zsang , zmod,            &
-         ztglx , ztgly ,          &
-         zt11, zt12, zt21, zt22 , &
-         zustm, zsfrld, zsfrldm4, &
+         ztairx, ztairy,           &  ! tempory scalars
+         zsang , zmod,             &
+         ztglx , ztgly ,           &
+         zt11, zt12, zt21, zt22 ,  &
+         zustm, zsfrld, zsfrldm4,  &
          zu_ice, zv_ice, ztair2
-
+      REAL(wp),DIMENSION(jpj) ::   &
+         zind,                     &  ! i-averaged indicator of sea-ice
+         zmsk                         ! i-averaged of tmask
       !!---------------------------------------------------------------------
-
 
       IF( numit == nstart  )   CALL lim_dyn_init   ! Initialization (first time-step only)
       
-      IF ( ldyn ) THEN
+      IF ( ln_limdyn ) THEN
 
          ! Mean ice and snow thicknesses.          
@@ -80,16 +82,75 @@
        
          !                                         ! Rheology (ice dynamics)
-         !--  Northern hemisphere                  ! ========
-         jhemis = +1
-         CALL lim_rhg( jhemis )
-
-         !--  Southern hemisphere.
-         jhemis = -1
-         CALL lim_rhg( jhemis )
-
-         u_ice(:,1) = 0.0       !ibug  est-ce vraiment necessaire?
-         v_ice(:,1) = 0.0
-
-         IF( l_ctl .AND. lwp ) THEN 
+         !                                         ! ========
+         
+         !  Define the j-limits where ice rheology is computed
+         ! ---------------------------------------------------
+         
+         IF( lk_mpp ) THEN                    ! mpp: compute over the whole domain
+            i_j1 = 1   
+            i_jpj = jpj
+            IF(l_ctl)   WRITE(numout,*) 'lim_rhg :    i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
+            CALL lim_rhg( i_j1, i_jpj )
+
+         ELSE                                 ! optimization of the computational area
+
+            DO jj = 1, jpj
+               zind(jj) = SUM( frld (:,jj  ) )   ! = FLOAT(jpj) if ocean everywhere on a j-line
+               zmsk(jj) = SUM( tmask(:,jj,1) )   ! = 0          if land  everywhere on a j-line
+   !!i         write(numout,*) narea, 'limrhg' , jj, zind(jj), zmsk(jj)
+            END DO
+
+            IF( l_jeq ) THEN                     ! local domain include both hemisphere
+               !                                 ! Rheology is computed in each hemisphere
+               !                                 ! only over the ice cover latitude strip
+               ! Northern hemisphere
+               i_j1  = njeq
+               i_jpj = jpj
+               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
+                  i_j1 = i_j1 + 1
+               END DO
+               i_j1 = MAX( 1, i_j1-1 )
+               IF(l_ctl)   WRITE(numout,*) 'lim_rhg : NH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
+    
+               CALL lim_rhg( i_j1, i_jpj )
+    
+               ! Southern hemisphere
+               i_j1  =  1 
+               i_jpj = njeq
+               DO WHILE ( i_jpj >= 1 .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
+                  i_jpj = i_jpj - 1
+               END DO
+               i_jpj = MIN( jpj, i_jpj+2 )
+               IF(l_ctl)   WRITE(numout,*) 'lim_rhg : SH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
+    
+               CALL lim_rhg( i_j1, i_jpj )
+    
+            ELSE                                 ! local domain extends over one hemisphere only
+               !                                 ! Rheology is computed only over the ice cover
+               !                                 ! latitude strip
+               i_j1  = 1
+               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
+                  i_j1 = i_j1 + 1
+               END DO
+               i_j1 = MAX( 1, i_j1-1 )
+    
+               i_jpj  = jpj
+               DO WHILE ( i_jpj >= 1  .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
+                  i_jpj = i_jpj - 1
+               END DO
+               i_jpj = MIN( jpj, i_jpj+2)
+    
+               IF(l_ctl)   WRITE(numout,*) 'lim_rhg : one hemisphere: i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
+    
+               CALL lim_rhg( i_j1, i_jpj )
+
+            ENDIF
+
+         ENDIF
+
+         u_ice(:,1) = 0.e0       !ibug  est-ce vraiment necessaire?
+         v_ice(:,1) = 0.e0
+
+         IF(l_ctl) THEN 
             WRITE(numout,*) ' lim_dyn  : u_oce ', SUM( u_oce ), ' v_oce ', SUM( v_oce )
             WRITE(numout,*) ' lim_dyn  : u_ice ', SUM( u_ice ), ' v_ice ', SUM( v_ice )
@@ -99,6 +160,5 @@
          !                                         ! ================
          DO jj = 2, jpjm1
-            jhemis = SIGN(1, jj - jeq )
-            zsang  = REAL(jhemis) * sangvg
+            zsang  = SIGN(1.e0, gphif(1,jj-1) ) * sangvg
             DO ji = 2, jpim1
                ! computation of wind stress over ocean in X and Y direction
@@ -155,5 +215,5 @@
          END DO
 
-       ELSE      ! If no ice dynamics  
+       ELSE      ! no ice dynamics : transmit directly the atmospheric stress to the ocean
                     
           DO jj = 2, jpjm1
@@ -182,5 +242,5 @@
       CALL lbc_lnk( tio_v, 'I', -1. )   ! I-point (i.e. ice U-V point)
 
-      IF( l_ctl .AND. lwp ) THEN 
+      IF(l_ctl) THEN 
          WRITE(numout,*) ' lim_dyn  : tio_u ', SUM( tio_u ), ' tio_v ', SUM( tio_v )
          WRITE(numout,*) ' lim_dyn  : ust2s ', SUM( ust2s )
@@ -189,5 +249,6 @@
    END SUBROUTINE lim_dyn
 
-    SUBROUTINE lim_dyn_init
+
+   SUBROUTINE lim_dyn_init
       !!-------------------------------------------------------------------
       !!                  ***  ROUTINE lim_dyn_init  ***
@@ -219,19 +280,19 @@
          WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics '
          WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*) '   tolerance parameter                              epsd   = ', epsd
-         WRITE(numout,*) '   coefficient for semi-implicit coriolis           alpha  = ', alpha
-         WRITE(numout,*) '   diffusion constant for dynamics                  dm     = ', dm
-         WRITE(numout,*) '   number of sub-time steps for relaxation          nbiter = ', nbiter
-         WRITE(numout,*) '   maximum number of iterations for relaxation      nbitdr = ', nbitdr
-         WRITE(numout,*) '   relaxation constant                              om     = ', om
-         WRITE(numout,*) '   maximum value for the residual of relaxation     resl   = ', resl
-         WRITE(numout,*) '   drag coefficient for oceanic stress              cw     = ', cw
-         WRITE(numout,*) '   turning angle for oceanic stress                 angvg  = ', angvg
-         WRITE(numout,*) '   first bulk-rheology parameter                    pstar  = ', pstar
-         WRITE(numout,*) '   second bulk-rhelogy parameter                    c_rhg  = ', c_rhg
-         WRITE(numout,*) '   minimun value for viscosity                      etamn  = ', etamn
-         WRITE(numout,*) '   creep limit                                      creepl = ', creepl
-         WRITE(numout,*) '   eccentricity of the elliptical yield curve       ecc    = ', ecc
-         WRITE(numout,*) '   horizontal diffusivity coeff. for sea-ice        ahi0   = ', ahi0
+         WRITE(numout,*) '       tolerance parameter                              epsd   = ', epsd
+         WRITE(numout,*) '       coefficient for semi-implicit coriolis           alpha  = ', alpha
+         WRITE(numout,*) '       diffusion constant for dynamics                  dm     = ', dm
+         WRITE(numout,*) '       number of sub-time steps for relaxation          nbiter = ', nbiter
+         WRITE(numout,*) '       maximum number of iterations for relaxation      nbitdr = ', nbitdr
+         WRITE(numout,*) '       relaxation constant                              om     = ', om
+         WRITE(numout,*) '       maximum value for the residual of relaxation     resl   = ', resl
+         WRITE(numout,*) '       drag coefficient for oceanic stress              cw     = ', cw
+         WRITE(numout,*) '       turning angle for oceanic stress                 angvg  = ', angvg
+         WRITE(numout,*) '       first bulk-rheology parameter                    pstar  = ', pstar
+         WRITE(numout,*) '       second bulk-rhelogy parameter                    c_rhg  = ', c_rhg
+         WRITE(numout,*) '       minimun value for viscosity                      etamn  = ', etamn
+         WRITE(numout,*) '       creep limit                                      creepl = ', creepl
+         WRITE(numout,*) '       eccentricity of the elliptical yield curve       ecc    = ', ecc
+         WRITE(numout,*) '       horizontal diffusivity coeff. for sea-ice        ahi0   = ', ahi0
       ENDIF
 
Index: /trunk/NEMO/LIM_SRC/limtrp.F90
===================================================================
--- /trunk/NEMO/LIM_SRC/limtrp.F90	(revision 75)
+++ /trunk/NEMO/LIM_SRC/limtrp.F90	(revision 76)
@@ -33,5 +33,5 @@
    !! * Shared module variables
    REAL(wp), PUBLIC  ::   &  !:
-      bound  = 0.0           !: boundary condit. (0.0 no-slip, 1.0 free-slip)
+      bound  = 0.e0          !: boundary condit. (0.0 no-slip, 1.0 free-slip)
 
    !! * Module variables
@@ -97,5 +97,5 @@
       zsm(:,:) = area(:,:)
       
-      IF( ldyn ) THEN
+      IF( ln_limdyn ) THEN
          !-------------------------------------!
          !   Advection of sea ice properties   !
@@ -228,5 +228,6 @@
          ! Up-dating and limitation of sea ice properties after transport.
          DO jj = 1, jpj
-            zindhe = REAL( MAX( 0, isign(1, jj - jeq ) ) )              !ibug mpp  !!bugmpp  jeq!
+!!!iii      zindhe = REAL( MAX( 0, isign(1, jj - njeq ) ) )              !ibug mpp  !!bugmpp  njeq!
+            zindhe = MAX( 0.e0, SIGN( 1.e0, fcor(1,jj) ) )              ! = 0 for SH, =1 for NH
             DO ji = 1, jpi
 
@@ -298,5 +299,5 @@
          WRITE(numout,*) 'lim_trp_init : Ice parameters for advection '
          WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*) '   boundary conditions (0.0 no-slip, 1.0 free-slip) bound  = ', bound
+         WRITE(numout,*) '   boundary conditions (0. no-slip, 1. free-slip) bound  = ', bound
       ENDIF
             
