Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn.F90	(revision 8516)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn.F90	(revision 8517)
@@ -34,11 +34,15 @@
    PUBLIC   ice_dyn_init   ! called by icestp.F90
    
-   INTEGER ::              nice_dyn   ! choice of the type of advection scheme
+   INTEGER ::              nice_dyn   ! choice of the type of dynamics
    !                                        ! associated indices:
-   INTEGER, PARAMETER ::   np_dynNO   = 0   ! no ice dynamics and ice advection
-   INTEGER, PARAMETER ::   np_dynFULL = 1   ! full ice dynamics  (rheology + advection + ridging/rafting + correction)
-   INTEGER, PARAMETER ::   np_dyn     = 2   ! no ridging/rafting (rheology + advection                   + correction)
-   INTEGER, PARAMETER ::   np_dynPURE = 3   ! pure dynamics      (rheology + advection) 
-
+   INTEGER, PARAMETER ::   np_dynFULL    = 1   ! full ice dynamics               (rheology + advection + ridging/rafting + correction)
+   INTEGER, PARAMETER ::   np_dynRHGADV1 = 2   ! no ridging/rafting              (rheology + advection                   + correction)
+   INTEGER, PARAMETER ::   np_dynRHGADV2 = 3   ! pure dynamics                   (rheology + advection) 
+   INTEGER, PARAMETER ::   np_dynADV     = 4   ! only advection w prescribed vel.(rn_uvice + advection)
+   !
+   ! ** namelist (namdyn) **
+   REAL(wp) ::   rn_uice          ! prescribed u-vel (case np_dynADV)
+   REAL(wp) ::   rn_vice          ! prescribed v-vel (case np_dynADV)
+   
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
@@ -61,5 +65,6 @@
       INTEGER, INTENT(in) ::   kt     ! ice time step
       !!
-      INTEGER  ::   jl   ! dummy loop indices
+      INTEGER ::   ji, jj, jl         ! dummy loop indices
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zhmax
       !!--------------------------------------------------------------------
       !
@@ -72,7 +77,6 @@
       ENDIF
 
-      CALL ice_var_agg(1)           ! -- aggregate ice categories
       !                      
-      IF( ln_landfast ) THEN        ! -- Landfast ice parameterization: define max bottom friction
+      IF( ln_landfast ) THEN            !-- Landfast ice parameterization: define max bottom friction
          tau_icebfr(:,:) = 0._wp
          DO jl = 1, jpl
@@ -81,26 +85,41 @@
       ENDIF
 
-      SELECT CASE( nice_dyn )       ! -- Set which dynamics is running
+      zhmax(:,:,:) = ht_i_b(:,:,:)      !-- Record max of the surrounding 9-pts ice thick. (for CALL Hbig)
+      DO jl = 1, jpl
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+!!gm use of MAXVAL here is very probably less efficient than expending the 9 values
+               zhmax(ji,jj,jl) = MAX( epsi20, MAXVAL( ht_i_b(ji-1:ji+1,jj-1:jj+1,jl) ) )
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( zhmax(:,:,:), 'T', 1. )
+      !
+      !
+      SELECT CASE( nice_dyn )           !-- Set which dynamics is running
 
       CASE ( np_dynFULL )          !==  all dynamical processes  ==!
