Index: NEMO/branches/UKMO/NEMO_4.0_mirror/src/ICE/icedyn.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_mirror/src/ICE/icedyn.F90	(revision 10888)
+++ NEMO/branches/UKMO/NEMO_4.0_mirror/src/ICE/icedyn.F90	(revision 11081)
@@ -74,8 +74,7 @@
       INTEGER, INTENT(in) ::   kt     ! ice time step
       !!
-      INTEGER  ::   ji, jj, jl        ! dummy loop indices
+      INTEGER  ::   ji, jj        ! dummy loop indices
       REAL(wp) ::   zcoefu, zcoefv
-      REAL(wp),              DIMENSION(jpi,jpj,jpl) ::   zhi_max, zhs_max, zhip_max
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)         ::   zdivu_i
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdivu_i
       !!--------------------------------------------------------------------
       !
@@ -89,48 +88,38 @@
       ENDIF
       !                      
-      IF( ln_landfast_home ) THEN      !-- Landfast ice parameterization
-         tau_icebfr(:,:) = 0._wp
-         DO jl = 1, jpl
-            WHERE( h_i_b(:,:,jl) > ht_n(:,:) * rn_depfra )   tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr
-         END DO
-      ENDIF
-      !
-      !                                !-- Record max of the surrounding 9-pts ice thick. (for CALL Hbig)
-      DO jl = 1, jpl
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1
-               zhip_max(ji,jj,jl) = MAX( epsi20, h_ip_b(ji,jj,jl), h_ip_b(ji+1,jj  ,jl), h_ip_b(ji  ,jj+1,jl), &
-                  &                                                h_ip_b(ji-1,jj  ,jl), h_ip_b(ji  ,jj-1,jl), &
-                  &                                                h_ip_b(ji+1,jj+1,jl), h_ip_b(ji-1,jj-1,jl), &
-                  &                                                h_ip_b(ji+1,jj-1,jl), h_ip_b(ji-1,jj+1,jl) )
-               zhi_max (ji,jj,jl) = MAX( epsi20, h_i_b (ji,jj,jl), h_i_b (ji+1,jj  ,jl), h_i_b (ji  ,jj+1,jl), &
-                  &                                                h_i_b (ji-1,jj  ,jl), h_i_b (ji  ,jj-1,jl), &
-                  &                                                h_i_b (ji+1,jj+1,jl), h_i_b (ji-1,jj-1,jl), &
-                  &                                                h_i_b (ji+1,jj-1,jl), h_i_b (ji-1,jj+1,jl) )
-               zhs_max (ji,jj,jl) = MAX( epsi20, h_s_b (ji,jj,jl), h_s_b (ji+1,jj  ,jl), h_s_b (ji  ,jj+1,jl), &
-                  &                                                h_s_b (ji-1,jj  ,jl), h_s_b (ji  ,jj-1,jl), &
-                  &                                                h_s_b (ji+1,jj+1,jl), h_s_b (ji-1,jj-1,jl), &
-                  &                                                h_s_b (ji+1,jj-1,jl), h_s_b (ji-1,jj+1,jl) )
-            END DO
-         END DO
-      END DO
-      CALL lbc_lnk_multi( 'icedyn', zhi_max(:,:,:), 'T', 1., zhs_max(:,:,:), 'T', 1., zhip_max(:,:,:), 'T', 1. )
-      !
-      !
-      SELECT CASE( nice_dyn )           !-- Set which dynamics is running
+      ! retrieve thickness from volume for landfast param. and UMx advection scheme
+      WHERE( a_i(:,:,:) >= epsi20 )
+         h_i(:,:,:) = v_i(:,:,:) / a_i_b(:,:,:)
+         h_s(:,:,:) = v_s(:,:,:) / a_i_b(:,:,:)
+      ELSEWHERE
+         h_i(:,:,:) = 0._wp
+         h_s(:,:,:) = 0._wp
+      END WHERE
+      !
+      WHERE( a_ip(:,:,:) >= epsi20 )
+         h_ip(:,:,:) = v_ip(:,:,:) / a_ip(:,:,:)
+      ELSEWHERE
+         h_ip(:,:,:) = 0._wp
+      END WHERE
+      !
+      !
+      SELECT CASE( nice_dyn )          !-- Set which dynamics is running
 
       CASE ( np_dynALL )           !==  all dynamical processes  ==!
