Index: NEMO/branches/UKMO/NEMO_4.0.4_ice_strength_EAP/src/ICE/icedyn_rdgrft.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.4_ice_strength_EAP/src/ICE/icedyn_rdgrft.F90	(revision 15647)
+++ NEMO/branches/UKMO/NEMO_4.0.4_ice_strength_EAP/src/ICE/icedyn_rdgrft.F90	(revision 15734)
@@ -56,4 +56,10 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ze_i_2d
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ze_s_2d
+   !
+   ! For ridging diagnostics
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)     ::   airdg1          ! Ridging ice area loss
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)     ::   airft1          ! Rafting ice area loss
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)     ::   airdg2          ! New ridged ice area gain
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)     ::   airft2          ! New rafted ice area gain
    !
    REAL(wp), PARAMETER ::   hrdg_hi_min = 1.1_wp    ! min ridge thickness multiplier: min(hrdg/hi)
@@ -100,5 +106,6 @@
          &      zaksum(jpij)       , hraft(jpij,jpl)   , hrexp(jpij,jpl)     ,  &
          &      hi_hrdg(jpij,jpl)  , araft(jpij,jpl)   , aridge(jpij,jpl)    , zasum(jpij), &
-         &      ze_i_2d(jpij,nlay_i,jpl), ze_s_2d(jpij,nlay_s,jpl), STAT=ice_dyn_rdgrft_alloc )
+         &      ze_i_2d(jpij,nlay_i,jpl), ze_s_2d(jpij,nlay_s,jpl), &
+         &      airdg1(jpij), airft1(jpij), airdg2(jpij), airft2(jpij), STAT=ice_dyn_rdgrft_alloc )
 
       CALL mpp_sum ( 'icedyn_rdgrft', ice_dyn_rdgrft_alloc )
@@ -150,4 +157,15 @@
       REAL(wp), DIMENSION(jpij) ::   zdivu, zdelt  ! 1D divu_i & delta_i
       REAL(wp), DIMENSION(jpij) ::   zconv         ! 1D rdg_conv (if EAP rheology)
+
+      REAL(wp), DIMENSION(jpi,jpj) ::   opning_2d  ! Lead opening diagnostic
+      REAL(wp), DIMENSION(jpi,jpj) ::   airdg1_2d  ! Ridging ice area loss diagnostic
+      REAL(wp), DIMENSION(jpi,jpj) ::   airft1_2d  ! Rafting ice area loss diagnostic
+      REAL(wp), DIMENSION(jpi,jpj) ::   airdg2_2d  ! New ridged ice area gain diagnostic
+      REAL(wp), DIMENSION(jpi,jpj) ::   airft2_2d  ! New rafted ice area gain diagnostic
+      REAL(wp), DIMENSION(jpi,jpj) ::   dairdg1dt  ! Ridging ice area loss rate diagnostic
+      REAL(wp), DIMENSION(jpi,jpj) ::   dairft1dt  ! Rafting ice area loss rate diagnostic
+      REAL(wp), DIMENSION(jpi,jpj) ::   dairdg2dt  ! New ridged ice area gain rate diagnostic
+      REAL(wp), DIMENSION(jpi,jpj) ::   dairft2dt  ! New rafted ice area gain rate diagnostic
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00     ! Mask for ice presence
       !
       INTEGER, PARAMETER ::   jp_itermax = 20    
@@ -162,5 +180,16 @@
          IF(lwp) WRITE(numout,*)'ice_dyn_rdgrft: ice ridging and rafting'
          IF(lwp) WRITE(numout,*)'~~~~~~~~~~~~~~'
-      ENDIF      
+      ENDIF     
+
+      ! Initialise ridging diagnostics
+      opning_2d(:,:) = 0.0_wp
+      dairdg1dt(:,:) = 0.0_wp
+      dairft1dt(:,:) = 0.0_wp
+      dairdg2dt(:,:) = 0.0_wp
+      dairft2dt(:,:) = 0.0_wp
+      airdg1_2d(:,:) = 0.0_wp
+      airft1_2d(:,:) = 0.0_wp
+      airdg2_2d(:,:) = 0.0_wp
+      airft2_2d(:,:) = 0.0_wp
 
       !--------------------------------
@@ -173,4 +202,5 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
+            zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06  ) ) ! 1 if ice, 0 if no ice. epsi06 used here rather than epsi10 as below for consistency with ice masking elsewhere
             IF ( at_i(ji,jj) > epsi10 ) THEN
                npti           = npti + 1
@@ -278,6 +308,60 @@
          CALL ice_dyn_1d2d( 2 )            ! --- Move to 2D arrays --- !
 