-         CALL ice_rhg   ( kt )          ! -- rheology  
-         CALL ice_adv   ( kt )          ! -- advection of ice
-         CALL ice_rdgrft( kt )          ! -- ridging/rafting 
-         CALL ice_cor   ( kt , 1 )      ! -- Corrections
-
-      CASE ( np_dyn )              !==  pure dynamics only ==!   (no ridging/rafting)   (nono cat. case 2)
-         CALL ice_rhg   ( kt )          ! -- rheology  
-         CALL ice_adv   ( kt )          ! -- advection of ice
-         CALL ice_cor   ( kt , 1 )      ! -- Corrections
-
-      CASE ( np_dynPURE )          !==  pure dynamics only ==!   (nn_icedyn= 1 )
-         CALL ice_rhg   ( kt )          ! -- rheology  
-         CALL ice_adv   ( kt )          ! -- advection of ice
-
-      CASE ( np_dynNO )            !==  prescribed ice velocities ==!   (nn_icedyn= 0 )
+         CALL ice_rhg   ( kt )                            ! -- rheology  
+         CALL ice_adv   ( kt )   ;   CALL Hbig( zhmax )   ! -- advection of ice + correction on ice thickness
+         CALL ice_rdgrft( kt )                            ! -- ridging/rafting 
+         CALL ice_cor   ( kt , 1 )                        ! -- Corrections
+
+      CASE ( np_dynRHGADV1 )       !==  no ridge/raft ==!   (mono cat. case 2)
+         CALL ice_rhg   ( kt )                            ! -- rheology  
+         CALL ice_adv   ( kt )                            ! -- advection of ice
+         CALL Hpiling                                     ! -- simple pile-up (replaces ridging/rafting)
+         CALL ice_cor   ( kt , 1 )                        ! -- Corrections
+
+      CASE ( np_dynRHGADV2 )       !==  no ridge/raft & no corrections ==!
+         CALL ice_rhg   ( kt )                            ! -- rheology  
+         CALL ice_adv   ( kt )                            ! -- advection of ice
+         CALL Hpiling                                     ! -- simple pile-up (replaces ridging/rafting)
+
+      CASE ( np_dynADV )           !==  pure advection ==!   (prescribed velocities)
          u_ice(:,:) = rn_uice * umask(:,:,1)
          v_ice(:,:) = rn_vice * vmask(:,:,1)
          !!CALL RANDOM_NUMBER(u_ice(:,:))
          !!CALL RANDOM_NUMBER(v_ice(:,:))
+         CALL ice_adv   ( kt )                            ! -- advection of ice
 
       END SELECT
@@ -110,4 +129,73 @@
    END SUBROUTINE ice_dyn
 
+   SUBROUTINE Hbig( phmax )
+      !!-------------------------------------------------------------------
+      !!                  ***  ROUTINE Hbig  ***
+      !!
+      !! ** Purpose : Thickness correction in case advection scheme creates
+      !!              abnormally tick ice
+      !!
+      !! ** Method  : 1- check whether ice thickness resulting from advection is
+      !!                 larger than the surrounding 9-points before advection
+      !!                 and reduce it if a) divergence or b) convergence & at_i>0.8
+      !!              2- bound ice thickness with hi_max (99m)
+      !!
+      !! ** input   : Max thickness of the surrounding 9-points
+      !!-------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   phmax   ! max ice thick from surrounding 9-pts
+      !
+      INTEGER  ::   ji, jj, jl         ! dummy loop indices
+      REAL(wp) ::   zh, zdv
+      !!-------------------------------------------------------------------
+      !
+      CALL ice_var_zapsmall                       !-- zap small areas
+      !
+      DO jl = 1, jpl
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF ( v_i(ji,jj,jl) > 0._wp ) THEN  !-- bound to hmax
+                  !
+                  zh  = v_i (ji,jj,jl) / a_i(ji,jj,jl)
+                  zdv = v_i(ji,jj,jl) - v_i_b(ji,jj,jl)  
+                  !
+                  IF ( ( zdv >  0.0 .AND. zh > phmax(ji,jj,jl) .AND. at_i_b(ji,jj) < 0.80 ) .OR. &
+                     & ( zdv <= 0.0 .AND. zh > phmax(ji,jj,jl) ) ) THEN
+                     a_i (ji,jj,jl) = v_i(ji,jj,jl) / MIN( phmax(ji,jj,jl), hi_max(jpl) )   !-- bound ht_i to hi_max (99 m)
+                  ENDIF
+                  !
+               ENDIF
+            END DO
+         END DO
+      END DO
+            
+      IF ( nn_pnd_scheme > 0 ) THEN               !-- correct pond fraction to avoid a_ip > a_i
+         WHERE( a_ip(:,:,:) > a_i(:,:,:) )   a_ip(:,:,:) = a_i(:,:,:)
+      ENDIF
+      !
+   END SUBROUTINE Hbig
+
+   SUBROUTINE Hpiling
+      !!-------------------------------------------------------------------
+      !!                  ***  ROUTINE Hpiling  ***
+      !!
+      !! ** Purpose : Simple conservative piling comparable with 1-cat models
+      !!
+      !! ** Method  : pile-up ice when no ridging/rafting
+      !!
+      !! ** input   : a_i
+      !!-------------------------------------------------------------------
+      INTEGER ::   jl         ! dummy loop indices
+      !!-------------------------------------------------------------------
+      !
+      CALL ice_var_zapsmall                       !-- zap small areas
+      !
+      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
+      DO jl = 1, jpl
+         WHERE( at_i(:,:) > epsi20 )
+            a_i(:,:,jl) = a_i(:,:,jl) * (  1._wp + MIN( rn_amax_2d(:,:) - at_i(:,:) , 0._wp ) / at_i(:,:)  )
+         END WHERE
+      END DO
+      !
+   END SUBROUTINE Hpiling
 
    SUBROUTINE ice_dyn_init
