Index: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90
===================================================================
--- branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90	(revision 5034)
+++ branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90	(revision 5600)
@@ -6,5 +6,5 @@
    !! history :  1.0  ! 2002-08  (C. Ethe, G. Madec)  original VP code 
    !!            3.0  ! 2007-03  (MA Morales Maqueda, S. Bouillon, M. Vancoppenolle)  LIM3: EVP-Cgrid
-   !!            4.0  ! 2011-02  (G. Madec) dynamical allocation
+   !!            3.5  ! 2011-02  (G. Madec) dynamical allocation
    !!----------------------------------------------------------------------
 #if defined key_lim3
@@ -20,5 +20,4 @@
    USE sbc_ice          ! Surface boundary condition: ice   fields
    USE ice              ! LIM-3 variables
-   USE par_ice          ! LIM-3 parameters
    USE dom_ice          ! LIM-3 domain
    USE limrhg           ! LIM-3 rheology
@@ -31,4 +30,5 @@
    USE timing          ! Timing
    USE limcons        ! conservation tests
+   USE limvar
 
    IMPLICIT NONE
@@ -76,4 +76,6 @@
       CALL wrk_alloc( jpj, zswitch, zmsk )
 
+      CALL lim_var_agg(1)             ! aggregate ice categories
+
       IF( kt == nit000 )   CALL lim_dyn_init   ! Initialization (first time-step only)
 
@@ -83,6 +85,6 @@
          IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
 
-         u_ice_b(:,:) = u_ice(:,:) * tmu(:,:)
-         v_ice_b(:,:) = v_ice(:,:) * tmv(:,:)
+         u_ice_b(:,:) = u_ice(:,:) * umask(:,:,1)
+         v_ice_b(:,:) = v_ice(:,:) * vmask(:,:,1)
 
          ! Rheology (ice dynamics)
@@ -101,5 +103,5 @@
             DO jj = 1, jpj
                zswitch(jj) = SUM( 1.0 - at_i(:,jj) )   ! = REAL(jpj) if ocean everywhere on a j-line
-               zmsk(jj) = SUM( tmask(:,jj,1)    )   ! = 0         if land  everywhere on a j-line
+               zmsk   (jj) = SUM( tmask(:,jj,1)    )   ! = 0         if land  everywhere on a j-line
             END DO
 
@@ -157,9 +159,9 @@
          zv_io(:,:) = v_ice(:,:) - ssv_m(:,:)
          ! frictional velocity at T-point
-         zcoef = 0.5_wp * cw
+         zcoef = 0.5_wp * rn_cio
          DO jj = 2, jpjm1 
             DO ji = fs_2, fs_jpim1   ! vector opt.
                ust2s(ji,jj) = zcoef * (  zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj)   &
-                  &                    + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1)   ) * tms(ji,jj)
+                  &                    + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1) ) * tmask(ji,jj,1)
             END DO
          END DO
@@ -170,9 +172,9 @@
       ELSE      ! no ice dynamics : transmit directly the atmospheric stress to the ocean
          !
-         zcoef = SQRT( 0.5_wp ) / rau0
+         zcoef = SQRT( 0.5_wp ) * r1_rau0
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
                ust2s(ji,jj) = zcoef * SQRT(  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   &
-                  &                        + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1)   ) * tms(ji,jj)
+                  &                        + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) * tmask(ji,jj,1)
             END DO
          END DO
@@ -189,5 +191,5 @@
          CALL prt_ctl(tab2d_1=delta_i   , clinfo1=' lim_dyn  : delta_i   :')
          CALL prt_ctl(tab2d_1=strength  , clinfo1=' lim_dyn  : strength  :')
-         CALL prt_ctl(tab2d_1=area      , clinfo1=' lim_dyn  : cell area :')
+         CALL prt_ctl(tab2d_1=e12t      , clinfo1=' lim_dyn  : cell area :')
          CALL prt_ctl(tab2d_1=at_i      , clinfo1=' lim_dyn  : at_i      :')
          CALL prt_ctl(tab2d_1=vt_i      , clinfo1=' lim_dyn  : vt_i      :')
@@ -241,7 +243,8 @@
       !!-------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namicedyn/ epsd, om, cw, pstar,   &
-         &                c_rhg, creepl, ecc, ahi0,     &
-         &                nevp, relast, alphaevp, hminrhg
+      NAMELIST/namicedyn/ nn_icestr, ln_icestr_bvf, rn_pe_rdg, rn_pstar, rn_crhg, rn_cio,  rn_creepl, rn_ecc, &
+         &                nn_nevp, rn_relast, nn_ahi0, rn_ahi0_ref
+      INTEGER  ::   ji, jj
+      REAL(wp) ::   za00, zd_max
       !!-------------------------------------------------------------------
 
@@ -259,28 +262,64 @@
          WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics '
          WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*) '   tolerance parameter                              epsd   = ', epsd
