Index: /NEMO/branches/2020/SI3_martin_ponds/cfgs/SHARED/field_def_nemo-ice.xml
===================================================================
--- /NEMO/branches/2020/SI3_martin_ponds/cfgs/SHARED/field_def_nemo-ice.xml	(revision 13907)
+++ /NEMO/branches/2020/SI3_martin_ponds/cfgs/SHARED/field_def_nemo-ice.xml	(revision 13908)
@@ -51,4 +51,8 @@
           <field id="icehlid"      long_name="melt pond lid depth"                                     standard_name="sea_ice_meltpondlid_depth"                 unit="m" /> 
           <field id="icevlid"      long_name="melt pond lid volume"                                    standard_name="sea_ice_meltpondlid_volume"                unit="m" /> 
+          <field id="dvpn_mlt"     long_name="pond volume tendency due to surface melt"                standard_name="sea_ice_pondvolume_tendency_melt"          unit="cm/d" />
+          <field id="dvpn_lid"     long_name="pond volume tendency due to exchanges with lid"          standard_name="sea_ice_pondvolume_tendency_lids"          unit="cm/d" />
+          <field id="dvpn_rnf"     long_name="pond volume tendency due to runoff"                      standard_name="sea_ice_pondvolume_tendency_runoff"        unit="cm/d" />
+          <field id="dvpn_drn"     long_name="pond volume tendency due to drainage"                    standard_name="sea_ice_pondvolume_tendency_drainage"      unit="cm/d" />  
 	  
 	  <!-- heat -->
@@ -297,8 +301,9 @@
           <field id="snwtemp_cat"  long_name="Snow temperature per category"                     unit="degC"    detect_missing_value="true" />
           <field id="icettop_cat"  long_name="Ice/snow surface temperature per category"         unit="degC"    detect_missing_value="true" />
-          <field id="iceapnd_cat"  long_name="Ice melt pond concentration per category"          unit=""        /> 
+          <field id="iceapnd_cat"  long_name="Ice melt pond grid fraction per category"          unit=""        /> 
+          <field id="icevpnd_cat"  long_name="Ice melt pond volume per grid area per category"   unit="m"       />  
           <field id="icehpnd_cat"  long_name="Ice melt pond thickness per category"              unit="m"       detect_missing_value="true" /> 
           <field id="icehlid_cat"  long_name="Ice melt pond lid thickness per category"          unit="m"       detect_missing_value="true" /> 
-          <field id="iceafpnd_cat" long_name="Ice melt pond fraction per category"               unit=""        /> 
+          <field id="iceafpnd_cat" long_name="Ice melt pond ice fraction per category"           unit=""        /> 
           <field id="iceaepnd_cat" long_name="Ice melt pond effective fraction per category"     unit=""        /> 
           <field id="icemask_cat"  long_name="Fraction of time step with sea ice (per category)" unit=""        />
Index: /NEMO/branches/2020/SI3_martin_ponds/cfgs/SHARED/namelist_ice_ref
===================================================================
--- /NEMO/branches/2020/SI3_martin_ponds/cfgs/SHARED/namelist_ice_ref	(revision 13907)
+++ /NEMO/branches/2020/SI3_martin_ponds/cfgs/SHARED/namelist_ice_ref	(revision 13908)
@@ -195,7 +195,9 @@
 !------------------------------------------------------------------------------
    ln_pnd            = .true.         !  activate melt ponds or not
-      ln_pnd_LEV     = .true.         !  level ice melt ponds (from Flocco et al 2007,2010 & Holland et al 2012)
-         rn_apnd_min =   0.15         !     minimum ice fraction that contributes to melt pond. range: 0.0 -- 0.15 ??
-         rn_apnd_max =   0.85         !     maximum ice fraction that contributes to melt pond. range: 0.7 -- 0.85 ??
+      ln_pnd_TOPO    = .false.        !  topographic melt ponds
+      ln_pnd_LEV     = .true.         !  level ice melt ponds
+      rn_apnd_min =   0.15            !     minimum meltwater fraction contributing to pond growth (TOPO and LEV)
+      rn_apnd_max =   0.85            !     maximum meltwater fraction contributing to pond growth (TOPO and LEV)
+      rn_pnd_flush=   0.1             !     pond flushing efficiency (tuning parameter) (LEV)
       ln_pnd_CST     = .false.        !  constant  melt ponds
          rn_apnd     =   0.2          !     prescribed pond fraction, at Tsu=0 degC
