Index: NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE/icethd_pnd.F90
===================================================================
--- NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE/icethd_pnd.F90	(revision 13890)
+++ NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE/icethd_pnd.F90	(revision 13891)
@@ -45,4 +45,31 @@
       zTd     = 0.15_wp       , & ! temperature difference for freeze-up (C)
       zvp_min = 1.e-4_wp          ! minimum pond volume (m)
+   ! 
+   ! Pond volume per area budget diags
+   !  
+   ! The idea of diags is the volume of ponds per grid cell area is
+   !
+   ! dV/dt = mlt + drn + lid + rnf
+   ! mlt   = input from surface melting
+   ! drn   = drainage through brine network
+   ! lid   = lid growth & melt
+   ! rnf   = runoff (water directly removed out of surface melting + overflow)
+   !
+   ! In topo mode, the pond water lost because it is in the snow is not included in the budget
+   !
+   ! In level mode
+
+   REAL(wp), DIMENSION(jpi,jpj) ::   ! pond volume budget diagnostics
+      diag_dvpn_mlt,  &     !! meltwater pond volume input      [m/day]
+      diag_dvpn_drn,  &     !! pond volume lost by drainage     [m/day]
+      diag_dvpn_lid,  &     !! pond volume lost by refreezing   [m/day]
+      diag_dvpn_rnf         !! meltwater pond lost to runoff    [m/day]
+      
+   REAL(wp), DIMENSION(jpij) ::   ! pond volume budget diagnostics (1d)
+      diag_dvpn_mlt_1d,  &  !! meltwater pond volume input      [m/day]
+      diag_dvpn_drn_1d,  &  !! pond volume lost by drainage     [m/day]
+      diag_dvpn_lid_1d,  &  !! pond volume lost by refreezing   [m/day]
+      diag_dvpn_rnf_1d      !! meltwater pond lost to runoff    [m/day]
+
 
    !! * Substitutions
@@ -76,4 +103,9 @@
       !!-------------------------------------------------------------------
       !
+      diag_dvpn_mlt(:,:) = 0._wp      ;   diag_dvpn_drn(:,:)    = 0._wp
+      diag_dvpn_lid(:,:) = 0._wp      ;   diag_dvpn_rnf(:,:)    = 0._wp
+      diag_dvpn_mlt_1d(:,:) = 0._wp   ;   diag_dvpn_drn_1d(:,:) = 0._wp
+      diag_dvpn_lid_1d(:,:) = 0._wp   ;   diag_dvpn_rnf_1d(:,:) = 0._wp
+
       SELECT CASE ( nice_pnd )
       !
@@ -82,9 +114,16 @@
       CASE (np_pndLEV)   ;   CALL pnd_LEV                              !==  Level ice melt ponds  ==!
          !
-      CASE (np_pndTOPO)  ;   CALL pnd_TOPO                  !&         !==  Topographic melt ponds  ==!
+      CASE (np_pndTOPO)  ;   CALL pnd_TOPO                             !==  Topographic melt ponds  ==!
          !
       END SELECT
       !
+
+      IF( iom_use('dvpn_mlt'  ) )   CALL iom_put( 'dvpn_mlt', dvpn_mlt * 100._wp ) ! input from melting
+      IF( iom_use('dvpn_lid'  ) )   CALL iom_put( 'dvpn_lid', dvpn_lid * 100._wp ) ! exchanges with lid
+      IF( iom_use('dvpn_drn'  ) )   CALL iom_put( 'dvpn_drn', dvpn_drn * 100._wp ) ! vertical drainage
+      IF( iom_use('dvpn_rnf'  ) )   CALL iom_put( 'dvpn_rnf', dvpn_rnf * 100._wp ) ! runoff + overflow
+      
    END SUBROUTINE ice_thd_pnd 
+
 
 
@@ -191,4 +230,6 @@
       REAL(wp) ::   zfac, zdum                        ! temporary arrays
       REAL(wp) ::   z1_rhow, z1_aspect, z1_Tp         ! inverse
+      REAL(wp) ::   z1_rdtice                         ! inverse time step
+      REAL(wp) ::   zvold                             !
       !!
       INTEGER  ::   ji, jj, jk, jl                    ! loop indices