+! --- Ridging diagnostics --- !
+
+         IF( iom_use('opning') ) THEN
+           CALL tab_1d_2d( npti, nptidx(1:npti), opning(1:npti), opning_2d(:,:) )
+         END IF
+
+         IF( iom_use('dairdg1dt') .OR. iom_use('dairft1dt') .OR. iom_use('dairdg2dt') .OR. iom_use('dairft2dt') ) THEN
+           !
+           CALL tab_1d_2d( npti, nptidx(1:npti), airdg1(1:npti), airdg1_2d(:,:) )
+           CALL tab_1d_2d( npti, nptidx(1:npti), airft1(1:npti), airft1_2d(:,:) )
+           CALL tab_1d_2d( npti, nptidx(1:npti), airdg2(1:npti), airdg2_2d(:,:) )
+           CALL tab_1d_2d( npti, nptidx(1:npti), airft2(1:npti), airft2_2d(:,:) )
+           !
+           DO jj = 1, jpj
+              DO ji = 1, jpi
+              !
+              dairdg1dt(ji,jj) = airdg1_2d(ji,jj) * r1_rdtice
+              dairft1dt(ji,jj) = airft1_2d(ji,jj) * r1_rdtice                      
+              dairdg2dt(ji,jj) = airdg2_2d(ji,jj) * r1_rdtice
+              dairft2dt(ji,jj) = airft2_2d(ji,jj) * r1_rdtice 
+              !
+              END DO
+           END DO    
+          !         
+         ENDIF
+        ! 
+      ENDIF  ! npti>0
+
+
+      IF( iom_use('opning') ) THEN
+        CALL lbc_lnk( 'icedyn_rdgrft', opning_2d(:,:), 'T', 1. )
+        CALL iom_put( 'opning', opning_2d(:,:) * zmsk00(:,:) )     ! Lead area opening rate
       ENDIF
-   
+
+      IF( iom_use('dairdg1dt') ) THEN
+        CALL lbc_lnk( 'icedyn_rdgrft', dairdg1dt(:,:), 'T', 1. )
+        CALL iom_put( 'dairdg1dt', dairdg1dt(:,:) * zmsk00(:,:) )  ! Ridging ice area loss rate
+      ENDIF
+
+      IF( iom_use('dairft1dt') ) THEN
+        CALL lbc_lnk( 'icedyn_rdgrft', dairft1dt(:,:), 'T', 1. )
+        CALL iom_put( 'dairft1dt', dairft1dt(:,:) * zmsk00(:,:) )  ! Rafting ice area loss rate
+      ENDIF
+
+      IF( iom_use('dairdg2dt') ) THEN
+        CALL lbc_lnk( 'icedyn_rdgrft', dairdg2dt(:,:), 'T', 1. )
+        CALL iom_put( 'dairdg2dt', dairdg2dt(:,:) * zmsk00(:,:) )  ! New ridged ice area gain rate
+      ENDIF
+
+      IF( iom_use('dairft2dt') ) THEN
+        CALL lbc_lnk( 'icedyn_rdgrft', dairft2dt(:,:), 'T', 1. )
+        CALL iom_put( 'dairft2dt', dairft2dt(:,:) * zmsk00(:,:) )  ! New rafted ice area gain rate
+      ENDIF
+
+! ------- !
+
       CALL ice_var_agg( 1 ) 
 
@@ -533,8 +617,8 @@
       REAL(wp) ::   vsw                        ! vol of water trapped into ridges
       REAL(wp) ::   afrdg, afrft               ! fraction of category area ridged/rafted 
-      REAL(wp)                  ::   airdg1, oirdg1, aprdg1, virdg1, sirdg1
-      REAL(wp)                  ::   airft1, oirft1, aprft1
-      REAL(wp), DIMENSION(jpij) ::   airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg, vlrdg  ! area etc of new ridges
-      REAL(wp), DIMENSION(jpij) ::   airft2, oirft2, aprft2, virft , sirft , vsrft, vprft, vlrft  ! area etc of rafted ice
+      REAL(wp)                  ::   oirdg1, aprdg1, virdg1, sirdg1
+      REAL(wp)                  ::   oirft1, aprft1
+      REAL(wp), DIMENSION(jpij) ::   oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg, vlrdg  ! area etc of new ridges
+      REAL(wp), DIMENSION(jpij) ::   oirft2, aprft2, virft , sirft , vsrft, vprft, vlrft  ! area etc of rafted ice
       !
       REAL(wp), DIMENSION(jpij) ::   ersw             ! enth of water trapped into ridges
@@ -576,13 +660,13 @@
                
                ! area of ridging / rafting ice (airdg1) and of new ridge (airdg2)
-               airdg1 = aridge(ji,jl1) * closing_gross(ji) * rdt_ice
-               airft1 = araft (ji,jl1) * closing_gross(ji) * rdt_ice
-
-               airdg2(ji) = airdg1 * hi_hrdg(ji,jl1)
-               airft2(ji) = airft1 * hi_hrft
+               airdg1(ji) = aridge(ji,jl1) * closing_gross(ji) * rdt_ice
+               airft1(ji) = araft (ji,jl1) * closing_gross(ji) * rdt_ice
+
+               airdg2(ji) = airdg1(ji) * hi_hrdg(ji,jl1)
+               airft2(ji) = airft1(ji) * hi_hrft
 
                ! ridging /rafting fractions
-               afrdg = airdg1 * z1_ai(ji)
-               afrft = airft1 * z1_ai(ji)
+               afrdg = airdg1(ji) * z1_ai(ji)
+               afrft = airft1(ji) * z1_ai(ji)
 
                ! volume and enthalpy (J/m2, >0) of seawater trapped into ridges
@@ -640,5 +724,5 @@
                ! Remove area, volume of new ridge to each category jl1
                !------------------------------------------------------
-               a_i_2d (ji,jl1) = a_i_2d (ji,jl1) - airdg1    - airft1
+               a_i_2d (ji,jl1) = a_i_2d (ji,jl1) - airdg1(ji)- airft1(ji)
                v_i_2d (ji,jl1) = v_i_2d (ji,jl1) - virdg1    - virft(ji)
                v_s_2d (ji,jl1) = v_s_2d (ji,jl1) - vsrdg(ji) - vsrft(ji)
