Index: NEMO/branches/2019/ENHANCE-03_closea/src/ICE/icedyn_rdgrft.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-03_closea/src/ICE/icedyn_rdgrft.F90	(revision 11207)
+++ NEMO/branches/2019/ENHANCE-03_closea/src/ICE/icedyn_rdgrft.F90	(revision 12149)
@@ -86,7 +86,7 @@
       !!                ***  ROUTINE ice_dyn_rdgrft_alloc ***
       !!-------------------------------------------------------------------
-      ALLOCATE( closing_net(jpij), opning(jpij)   , closing_gross(jpij),   &
-         &      apartf(jpij,0:jpl), hrmin(jpij,jpl), hraft(jpij,jpl)    , aridge(jpij,jpl),   &
-         &      hrmax(jpij,jpl), hi_hrdg(jpij,jpl)  , araft (jpij,jpl),  &
+      ALLOCATE( closing_net(jpij)  , opning(jpij)      , closing_gross(jpij) ,               &
+         &      apartf(jpij,0:jpl) , hrmin  (jpij,jpl) , hraft(jpij,jpl) , aridge(jpij,jpl), &
+         &      hrmax (jpij,jpl)   , hi_hrdg(jpij,jpl) , araft(jpij,jpl) ,                   &
          &      ze_i_2d(jpij,nlay_i,jpl), ze_s_2d(jpij,nlay_s,jpl), STAT=ice_dyn_rdgrft_alloc )
 
@@ -137,5 +137,4 @@
       REAL(wp) ::   zfac                       ! local scalar
       INTEGER , DIMENSION(jpij) ::   iptidx        ! compute ridge/raft or not
-      REAL(wp), DIMENSION(jpij) ::   zdivu_adv     ! divu as implied by transport scheme  (1/s)
       REAL(wp), DIMENSION(jpij) ::   zdivu, zdelt  ! 1D divu_i & delta_i
       !
@@ -145,4 +144,5 @@
       IF( ln_timing    )   CALL timing_start('icedyn_rdgrft')                                                             ! timing
       IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
+      IF( ln_icediachk )   CALL ice_cons2D  (0, 'icedyn_rdgrft',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation
 
       IF( kt == nit000 ) THEN
@@ -174,7 +174,7 @@
         
          ! just needed here
-         CALL tab_2d_1d( npti, nptidx(1:npti), zdivu   (1:npti)      , divu_i  )
          CALL tab_2d_1d( npti, nptidx(1:npti), zdelt   (1:npti)      , delta_i )
          ! needed here and in the iteration loop
+         CALL tab_2d_1d( npti, nptidx(1:npti), zdivu   (1:npti)      , divu_i) ! zdivu is used as a work array here (no change in divu_i)
          CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d  (1:npti,1:jpl), a_i   )
          CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d  (1:npti,1:jpl), v_i   )
@@ -186,16 +186,8 @@
             closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp )
             !
-            ! divergence given by the advection scheme
-            !   (which may not be equal to divu as computed from the velocity field)
-            IF    ( ln_adv_Pra ) THEN
-               zdivu_adv(ji) = ( 1._wp - ato_i_1d(ji) - SUM( a_i_2d(ji,:) ) ) * r1_rdtice
-            ELSEIF( ln_adv_UMx ) THEN
-               zdivu_adv(ji) = zdivu(ji)
-            ENDIF
-            !
-            IF( zdivu_adv(ji) < 0._wp )   closing_net(ji) = MAX( closing_net(ji), -zdivu_adv(ji) )   ! make sure the closing rate is large enough
-            !                                                                                        ! to give asum = 1.0 after ridging
+            IF( zdivu(ji) < 0._wp )   closing_net(ji) = MAX( closing_net(ji), -zdivu(ji) )   ! make sure the closing rate is large enough
+            !                                                                                ! to give asum = 1.0 after ridging
             ! Opening rate (non-negative) that will give asum = 1.0 after ridging.
-            opning(ji) = closing_net(ji) + zdivu_adv(ji)
+            opning(ji) = closing_net(ji) + zdivu(ji)
          END DO
          !
@@ -214,5 +206,5 @@
                ato_i_1d   (ipti)   = ato_i_1d   (ji)
                closing_net(ipti)   = closing_net(ji)
-               zdivu_adv  (ipti)   = zdivu_adv  (ji)
+               zdivu      (ipti)   = zdivu      (ji)
                opning     (ipti)   = opning     (ji)
             ENDIF