@@ -199,4 +240,5 @@
       z1_aspect = 1._wp / zaspect
       z1_Tp     = 1._wp / zTp 
+      z1_rdtice = 1._wp / rdt_ice
       
       !-----------------------------------------------------------------------------------------------
@@ -221,7 +263,12 @@
       IF( npti > 0 ) THEN
       
-         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum   )
-         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sum_1d (1:npti), wfx_sum          )
-         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd          )
+         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum      )
+         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sum_1d (1:npti)   , wfx_sum          )
+         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti)   , wfx_pnd          )
+         
+         CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt )
+         CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn )
+         CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_lid_1d (1:npti), diag_dvpn_lid )
+         CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_rnf_1d (1:npti), diag_dvpn_rnf )
       
          DO jl = 1, jpl
@@ -268,17 +315,32 @@
                   zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i_1d(ji) !  = ( 1 - r ) = fraction of melt water that is not flushed
                   zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! max for roundoff errors? 
+                  
+                  diag_dvpn_mlt_1d(ji) = diag_dvpn_mlt_1d(ji) + zdum * z1_rdtice                      ! surface melt input diag
+                  
+                  diag_dvpn_rnf_1d(ji) = diag_dvpn_rnf_1d(ji) + ( 1. - zfr_mlt ) * zdum * z1_rdtice   ! runoff diag
+                  
                   !
                   !--- overflow ---!
-                  ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume
+                  !
+                  ! 1) area driven overflow
+                  !
+                  ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond water volume
                   !    a_ip_max = zfr_mlt * a_i
                   !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 
                   zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect
-                  zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) )
-
+                  zvold     = zdv_mlt
+                  zdv_mlt   = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) )                     
+                  
+                  ! 
+                  ! 2) depth driven overflow
+                  !
                   ! If pond depth exceeds half the ice thickness then reduce the pond volume
                   !    h_ip_max = 0.5 * h_i
                   !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 
                   zv_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji) ! MV dimensions are wrong here or comment is unclear
+                                   
                   zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) )
+                  
+                  diag_dvpn_rnf_1d(ji) = diag_dvpn_rnf_1d(ji) + ( zdv_mlt - zvold ) * z1_rdtice       ! runoff diag - overflow contribution
             
                   !--- Pond growing ---!
@@ -308,4 +370,6 @@
                   !
                   !--- Pond contraction (due to refreezing) ---!
+                  zvold       = v_ip_1d(ji) ! for diag
+
                   IF( ln_pnd_lids ) THEN
                      !
@@ -319,9 +383,14 @@
                      ! Pond shrinking
                      v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) - zdv_frz )
-
+                     
                   ELSE
+                  
                      ! Pond shrinking
                      v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * zdT * z1_Tp ) ! Holland 2012 (eq. 6)
+                     
                   ENDIF
+
+                  diag_dvpn_lid_1d(ji) = diag_dvpn_lid_1d(ji) + ( v_ip_1d(ji) - zvold ) * z1_rdtice   ! shrinking counted as lid diagnostic
+
                   !
                   !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac
@@ -331,7 +400,7 @@
                   h_ip_1d(ji)      = zaspect * a_ip_1d(ji) / a_i_1d(ji)
 
-                  !---------------!            
-                  ! Pond flushing !
-                  !---------------!
+                  !------------------------------------------------!            
+                  ! Pond drainage through brine network (flushing) !
+                  !------------------------------------------------!
                   ! height of top of the pond above sea-level
                   zhp = ( h_i_1d(ji) * ( rau0 - rhoi ) + h_ip_1d(ji) * ( rau0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rau0
@@ -352,9 +421,10 @@
                   zdv_flush   = -zperm * rau0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji)
                   zdv_flush   = MAX( zdv_flush, -v_ip_1d(ji) )
-                  zdv_flush   = 0._wp ! MV remove pond drainage for now
+                  ! zdv_flush   = 0._wp ! MV remove pond drainage for now
                   v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush
+                  
+                  diag_dvpn_drn_1d(ji) = diag_dvpn_drn_1d(ji) + zdv_flush * z1_rdtice   ! shrinking counted as lid diagnostic
             
-                  ! MV --- why pond drainage does not give back water into freshwater flux ?
-            
+                  ! MV --- why pond drainage does not give back water into freshwater flux ?            
                   !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac
                   a_ip_1d(ji)      = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i
@@ -385,4 +455,5 @@
             v_ip_1d(1:npti) = h_ip_1d(1:npti) * a_ip_1d(1:npti)
             v_il_1d(1:npti) = h_il_1d(1:npti) * a_ip_1d(1:npti)
+            
             CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d     (1:npti), a_ip     (:,:,jl) )
             CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,jl) )