@@ -123,7 +211,7 @@
       !! ** input   :   Namelist namice_dyn
       !!-------------------------------------------------------------------
-      INTEGER ::   ios   ! Local integer output status for namelist read
-      !!
-      NAMELIST/namice_dyn/ ln_icedyn  , nn_icedyn, rn_uice  , rn_vice ,    &
+      INTEGER ::   ios, ioptio   ! Local integer output status for namelist read
+      !!
+      NAMELIST/namice_dyn/ ln_dynFULL, ln_dynRHGADV, ln_dynADV, rn_uice, rn_vice,  &
          &                 rn_ishlat  , rn_cio   ,                         &
          &                 ln_landfast, rn_gamma , rn_icebfr, rn_lfrelax
@@ -144,9 +232,8 @@
          WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namice_dyn'
-         WRITE(numout,*) '      Ice dynamics       (T) or not (F)                         ln_icedyn  = ', ln_icedyn
-         WRITE(numout,*) '         associated switch                                      nn_icedyn  = ', nn_icedyn
-         WRITE(numout,*) '            =2 all processes (default option)'
-         WRITE(numout,*) '            =1 advection only (no ridging/rafting)'
-         WRITE(numout,*) '            =0 advection only with prescribed velocity given by '
+         WRITE(numout,*) '      Full ice dynamics      (rhg + adv + ridge/raft + corr)  ln_dynFULL   = ', ln_dynFULL
+         WRITE(numout,*) '      No ridge/raft & No cor (rhg + adv)                      ln_dynRHGADV = ', ln_dynRHGADV
+         WRITE(numout,*) '      Advection only         (rn_uvice + adv)                 ln_dynADV    = ', ln_dynADV
+         WRITE(numout,*) '           with prescribed velocity given by '
          WRITE(numout,*) '               a uniform field               (u,v)_ice = (rn_uice,rn_vice) = (', rn_uice,',', rn_vice,')'
          WRITE(numout,*) '      lateral boundary condition for sea ice dynamics        rn_ishlat     = ', rn_ishlat
@@ -158,16 +245,15 @@
       ENDIF
       !                             !== set the choice of ice dynamics ==!
-      SELECT CASE( nn_icedyn )
-      CASE( 2 )                    
-         IF( nn_monocat /= 2 ) THEN          !--- full dynamics (rheology + advection + ridging/rafting + correction)
-            nice_dyn = np_dynFULL
-         ELSE
-            nice_dyn = np_dyn                !--- dynamics without ridging/rafting
-         ENDIF
-      CASE( 1 )                              !--- dynamics without ridging/rafting and correction 
-         nice_dyn = np_dynPURE
-      CASE( 0 )                              !--- prescribed ice velocities (from namelist) 
-         nice_dyn = np_dynNO
-      END SELECT
+      ioptio = 0 
+      !      !--- full dynamics                               (rheology + advection + ridging/rafting + correction)
+      IF( ln_dynFULL .AND. nn_monocat /= 2 ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynFULL      ;   ENDIF
+      !      !--- dynamics without ridging/rafting            (rheology + advection                   + correction)
+      IF( ln_dynFULL .AND. nn_monocat == 2 ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynRHGADV1   ;   ENDIF
+      !      !--- dynamics without ridging/rafting and corr   (rheology + advection)
+      IF( ln_dynRHGADV                     ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynRHGADV2   ;   ENDIF
+      !      !--- advection only with prescribed ice velocities (from namelist)
+      IF( ln_dynADV                        ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynADV       ;   ENDIF
+      !
+      IF( ioptio /= 1 )    CALL ctl_stop( 'ice_dyn_init: one and only one ice dynamics option has to be defined ' )
       !
       !                                      !--- Lateral boundary conditions
@@ -177,11 +263,7 @@
       ELSEIF ( 2. < rn_ishlat                      ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  strong-slip'
       ENDIF
-      !
       !                                      !--- NO Landfast ice : set to zero once for all
       IF( .NOT. ln_landfast )   tau_icebfr(:,:) = 0._wp 
       !
-      !                                      !--- simple conservative piling, comparable with LIM2
-      l_piling = nn_icedyn == 1 .OR. ( nn_monocat == 2  .AND.  jpl == 1 )
-      !
    END SUBROUTINE ice_dyn_init
 