@@ -258,7 +250,7 @@
                ELSE
                   iterate_ridging  = 1
-                  zdivu_adv  (ji) = zfac * r1_rdtice
-                  closing_net(ji) = MAX( 0._wp, -zdivu_adv(ji) )
-                  opning     (ji) = MAX( 0._wp,  zdivu_adv(ji) )
+                  zdivu      (ji) = zfac * r1_rdtice
+                  closing_net(ji) = MAX( 0._wp, -zdivu(ji) )
+                  opning     (ji) = MAX( 0._wp,  zdivu(ji) )
                ENDIF
             END DO
@@ -276,6 +268,8 @@
 
       ! controls
+      IF( ln_ctl       )   CALL ice_prt3D   ('icedyn_rdgrft')                                                             ! prints
+      IF( ln_icectl    )   CALL ice_prt     (kt, iiceprt, jiceprt,-1, ' - ice dyn rdgrft - ')                             ! prints
       IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
-      IF( ln_ctl       )   CALL ice_prt3D   ('icedyn_rdgrft')                                                             ! prints
+      IF( ln_icediachk )   CALL ice_cons2D  (1, 'icedyn_rdgrft',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation
       IF( ln_timing    )   CALL timing_stop ('icedyn_rdgrft')                                                             ! timing
       !
@@ -306,5 +300,5 @@
 
       !                       ! Ice thickness needed for rafting
-      WHERE( pa_i(1:npti,:) > epsi20 )   ;   zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:)
+      WHERE( pa_i(1:npti,:) > epsi10 )   ;   zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:)
       ELSEWHERE                          ;   zhi(1:npti,:) = 0._wp
       END WHERE
@@ -325,5 +319,5 @@
       zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 )
       !
-      WHERE( zasum(1:npti) > epsi20 )   ;   z1_asum(1:npti) = 1._wp / zasum(1:npti)
+      WHERE( zasum(1:npti) > epsi10 )   ;   z1_asum(1:npti) = 1._wp / zasum(1:npti)
       ELSEWHERE                         ;   z1_asum(1:npti) = 0._wp
       END WHERE
@@ -451,5 +445,5 @@
       ! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate.  
       ! NOTE: 0 < aksum <= 1
-      WHERE( zaksum(1:npti) > epsi20 )   ;   closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti)
+      WHERE( zaksum(1:npti) > epsi10 )   ;   closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti)
       ELSEWHERE                          ;   closing_gross(1:npti) = 0._wp
       END WHERE
@@ -534,5 +528,5 @@
             IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN   ! only if ice is ridging
 
-               IF( a_i_2d(ji,jl1) > epsi20 ) THEN   ;   z1_ai(ji) = 1._wp / a_i_2d(ji,jl1)
+               IF( a_i_2d(ji,jl1) > epsi10 ) THEN   ;   z1_ai(ji) = 1._wp / a_i_2d(ji,jl1)
                ELSE                                 ;   z1_ai(ji) = 0._wp
                ENDIF
@@ -592,5 +586,5 @@
                ! virtual salt flux to keep salinity constant
                IF( nn_icesal /= 2 )  THEN
-                  sirdg2(ji)     = sirdg2(ji)     - vsw * ( sss_1d(ji) - s_i_1d(ji) )        ! ridge salinity = s_i
+                  sirdg2(ji)     = sirdg2(ji)     - vsw * ( sss_1d(ji) - s_i_1d(ji) )       ! ridge salinity = s_i
                   sfx_bri_1d(ji) = sfx_bri_1d(ji) + sss_1d(ji) * vsw * rhoi * r1_rdtice  &  ! put back sss_m into the ocean
                      &                            - s_i_1d(ji) * vsw * rhoi * r1_rdtice     ! and get  s_i  from the ocean 
@@ -916,8 +910,8 @@
       REWIND( numnam_ice_ref )              ! Namelist namicetdme in reference namelist : Ice mechanical ice redistribution
       READ  ( numnam_ice_ref, namdyn_rdgrft, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_rdgrft in reference namelist', lwp )
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_rdgrft in reference namelist' )
       REWIND( numnam_ice_cfg )              ! Namelist namdyn_rdgrft in configuration namelist : Ice mechanical ice redistribution
       READ  ( numnam_ice_cfg, namdyn_rdgrft, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_rdgrft in configuration namelist', lwp )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_rdgrft in configuration namelist' )
       IF(lwm) WRITE ( numoni, namdyn_rdgrft )
       !