-         CALL ice_dyn_rhg   ( kt )                                                 ! -- rheology  
-         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max, zhip_max )   ! -- advection of ice + correction on ice thickness
-         CALL ice_dyn_rdgrft( kt )                                                 ! -- ridging/rafting 
-         CALL ice_cor       ( kt , 1 )                                             ! -- Corrections
-
+         !
+         CALL ice_dyn_rhg   ( kt )                                          ! -- rheology  
+         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice
+         CALL ice_dyn_rdgrft( kt )                                          ! -- ridging/rafting 
+         CALL ice_cor       ( kt , 1 )                                      ! -- Corrections
+         !
       CASE ( np_dynRHGADV  )       !==  no ridge/raft & no corrections ==!
-         CALL ice_dyn_rhg   ( kt )                                                 ! -- rheology  
-         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max, zhip_max )   ! -- advection of ice + correction on ice thickness
-         CALL Hpiling                                                              ! -- simple pile-up (replaces ridging/rafting)
-
+         !
+         CALL ice_dyn_rhg   ( kt )                                          ! -- rheology  
+         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice
+         CALL Hpiling                                                       ! -- simple pile-up (replaces ridging/rafting)
+         CALL ice_var_zapsmall                                              ! -- zap small areas
+         !
       CASE ( np_dynADV1D )         !==  pure advection ==!   (1D)
-         ALLOCATE( zdivu_i(jpi,jpj) )
+         !
          ! --- monotonicity test from Schar & Smolarkiewicz 1996 --- !
          ! CFL = 0.5 at a distance from the bound of 1/6 of the basin length
@@ -145,20 +134,8 @@
          END DO
          ! ---
-         CALL ice_dyn_adv   ( kt )                                       ! -- advection of ice
-
-         ! diagnostics: divergence at T points 
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-               zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
-                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj)
-            END DO
-         END DO
-         CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. )
-         IF( iom_use('icediv') )   CALL iom_put( "icediv" , zdivu_i(:,:) )
-
-         DEALLOCATE( zdivu_i )
-
+         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice
+         !
       CASE ( np_dynADV2D )         !==  pure advection ==!   (2D w prescribed velocities)
-         ALLOCATE( zdivu_i(jpi,jpj) )
+         !
          u_ice(:,:) = rn_uice * umask(:,:,1)
          v_ice(:,:) = rn_vice * vmask(:,:,1)
@@ -166,19 +143,30 @@
          !CALL RANDOM_NUMBER(v_ice(:,:)) ; v_ice(:,:) = v_ice(:,:) * 0.1 + rn_vice * 0.9 * vmask(:,:,1)
          ! ---
-         CALL ice_dyn_adv   ( kt )                                       ! -- advection of ice
-
-         ! diagnostics: divergence at T points 
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-               zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
-                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj)
+         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice
+
+      END SELECT
+      !
+      !
+      ! diagnostics: divergence at T points 
+      IF( iom_use('icediv') ) THEN
+         !
+         SELECT CASE( nice_dyn )
+
+         CASE ( np_dynADV1D , np_dynADV2D )
+
+            ALLOCATE( zdivu_i(jpi,jpj) )
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1
+                  zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
+                     &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj)
+               END DO
             END DO
-         END DO
-         CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. )
-         IF( iom_use('icediv') )   CALL iom_put( "icediv" , zdivu_i(:,:) )
-
-         DEALLOCATE( zdivu_i )
-
-      END SELECT
+            CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. )
+            CALL iom_put( "icediv" , zdivu_i(:,:) )
+            DEALLOCATE( zdivu_i )
+
+         END SELECT
+         !
+      ENDIF
       !
       ! controls