-         WRITE(numout,*) '   relaxation constant                              om     = ', om
-         WRITE(numout,*) '   drag coefficient for oceanic stress              cw     = ', cw
-         WRITE(numout,*) '   first bulk-rheology parameter                    pstar  = ', pstar
-         WRITE(numout,*) '   second bulk-rhelogy parameter                    c_rhg  = ', c_rhg
-         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,*) '   number of iterations for subcycling              nevp   = ', nevp
-         WRITE(numout,*) '   ratio of elastic timescale over ice time step    relast = ', relast
-         WRITE(numout,*) '   coefficient for the solution of int. stresses  alphaevp = ', alphaevp
-         WRITE(numout,*) '   min ice thickness for rheology calculations     hminrhg = ', hminrhg
+         WRITE(numout,*)'    ice strength parameterization (0=Hibler 1=Rothrock)  nn_icestr     = ', nn_icestr 
+         WRITE(numout,*)'    Including brine volume in ice strength comp.         ln_icestr_bvf = ', ln_icestr_bvf
+         WRITE(numout,*)'   Ratio of ridging work to PotEner change in ridging    rn_pe_rdg     = ', rn_pe_rdg 
+         WRITE(numout,*) '   drag coefficient for oceanic stress                  rn_cio        = ', rn_cio
+         WRITE(numout,*) '   first bulk-rheology parameter                        rn_pstar      = ', rn_pstar
+         WRITE(numout,*) '   second bulk-rhelogy parameter                        rn_crhg       = ', rn_crhg
+         WRITE(numout,*) '   creep limit                                          rn_creepl     = ', rn_creepl
+         WRITE(numout,*) '   eccentricity of the elliptical yield curve           rn_ecc        = ', rn_ecc
+         WRITE(numout,*) '   number of iterations for subcycling                  nn_nevp       = ', nn_nevp
+         WRITE(numout,*) '   ratio of elastic timescale over ice time step        rn_relast     = ', rn_relast
+         WRITE(numout,*) '   horizontal diffusivity calculation                   nn_ahi0       = ', nn_ahi0
+         WRITE(numout,*) '   horizontal diffusivity coeff. (orca2 grid)           rn_ahi0_ref   = ', rn_ahi0_ref
       ENDIF
       !
-      usecc2 = 1._wp / ( ecc * ecc )
-      rhoco  = rau0  * cw
-
-      ! elastic damping
-      telast = relast * rdt_ice
-
-      !  Diffusion coefficients.
-      ahiu(:,:) = ahi0 * umask(:,:,1)
-      ahiv(:,:) = ahi0 * vmask(:,:,1)
-      !
+      usecc2 = 1._wp / ( rn_ecc * rn_ecc )
+      rhoco  = rau0  * rn_cio
+      !
+      !  Diffusion coefficients
+      SELECT CASE( nn_ahi0 )
+
+      CASE( 0 )
+         ahiu(:,:) = rn_ahi0_ref
+         ahiv(:,:) = rn_ahi0_ref
+
+         IF(lwp) WRITE(numout,*) ''
+         IF(lwp) WRITE(numout,*) '   laplacian operator: ahim constant = rn_ahi0_ref'
+
+      CASE( 1 ) 
+
+         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )
+         IF( lk_mpp )   CALL mpp_max( zd_max )          ! max over the global domain
+         
+         ahiu(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp   ! 1.e05 = 100km = max grid space at 60° latitude in orca2
+                                                        !                    (60° = min latitude for ice cover)  
+         ahiv(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp
+
+         IF(lwp) WRITE(numout,*) ''
+         IF(lwp) WRITE(numout,*) '   laplacian operator: ahim proportional to max of e1 e2 over the domain (', zd_max, ')'
+         IF(lwp) WRITE(numout,*) '   value for ahim = ', rn_ahi0_ref * zd_max * 1.e-05_wp 
+         
+      CASE( 2 ) 
+
+         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )
+         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain
+         
+         za00 = rn_ahi0_ref * 1.e-05_wp          ! 1.e05 = 100km = max grid space at 60° latitude in orca2
+                                                 !                    (60° = min latitude for ice cover)  
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ahiu(ji,jj) = za00 * MAX( e1t(ji,jj), e2t(ji,jj) ) * umask(ji,jj,1)
+               ahiv(ji,jj) = za00 * MAX( e1f(ji,jj), e2f(ji,jj) ) * vmask(ji,jj,1)
+            END DO
+         END DO
+         !
+         IF(lwp) WRITE(numout,*) ''
+         IF(lwp) WRITE(numout,*) '   laplacian operator: ahim proportional to e1'
+         IF(lwp) WRITE(numout,*) '   maximum grid-spacing = ', zd_max, ' maximum value for ahim = ', za00*zd_max
+         
+      END SELECT
+
    END SUBROUTINE lim_dyn_init
 