@@ -399,5 +470,10 @@
          CALL tab_1d_2d( npti, nptidx(1:npti), wfx_sum_1d (1:npti), wfx_sum        )
          CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd        )
-                  
+         
+         CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt        )
+         CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn        )
+         CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_lid_1d (1:npti), diag_dvpn_lid        )
+         CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_rnf_1d (1:npti), diag_dvpn_rnf        )                  
+                           
       !
       ENDIF
@@ -438,20 +514,21 @@
       ! local variables
       REAL(wp) :: &
-         zdHui,   &   ! change in thickness of ice lid (m)
-         zomega,  &   ! conduction
-         zdTice,  &   ! temperature difference across ice lid (C)
-         zdvice,  &   ! change in ice volume (m)
-         zTavg,   &   ! mean surface temperature across categories (C)
-         zfsurf,  &   ! net heat flux, excluding conduction and transmitted radiation (W/m2)
-         zTp,     &   ! pond freezing temperature (C)
-         zrhoi_L, &   ! volumetric latent heat of sea ice (J/m^3)
-         zfr_mlt, &   ! fraction and volume of available meltwater retained for melt ponding
-         z1_rhow, &   ! inverse water density
-         zpond  , &   ! dummy variable
-         zdum         ! dummy variable
-
-      REAL(wp), DIMENSION(jpi,jpj) ::   zvolp, & !! total melt pond water available before redistribution and drainage
-                                        zvolp_res
-
+         zdHui,   &      ! change in thickness of ice lid (m)
+         zomega,  &      ! conduction
+         zdTice,  &      ! temperature difference across ice lid (C)
+         zdvice,  &      ! change in ice volume (m)
+         zTavg,   &      ! mean surface temperature across categories (C)
+         zfsurf,  &      ! net heat flux, excluding conduction and transmitted radiation (W/m2)
+         zTp,     &      ! pond freezing temperature (C)
+         zrhoi_L, &      ! volumetric latent heat of sea ice (J/m^3)
+         zfr_mlt, &      ! fraction and volume of available meltwater retained for melt ponding
+         z1_rhow, &      ! inverse water density
+         z1_rdtice, &    ! inverse time step
+         zv_pnd  , &     ! volume of meltwater contributing to ponds
+         zv_mlt          ! total amount of meltwater produced
+
+      REAL(wp), DIMENSION(jpi,jpj) ::   zvolp, &     !! total melt pond water available before redistribution and drainage
+                                        zvolp_res    !! remaining melt pond water available after drainage
+                                        
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_a_i
 
@@ -473,4 +550,5 @@
       zTp       = rt0 - 0.15_wp          ! pond freezing point, slightly below 0C (ponds are bid saline)
       z1_rhow   = 1._wp / rhow 
+      z1_rdtice = 1._wp / rdt_ice
 
       ! Set required ice variables (hard-coded here for now)
@@ -491,9 +569,17 @@
       ! Change 2D to 1D
       !---------------------------------------------------------------
+      
+      ! use what we have in iceitd.F90 (incremental remapping)
 
       !--------------------------------------------------------------
-      ! Collect total available pond water
+      ! Collect total available pond water volume
       !--------------------------------------------------------------
-
+      ! Assuming that meltwater (+rain in principle) runsoff the surface
+      ! Holland et al (2012) suggest that the fraction of runoff decreases with total ice fraction
+      ! I cite her words, they are very talkative
+      ! "grid cells with very little ice cover (and hence more open water area) 
+      ! have a higher runoff fraction to rep- resent the greater proximity of ice to open water."
+      ! "This results in the same runoff fraction r for each ice category within a grid cell"
+      
       zvolp(:,:) = 0._wp
 