Index: /NEMO/branches/2020/SI3_martin_ponds/src/ICE/ice.F90
===================================================================
--- /NEMO/branches/2020/SI3_martin_ponds/src/ICE/ice.F90	(revision 13907)
+++ /NEMO/branches/2020/SI3_martin_ponds/src/ICE/ice.F90	(revision 13908)
@@ -208,7 +208,9 @@
    !                                     !!** ice-ponds namelist (namthd_pnd)
    LOGICAL , PUBLIC ::   ln_pnd           !: Melt ponds (T) or not (F)
-   LOGICAL , PUBLIC ::   ln_pnd_LEV       !: Melt ponds scheme from Holland et al (2012), Flocco et al (2007, 2010)
-   REAL(wp), PUBLIC ::   rn_apnd_min      !: Minimum ice fraction that contributes to melt ponds
-   REAL(wp), PUBLIC ::   rn_apnd_max      !: Maximum ice fraction that contributes to melt ponds
+   LOGICAL , PUBLIC ::   ln_pnd_TOPO      !: Topographic Melt ponds scheme (Flocco et al 2007, 2010)
+   LOGICAL , PUBLIC ::   ln_pnd_LEV       !: Simple melt pond scheme
+   REAL(wp), PUBLIC ::   rn_apnd_min      !: Minimum fraction of melt water contributing to ponds
+   REAL(wp), PUBLIC ::   rn_apnd_max      !: Maximum fraction of melt water contributing to ponds
+   REAL(wp), PUBLIC ::   rn_pnd_flush     !: Pond flushing efficiency (tuning parameter)
    LOGICAL , PUBLIC ::   ln_pnd_CST       !: Melt ponds scheme with constant fraction and depth
    REAL(wp), PUBLIC ::   rn_apnd          !: prescribed pond fraction (0<rn_apnd<1)