@@ -188,92 +176,4 @@
 
 
-   SUBROUTINE Hbig( phi_max, phs_max, phip_max )
-      !!-------------------------------------------------------------------
-      !!                  ***  ROUTINE Hbig  ***
-      !!
-      !! ** Purpose : Thickness correction in case advection scheme creates
-      !!              abnormally tick ice or snow
-      !!
-      !! ** Method  : 1- check whether ice thickness is larger than the surrounding 9-points
-      !!                 (before advection) and reduce it by adapting ice concentration
-      !!              2- check whether snow thickness is larger than the surrounding 9-points
-      !!                 (before advection) and reduce it by sending the excess in the ocean
-      !!              3- check whether snow load deplets the snow-ice interface below sea level$
-      !!                 and reduce it by sending the excess in the ocean
-      !!              4- correct pond fraction to avoid a_ip > a_i
-      !!
-      !! ** input   : Max thickness of the surrounding 9-points
-      !!-------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   phi_max, phs_max, phip_max   ! max ice thick from surrounding 9-pts
-      !
-      INTEGER  ::   ji, jj, jl         ! dummy loop indices
-      REAL(wp) ::   zhip, zhi, zhs, zvs_excess, zfra
-      !!-------------------------------------------------------------------
-      ! controls
-      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'Hbig', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
-      !
-      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
-                  !
-                  !                               ! -- check h_ip -- !
-                  ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip
-                  IF( ln_pnd_H12 .AND. v_ip(ji,jj,jl) > 0._wp ) THEN
-                     zhip = v_ip(ji,jj,jl) / MAX( epsi20, a_ip(ji,jj,jl) )
-                     IF( zhip > phip_max(ji,jj,jl) .AND. a_ip(ji,jj,jl) < 0.15 ) THEN
-                        a_ip(ji,jj,jl) = v_ip(ji,jj,jl) / phip_max(ji,jj,jl)
-                     ENDIF
-                  ENDIF
-                  !
-                  !                               ! -- check h_i -- !
-                  ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i
-                  zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
-                  IF( zhi > phi_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN
-                     a_i(ji,jj,jl) = v_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) )   !-- bound h_i to hi_max (99 m)
-                  ENDIF
-                  !
-                  !                               ! -- check h_s -- !
-                  ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean
-                  zhs = v_s(ji,jj,jl) / a_i(ji,jj,jl)
-                  IF( v_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN
-                     zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 )
-                     !
-                     wfx_res(ji,jj) = wfx_res(ji,jj) + ( v_s(ji,jj,jl) - a_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * r1_rdtice
-                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( e_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0
-                     !
-                     e_s(ji,jj,1:nlay_s,jl) = e_s(ji,jj,1:nlay_s,jl) * zfra
-                     v_s(ji,jj,jl)          = a_i(ji,jj,jl) * phs_max(ji,jj,jl)
-                  ENDIF           
-                  !
-                  !                               ! -- check snow load -- !
-                  ! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean
-                  !    this correction is crucial because of the call to routine icecor afterwards which imposes a mini of ice thick. (rn_himin)
-                  !    this imposed mini can artificially make the snow very thick (if concentration decreases drastically)
-                  zvs_excess = MAX( 0._wp, v_s(ji,jj,jl) - v_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos )
-                  IF( zvs_excess > 0._wp ) THEN
-                     zfra = ( v_s(ji,jj,jl) - zvs_excess ) / MAX( v_s(ji,jj,jl), epsi20 )
-                     wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * r1_rdtice
-                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( e_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0
-                     !
-                     e_s(ji,jj,1:nlay_s,jl) = e_s(ji,jj,1:nlay_s,jl) * zfra
-                     v_s(ji,jj,jl)          = v_s(ji,jj,jl) - zvs_excess
-                  ENDIF
-                  
-               ENDIF
-            END DO
-         END DO
-      END DO 
-      !                                           !-- correct pond fraction to avoid a_ip > a_i
-      WHERE( a_ip(:,:,:) > a_i(:,:,:) )   a_ip(:,:,:) = a_i(:,:,:)
-      !
-      ! controls
-      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'Hbig', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
-      !
-   END SUBROUTINE Hbig
-
-
    SUBROUTINE Hpiling
       !!-------------------------------------------------------------------
@@ -290,6 +190,4 @@
       ! controls
       IF( ln_icediachk )   CALL ice_cons_hsm(0, 'Hpiling', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
-      !
-      CALL ice_var_zapsmall                       !-- zap small areas
       !
       at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
@@ -337,16 +235,16 @@
          WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namdyn:'
-         WRITE(numout,*) '      Full ice dynamics      (rhg + adv + ridge/raft + corr)  ln_dynALL       = ', ln_dynALL
-         WRITE(numout,*) '      No ridge/raft & No cor (rhg + adv)                      ln_dynRHGADV    = ', ln_dynRHGADV
-         WRITE(numout,*) '      Advection 1D only      (Schar & Smolarkiewicz 1996)     ln_dynADV1D     = ', ln_dynADV1D
-         WRITE(numout,*) '      Advection 2D only      (rn_uvice + adv)                 ln_dynADV2D     = ', ln_dynADV2D
-         WRITE(numout,*) '         with prescribed velocity given by   (u,v)_ice = (rn_uice,rn_vice)    = (', rn_uice,',', rn_vice,')'
-         WRITE(numout,*) '      lateral boundary condition for sea ice dynamics         rn_ishlat       = ', rn_ishlat
-         WRITE(numout,*) '      Landfast: param from Lemieux 2016                       ln_landfast_L16 = ', ln_landfast_L16
-         WRITE(numout,*) '      Landfast: param from home made                          ln_landfast_home= ', ln_landfast_home
-         WRITE(numout,*) '         fraction of ocean depth that ice must reach          rn_depfra       = ', rn_depfra
-         WRITE(numout,*) '         maximum bottom stress per unit area of contact       rn_icebfr       = ', rn_icebfr
-         WRITE(numout,*) '         relax time scale (s-1) to reach static friction      rn_lfrelax      = ', rn_lfrelax
-         WRITE(numout,*) '         isotropic tensile strength                           rn_tensile      = ', rn_tensile
+         WRITE(numout,*) '      Full ice dynamics      (rhg + adv + ridge/raft + corr) ln_dynALL       = ', ln_dynALL
+         WRITE(numout,*) '      No ridge/raft & No cor (rhg + adv)                     ln_dynRHGADV    = ', ln_dynRHGADV
+         WRITE(numout,*) '      Advection 1D only      (Schar & Smolarkiewicz 1996)    ln_dynADV1D     = ', ln_dynADV1D
+         WRITE(numout,*) '      Advection 2D only      (rn_uvice + adv)                ln_dynADV2D     = ', ln_dynADV2D
+         WRITE(numout,*) '         with prescribed velocity given by   (u,v)_ice = (rn_uice,rn_vice)   = (', rn_uice,',',rn_vice,')'
+         WRITE(numout,*) '      lateral boundary condition for sea ice dynamics        rn_ishlat       = ', rn_ishlat
+         WRITE(numout,*) '      Landfast: param from Lemieux 2016                      ln_landfast_L16 = ', ln_landfast_L16
+         WRITE(numout,*) '      Landfast: param from home made                         ln_landfast_home= ', ln_landfast_home
+         WRITE(numout,*) '         fraction of ocean depth that ice must reach         rn_depfra       = ', rn_depfra
+         WRITE(numout,*) '         maximum bottom stress per unit area of contact      rn_icebfr       = ', rn_icebfr
+         WRITE(numout,*) '         relax time scale (s-1) to reach static friction     rn_lfrelax      = ', rn_lfrelax
+         WRITE(numout,*) '         isotropic tensile strength                          rn_tensile      = ', rn_tensile
          WRITE(numout,*)
       ENDIF