@@ -505,7 +591,13 @@
             
                   !--- Available meltwater for melt ponding ---!
-                  zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * a_i(ji,jj,jl) !  = ( 1 - r ) = fraction of melt water that is not flushed
-                  zdum    = -( dh_i_sum_2d(ji,jj,jl)*rhoi + dh_s_mlt_2d(ji,jj,jl)*rhos ) * z1_rhow * a_i(ji,jj,jl)
-                  zpond   = zfr_mlt * zdum ! check
+                  zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i(ji,jj)                  ! fraction of surface meltwater going to ponds
+
+                  zv_mlt  = - ( dh_i_sum_2d(ji,jj,jl) * rhoi + dh_s_mlt_2d(ji,jj,jl) * rhos ) &        ! total volume of surface melt water
+                     &    * z1_rhow * a_i(ji,jj,jl)
+                  zv_pnd  = zfr_mlt * zv_mlt                                                           ! volume of meltwater contributing to ponds
+                  
+                  diag_dvpn_mlt(ji,jj) = diag_dvpn_mlt(ji,jj) + zv_mlt * z1_rdtice                     ! diagnostics
+                  
+                  diag_dvpn_rnf(ji,jj) = diag_dvpn_rnf(ji,jj) + ( 1. - zfr_mlt ) * zv_mlt * z1_rdtice   
 
                   !--- Create possible new ponds
@@ -513,11 +605,17 @@
                   IF ( a_ip_frac(ji,jj,jl) < epsi10 ) THEN
                      h_ip(ji,jj,jl)       = 0._wp
-                     a_ip_frac(ji,jj,jl) = 1.0_wp    ! pond fraction of sea ice (apnd for CICE)
+                     a_ip_frac(ji,jj,jl)  = 1.0_wp    ! pond fraction of sea ice (apnd for CICE)
+                     a_ip(ji,jj,jl)       = a_i(ji,jj,jl)
                   ENDIF
                   
-                  !--- Deepen existing ponds before redistribution and drainage(later on)
-                  h_ip(ji,jj,jl) = ( zpond + h_ip(ji,jj,jl) * a_ip_frac(ji,jj,jl) ) / a_ip_frac(ji,jj,jl)
-                  zvolp(ji,jj)   = zvolp(ji,jj) + h_ip(ji,jj,jl) * a_ip_frac(ji,jj,jl) * a_i(ji,jj,jl) 
-!                 zfpond(ji,jj) = zfpond(ji,jj) + zpond * a_ip_frac(ji,jj,jl) !useless
+                  !--- Deepen existing ponds before redistribution and drainage
+                  ! without pond fraction
+                  v_ip(ji,jj,jl) = v_i_p(ji,jj,jl) +  zv_pnd                                            ! use pond water to increase thickness
+                     
+                  h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl)
+                  
+                  zvolp(ji,jj)   = zvolp(ji,jj) + v_ip(ji,jj,jl)
+                  
+!                 zfpond(ji,jj) = zfpond(ji,jj) + zpond * a_ip_frac(ji,jj,jl) ! useless for now
                    
                ENDIF ! a_i
@@ -526,8 +624,5 @@
          END DO ! jj
       END DO ! ji
-         
-      h_ip(:,:,:)  = 0._wp ! pond thickness
-      a_ip(:,:,:)  = 0._wp ! pond fraction per grid area  
-         
+                  
       !--------------------------------------------------------------
       ! Redistribute and drain water from ponds
@@ -564,7 +659,10 @@
 
                         zdvice = MIN( dh_i_sum_2d(ji,jj,jl)*a_ip(ji,jj,jl), v_il(ji,jj,jl) )
-                        IF ( zdvice > epsi10 ) then
+                        
+                        IF ( zdvice > epsi10 ) THEN
+                        
                            v_il (ji,jj,jl) = v_il (ji,jj,jl)   - zdvice