@@ -308,4 +310,8 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   t1_ice          !: temperature of the first layer          (ln_cndflx=T) [K]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   cnd_ice         !: effective conductivity of the 1st layer (ln_cndflx=T) [W.m-2.K-1]
+
+   ! meltwater arrays to save for melt ponds (mv - could be grouped in a single meltwater volume array)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dh_i_sum_2d   !: surface melt (2d arrays for ponds)       [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dh_s_mlt_2d   !: snow surf melt (2d arrays for ponds)     [m]
 
    !!----------------------------------------------------------------------
@@ -458,4 +464,5 @@
       ii = ii + 1
       ALLOCATE( qtr_ice_bot(jpi,jpj,jpl) , cnd_ice(jpi,jpj,jpl) , t1_ice(jpi,jpj,jpl) ,  &
+         &      dh_i_sum_2d(jpi,jpj,jpl) , dh_s_mlt_2d(jpi,jpj,jpl) ,                    &
          &      h_i        (jpi,jpj,jpl) , a_i    (jpi,jpj,jpl) , v_i   (jpi,jpj,jpl) ,  &
          &      v_s        (jpi,jpj,jpl) , h_s    (jpi,jpj,jpl) , t_su  (jpi,jpj,jpl) ,  &
Index: /NEMO/branches/2020/SI3_martin_ponds/src/ICE/icedyn_adv_pra.F90
===================================================================
--- /NEMO/branches/2020/SI3_martin_ponds/src/ICE/icedyn_adv_pra.F90	(revision 13907)
+++ /NEMO/branches/2020/SI3_martin_ponds/src/ICE/icedyn_adv_pra.F90	(revision 13908)
@@ -178,5 +178,5 @@
                z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content
             END DO
-            IF ( ln_pnd_LEV ) THEN
+            IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
                z0ap(:,:,jl) = pa_ip(:,:,jl) * e1e2t(:,:)      ! Melt pond fraction
                z0vp(:,:,jl) = pv_ip(:,:,jl) * e1e2t(:,:)      ! Melt pond volume
@@ -214,5 +214,5 @@
             END DO
             !
-            IF ( ln_pnd_LEV ) THEN
+            IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
                CALL adv_x( zdt , zudy , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction
                CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) 
@@ -249,5 +249,5 @@
                   &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
             END DO
-            IF ( ln_pnd_LEV ) THEN
+            IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
                CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction
                CALL adv_x( zdt , zudy , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )
@@ -278,5 +278,5 @@
          CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ei  , 'T', 1._wp, sxe   , 'T', -1._wp, sye   , 'T', -1._wp  & ! ice enthalpy
             &                                , sxxe  , 'T', 1._wp, syye  , 'T',  1._wp, sxye  , 'T',  1._wp  )
-         IF ( ln_pnd_LEV ) THEN
+         IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
             CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ap , 'T', 1._wp, sxap , 'T', -1._wp, syap , 'T', -1._wp  & ! melt pond fraction
                &                                , sxxap, 'T', 1._wp, syyap, 'T',  1._wp, sxyap, 'T',  1._wp  &
@@ -302,5 +302,5 @@
                pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
             END DO
-            IF ( ln_pnd_LEV ) THEN
+            IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
                pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
                pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
@@ -769,5 +769,5 @@
                !                               ! -- check h_ip -- !
                ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip
-               IF( ln_pnd_LEV .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN
+               IF( ln_pnd_LEV .OR. ln_pnd_TOPO .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN
                   zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) )
                   IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN
@@ -1005,5 +1005,5 @@
             END DO
             !
-            IF( ln_pnd_LEV ) THEN                                    ! melt pond fraction
+            IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN                                    ! melt pond fraction
                IF( iom_varid( numrir, 'sxap', ldstop = .FALSE. ) > 0 ) THEN
                   CALL iom_get( numrir, jpdom_auto, 'sxap' , sxap , psgn = -1._wp )
@@ -1047,5 +1047,5 @@
             sxc0  = 0._wp   ;   syc0  = 0._wp   ;   sxxc0  = 0._wp   ;   syyc0  = 0._wp   ;   sxyc0  = 0._wp      ! snow layers heat content
             sxe   = 0._wp   ;   sye   = 0._wp   ;   sxxe   = 0._wp   ;   syye   = 0._wp   ;   sxye   = 0._wp      ! ice layers heat content
-            IF( ln_pnd_LEV ) THEN
+            IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
                sxap = 0._wp ;   syap = 0._wp    ;   sxxap = 0._wp    ;   syyap = 0._wp    ;   sxyap = 0._wp       ! melt pond fraction
                sxvp = 0._wp ;   syvp = 0._wp    ;   sxxvp = 0._wp    ;   syyvp = 0._wp    ;   sxyvp = 0._wp       ! melt pond volume
@@ -1115,5 +1115,5 @@
          END DO
          !
-         IF( ln_pnd_LEV ) THEN                                       ! melt pond fraction
+         IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN                                       ! melt pond fraction
             CALL iom_rstput( iter, nitrst, numriw, 'sxap' , sxap  )
             CALL iom_rstput( iter, nitrst, numriw, 'syap' , syap  )
Index: /NEMO/branches/2020/SI3_martin_ponds/src/ICE/icedyn_adv_umx.F90
===================================================================
--- /NEMO/branches/2020/SI3_martin_ponds/src/ICE/icedyn_adv_umx.F90	(revision 13907)
+++ /NEMO/branches/2020/SI3_martin_ponds/src/ICE/icedyn_adv_umx.F90	(revision 13908)
@@ -338,5 +338,5 @@
          !
          !== melt ponds ==!
-         IF ( ln_pnd_LEV ) THEN
+         IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
             ! concentration
             zamsk = 1._wp
@@ -358,8 +358,8 @@
 
          ! --- Lateral boundary conditions --- !
-         IF    ( ln_pnd_LEV .AND. ln_pnd_lids ) THEN
+         IF    ( ( ln_pnd_LEV .OR. ln_pnd_TOPO ) .AND. ln_pnd_lids ) THEN
             CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp &
                &                                , pa_ip,'T',1._wp, pv_ip,'T',1._wp, pv_il,'T',1._wp )
-         ELSEIF( ln_pnd_LEV .AND. .NOT.ln_pnd_lids ) THEN
+         ELSEIF( ln_pnd_LEV .OR. ln_pnd_TOPO ) .AND. .NOT.ln_pnd_lids ) THEN
             CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp &
                &                                , pa_ip,'T',1._wp, pv_ip,'T',1._wp )
@@ -1497,5 +1497,5 @@
                !                               ! -- check h_ip -- !
                ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip
-               IF( ln_pnd_LEV .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN
+               IF( ( ln_pnd_LEV .OR. ln_pnd_TOPO ) .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN
                   zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) )
                   IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN
Index: /NEMO/branches/2020/SI3_martin_ponds/src/ICE/icedyn_rdgrft.F90
===================================================================
--- /NEMO/branches/2020/SI3_martin_ponds/src/ICE/icedyn_rdgrft.F90	(revision 13907)
+++ /NEMO/branches/2020/SI3_martin_ponds/src/ICE/icedyn_rdgrft.F90	(revision 13908)
@@ -575,5 +575,5 @@
                oirft2(ji) = oa_i_2d(ji,jl1)   * afrft * hi_hrft 
 
-               IF ( ln_pnd_LEV ) THEN
+               IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
                   aprdg1     = a_ip_2d(ji,jl1) * afrdg
                   aprdg2(ji) = a_ip_2d(ji,jl1) * afrdg * hi_hrdg(ji,jl1)
@@ -612,5 +612,5 @@
                sv_i_2d(ji,jl1) = sv_i_2d(ji,jl1) - sirdg1    - sirft(ji)
                oa_i_2d(ji,jl1) = oa_i_2d(ji,jl1) - oirdg1    - oirft1
-               IF ( ln_pnd_LEV ) THEN
+               IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
                   a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - aprdg1    - aprft1
                   v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - vprdg(ji) - vprft(ji)
@@ -709,5 +709,5 @@
                   v_s_2d (ji,jl2) = v_s_2d (ji,jl2) + ( vsrdg (ji) * rn_fsnwrdg * fvol(ji)  +  &
                      &                                  vsrft (ji) * rn_fsnwrft * zswitch(ji) )
-                  IF ( ln_pnd_LEV ) THEN
+                  IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
                      v_ip_2d (ji,jl2) = v_ip_2d(ji,jl2) + (   vprdg (ji) * rn_fpndrdg * fvol   (ji)   &
                         &                                   + vprft (ji) * rn_fpndrft * zswitch(ji)   )
Index: /NEMO/branches/2020/SI3_martin_ponds/src/ICE/iceitd.F90
===================================================================
--- /NEMO/branches/2020/SI3_martin_ponds/src/ICE/iceitd.F90	(revision 13907)
+++ /NEMO/branches/2020/SI3_martin_ponds/src/ICE/iceitd.F90	(revision 13908)
@@ -315,5 +315,5 @@
             IF ( a_i_1d(ji) > epsi10 .AND. h_i_1d(ji) < rn_himin ) THEN
                a_i_1d(ji) = a_i_1d(ji) * h_i_1d(ji) / rn_himin 
-               IF( ln_pnd_LEV )   a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin
+               IF( ln_pnd_LEV .OR. ln_pnd_TOPO )   a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin
                h_i_1d(ji) = rn_himin
             ENDIF
@@ -486,5 +486,5 @@
                zaTsfn(ji,jl2)  = zaTsfn(ji,jl2) + ztrans
                !  
-               IF ( ln_pnd_LEV ) THEN
+               IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
                   ztrans          = a_ip_2d(ji,jl1) * zworka(ji)     ! Pond fraction
                   a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - ztrans
Index: /NEMO/branches/2020/SI3_martin_ponds/src/ICE/icestp.F90
===================================================================
--- /NEMO/branches/2020/SI3_martin_ponds/src/ICE/icestp.F90	(revision 13907)
+++ /NEMO/branches/2020/SI3_martin_ponds/src/ICE/icestp.F90	(revision 13908)
@@ -454,4 +454,7 @@
             qcn_ice    (ji,jj,jl) = 0._wp   ! conductive flux (ln_cndflx=T & ln_cndemule=T)
             qtr_ice_bot(ji,jj,jl) = 0._wp   ! part of solar radiation transmitted through the ice needed at least for outputs
+            ! Melt pond surface melt diagnostics (mv - more efficient: grouped into one water volume flux)
+            dh_i_sum_2d(ji,jj,jl) = 0._wp
+            dh_s_mlt_2d(ji,jj,jl) = 0._wp
          END_2D
       ENDDO
Index: /NEMO/branches/2020/SI3_martin_ponds/src/ICE/icethd.F90
===================================================================
--- /NEMO/branches/2020/SI3_martin_ponds/src/ICE/icethd.F90	(revision 13907)
+++ /NEMO/branches/2020/SI3_martin_ponds/src/ICE/icethd.F90	(revision 13908)
@@ -239,5 +239,4 @@
             IF( ln_icedH ) THEN                                     ! --- Growing/Melting --- !
                               CALL ice_thd_dh                           ! Ice-Snow thickness   
-                              CALL ice_thd_pnd                          ! Melt ponds formation
                               CALL ice_thd_ent( e_i_1d(1:npti,:) )      ! Ice enthalpy remapping
             ENDIF