-                           v_ip(ji,jj,jl)  = v_ip(ji,jj,jl)    + zdvice
+                           v_ip(ji,jj,jl)  = v_ip(ji,jj,jl)    + zdvice ! MV: not sure i understand dh_i_sum seems counted twice - 
+                                                                        ! as it is already counted in surface melt
 !                          zvolp(ji,jj)     = zvolp(ji,jj)     + zdvice ! pointless to calculate total volume (done in icevar)
 !                          zfpond(ji,jj)    = fpond(ji,jj)     + zdvice ! pointless to follow fw budget (ponds have no fw)
@@ -577,6 +675,9 @@
                            ENDIF
                            h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) !!! could get a division by zero here
+                           
+                           diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) + zdvice   ! diag
+                           
                         ENDIF
-                  
+                        
                      !----------------------------------------------------------------
                      ! Freeze pre-existing lid 
@@ -598,4 +699,7 @@
 !                          zfpond(ji,jj)   = zfpond(ji,jj)    - zdvice
                            h_ip(ji,jj,jl)   = v_ip(ji,jj,jl) / a_ip(ji,jj,jl)
+                           
+                           diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice    ! diag
+                           
                         ENDIF
                   
@@ -624,4 +728,6 @@
                         v_il (ji,jj,jl)  = v_il(ji,jj,jl)   + zdvice
                         v_ip(ji,jj,jl)   = v_ip(ji,jj,jl)   - zdvice
+                        
+                        diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice      ! diag
 !                       zvolp(ji,jj)     = zvolp(ji,jj)     - zdvice
 !                       zfpond(ji,jj)    = zfpond(ji,jj)    - zdvice
@@ -849,5 +955,5 @@
              DO ns = 1, jl
                 cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl) &
-                   - rhos/rhow * &   ! free air fraction that can be filled by water
+                   - rhos/rhow * &     ! free air fraction that can be filled by water
                      asnon(ns)  * &    ! effective areal fraction of snow in that category
                      max(min(hsnon(ns)+alfan(ns)-alfan(jl), alfan(jl+1)-alfan(jl)), 0._wp)
@@ -870,4 +976,7 @@
           drain = zvolp(ji,jj) - cum_max_vol(jpl) + epsi10
           zvolp(ji,jj) = zvolp(ji,jj) - drain ! update meltwater volume available
+
+          diag_dvpn_rnf(ji,jj) = - drain      ! diag - overflow counted in the runoff part (arbitrary choice)
+          
           zdvolp(ji,jj) = drain         ! this is the drained water
           IF (zvolp(ji,jj) < epsi10) THEN
@@ -878,4 +987,5 @@
 
        ! height and area corresponding to the remaining volume
+       ! routine leaves zvolp unchanged
        CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index)
 
@@ -917,5 +1027,5 @@
 
                 perm = 0._wp ! MV ugly dummy patch
-                ! CALL ice_thd_pnd_perm(t_i(ji,jj,:,jl),  sz_i(ji,jj,:,jl), perm) ! bof 
+                CALL ice_thd_pnd_perm(t_i(ji,jj,:,jl),  sz_i(ji,jj,:,jl), perm) ! bof 
                 IF (perm > 0._wp) permflag = 1
 
@@ -924,4 +1034,7 @@
                 zdvolp(ji,jj) = zdvolp(ji,jj) + min(drain, zvolp(ji,jj))
                 zvolp(ji,jj) = max(zvolp(ji,jj) - drain, 0._wp)
+
+                diag_dvpn_drn(ji,jj) = - drain ! diag (could be better coded)
+                
                 IF (zvolp(ji,jj) < epsi10) THEN
                    zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj)
@@ -950,9 +1063,14 @@
        ! snow in melt ponds is not melted
        !------------------------------------------------------------------------
+       
+       ! MV here, it seems that we remove some meltwater from the ponds, but I can't really tell
+       ! how much, so I did not diagnose it
+       ! so if there is a problem here, nobody is going to see it...
+       
 
        ! Calculate pond volume for lower categories
        DO jl = 1,m_index-1
           v_ip(ji,jj,jl) = a_ip(ji,jj,jl) * h_ip(ji,jj,jl) & ! what is not in the snow
-                     - (rhos/rhow) * asnon(jl) * min(hsnon(jl), h_ip(ji,jj,jl))
+                         - (rhos/rhow) * asnon(jl) * min(hsnon(jl), h_ip(ji,jj,jl))
        END DO
 