@@ -260,4 +259,6 @@
       IF( ln_icediachk )   CALL ice_cons2D  (1, 'icethd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft)
       !                   
+      IF ( ln_pnd  )          CALL ice_thd_pnd                      ! --- Melt ponds 
+      !
       IF( jpl > 1  )          CALL ice_itd_rem( kt )                ! --- Transport ice between thickness categories --- !
       !
@@ -529,4 +530,7 @@
          CALL tab_1d_2d( npti, nptidx(1:npti), cnd_ice_1d(1:npti), cnd_ice(:,:,kl) )
          CALL tab_1d_2d( npti, nptidx(1:npti), t1_ice_1d (1:npti), t1_ice (:,:,kl) )
+         ! Melt ponds
+         CALL tab_1d_2d( npti, nptidx(1:npti), dh_i_sum  (1:npti) , dh_i_sum_2d(:,:,kl) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), dh_s_mlt  (1:npti) , dh_s_mlt_2d(:,:,kl) )
          ! SIMIP diagnostics         
          CALL tab_1d_2d( npti, nptidx(1:npti), t_si_1d       (1:npti), t_si       (:,:,kl) )
Index: /NEMO/branches/2020/SI3_martin_ponds/src/ICE/icevar.F90
===================================================================
--- /NEMO/branches/2020/SI3_martin_ponds/src/ICE/icevar.F90	(revision 13907)
+++ /NEMO/branches/2020/SI3_martin_ponds/src/ICE/icevar.F90	(revision 13908)
@@ -554,4 +554,19 @@
          END_2D
          !
+         !-----------------------------------------------------------------
+         ! zap small ponds
+         !-----------------------------------------------------------------
+         IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
+            DO_2D( 1, 1, 1, 1 )
+               IF ( v_ip(ji,jj,jl) <= epsi10 ) THEN
+                  a_ip(ji,jj,jl)      = 0._wp
+                  a_ip_frac(ji,jj,jl) = 0._wp
+                  v_ip(ji,jj,jl)      = 0._wp
+                  h_ip(ji,jj,jl)      = 0._wp
+                  v_il(ji,jj,jl)      = 0._wp
+                  h_il(ji,jj,jl)      = 0._wp
+               ENDIF
+            END_2D
+         ENDIF
       END DO 
 
@@ -675,5 +690,5 @@
       WHERE( pe_i (1:npti,:,:) < 0._wp )   pe_i (1:npti,:,:) = 0._wp   !  e_i must be >= 0
       WHERE( pe_s (1:npti,:,:) < 0._wp )   pe_s (1:npti,:,:) = 0._wp   !  e_s must be >= 0
-      IF( ln_pnd_LEV ) THEN
+      IF( ln_pnd_LEV .OR ln_pnd_TOPO ) THEN
          WHERE( pa_ip(1:npti,:) < 0._wp )    pa_ip(1:npti,:)   = 0._wp   ! a_ip must be >= 0
          WHERE( pv_ip(1:npti,:) < 0._wp )    pv_ip(1:npti,:)   = 0._wp   ! v_ip must be >= 0
Index: /NEMO/branches/2020/SI3_martin_ponds/src/ICE/icewri.F90
===================================================================
--- /NEMO/branches/2020/SI3_martin_ponds/src/ICE/icewri.F90	(revision 13907)
+++ /NEMO/branches/2020/SI3_martin_ponds/src/ICE/icewri.F90	(revision 13908)
@@ -160,7 +160,8 @@
       IF( iom_use('icebrv_cat'  ) )   CALL iom_put( 'icebrv_cat'  ,   bv_i * 100.  * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! brine volume
       IF( iom_use('iceapnd_cat' ) )   CALL iom_put( 'iceapnd_cat' ,   a_ip         * zmsk00l                                   ) ! melt pond frac for categories
+      IF( iom_use('icevpnd_cat' ) )   CALL iom_put( 'icevpnd_cat' ,   v_ip         * zmsk00l                                   ) ! melt pond volume for categories
       IF( iom_use('icehpnd_cat' ) )   CALL iom_put( 'icehpnd_cat' ,   h_ip         * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! melt pond thickness for categories
       IF( iom_use('icehlid_cat' ) )   CALL iom_put( 'icehlid_cat' ,   h_il         * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! melt pond lid thickness for categories
-      IF( iom_use('iceafpnd_cat') )   CALL iom_put( 'iceafpnd_cat',   a_ip_frac    * zmsk00l                                   ) ! melt pond frac for categories
+      IF( iom_use('iceafpnd_cat') )   CALL iom_put( 'iceafpnd_cat',   a_ip_frac    * zmsk00l                                   ) ! melt pond frac per ice area for categories
       IF( iom_use('iceaepnd_cat') )   CALL iom_put( 'iceaepnd_cat',   a_ip_eff     * zmsk00l                                   ) ! melt pond effective frac for categories
       IF( iom_use('icealb_cat'  ) )   CALL iom_put( 'icealb_cat'  ,   alb_ice      * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! ice albedo for categories
