diff --git a/src/fesom_module.F90 b/src/fesom_module.F90
index 4f97013fb..9ea20678a 100755
--- a/src/fesom_module.F90
+++ b/src/fesom_module.F90
@@ -650,6 +650,8 @@ subroutine fesom_runloop(current_nsteps)
         ! LA - this causes the blowup !
         !        f%ice%data(size(f%ice%data))      = f%ice%data(2)
         !        f%ice%data(size(f%ice%data)-1)    = f%ice%data(1)
+                f%ice%data(size(f%ice%data))%values   = f%ice%data(2)%values   ! m_ice -> m_ice_ib
+                f%ice%data(size(f%ice%data)-1)%values = f%ice%data(1)%values   ! a_ice -> a_ice_ib
         !!!!!!!!!!!!!!!!!!
         
         
diff --git a/src/icb_allocate.F90 b/src/icb_allocate.F90
index bb84c5afd..6353c85c5 100644
--- a/src/icb_allocate.F90
+++ b/src/icb_allocate.F90
@@ -128,7 +128,7 @@ subroutine allocate_icb(partit, mesh)
   hfb_flux_ib = 0.0
   hfbv_flux_ib = 0.0
   lhfb_flux_ib = 0.0
-  allocate(arr_block(15*ib_num))
+  allocate(arr_block(16*ib_num))
   allocate(elem_block(ib_num))
   allocate(pe_block(ib_num))
 
@@ -138,7 +138,7 @@ subroutine allocate_icb(partit, mesh)
   call MPI_Bcast(elem_area_glob, elem2D, MPI_DOUBLE, 0, MPI_COMM_FESOM, MPIERR)
 
   allocate(vl_block(4*ib_num))
-  allocate(buoy_props(ib_num,13))
+  allocate(buoy_props(ib_num,14))
   buoy_props = 0.0
   allocate(melted(ib_num))
   melted = .false.
diff --git a/src/icb_coupling.F90 b/src/icb_coupling.F90
index 26a9b323d..6cac81eb1 100644
--- a/src/icb_coupling.F90
+++ b/src/icb_coupling.F90
@@ -60,122 +60,111 @@ subroutine prepare_icb2fesom(mesh, partit, ib,i_have_element,localelement,depth_
 #include "associate_part_ass.h"
 #include "associate_mesh_ass.h"
 
-    if(i_have_element) then                                     ! PE has iceberg
-        dz = 0.0                                                ! vertical distance over which heat flux is applied 
+    ! Also guard localelement: when an iceberg crosses a PE boundary during
+    ! iceberg_step1, i_have_element can remain .true. while local_idx_of()
+    ! returned 0 for the new element.  elem2d_nodes(i, 0) is an out-of-bounds
+    ! read that produces a garbage node index and may divide by zero area -> NaN.
+    if(i_have_element .and. localelement > 0) then
+        dz = 0.0
         allocate(tot_area_nods_in_ib_elem(mesh%nl))
 
-        num_ib_nods_in_ib_elem=0                                ! number of nodes in this element ???
-        tot_area_nods_in_ib_elem=0.0                            ! total area associated to nodes of iceberg element
-        idx_d = 0                                               ! index of level directly below or at iceberg base
+        num_ib_nods_in_ib_elem=0
+        tot_area_nods_in_ib_elem=0.0
+        idx_d = 0
 
-        ! loop over all three nodes of element
         do i=1,3
-            
-            ! assign node to iceberg_node
             iceberg_node=elem2d_nodes(i,localelement)
 
-            ! loop over all levels of iceberg_node
+            ! LOOP: consider all neighboring pairs (n_up,n_low) of 3D nodes
+            ! below n2..
+            !innerloop: do k=1, nl+1
             do k=1, nlevels_nod2D(iceberg_node)
-                
-                idx_d(i) = k                                    ! idx_d holds current level index until ...
-                                                                ! ... end of iceberg or bottom topography is reached
-                lev_up  = mesh%zbar_3d_n(k, iceberg_node)       ! upper level
+                idx_d(i) = k
+                lev_up  = mesh%zbar_3d_n(k, iceberg_node)
 
-                if( k==nlevels_nod2D(iceberg_node) ) then       ! if lower most level is reached ...
-                    lev_low = mesh%zbar_n_bot(iceberg_node)     ! ... lower level is equal to bottom topography
+                if( k==nlevels_nod2D(iceberg_node) ) then
+                    lev_low = mesh%zbar_n_bot(iceberg_node)
                 else
-                    lev_low = mesh%zbar_3d_n(k+1, iceberg_node) ! ... otherwise, lower level is one level below upper one
+                    lev_low = mesh%zbar_3d_n(k+1, iceberg_node)
                 end if
 
-                if( abs(lev_low)==abs(lev_up) ) then            ! if upper level is equal to lower level - when should this happen?
-                    idx_d(i) = idx_d(i) - 1                     ! level index is set back by one and exit loop
+                if( abs(lev_low)==abs(lev_up) ) then
+                    idx_d(i) = idx_d(i) - 1
                     exit
-                else if( abs(lev_low)>=abs(depth_ib) ) then     ! if lower level is below iceberg ...
-                                                                ! ... i.e. end of iceberg is reached, exit loop
+                else if( abs(lev_low)>=abs(depth_ib) ) then
                     exit
                 else
                     cycle
                 end if
             end do
-            ! at the end of the loop, idx_n holds the index of the level ...
-            ! ... directly below the iceberg or at iceberg base
-            ! -------------------------------------------------------------------
-
-            if (iceberg_node<=mydim_nod2d) then                                 ! if iceberg node on PE ...
-                ib_nods_in_ib_elem(i)           = iceberg_node                  ! ... add to list ib_nods_in_ib_elem
-                num_ib_nods_in_ib_elem          = num_ib_nods_in_ib_elem + 1    ! ... increase num_ib_nods_in_ib_elem by one
-                tot_area_nods_in_ib_elem        = tot_area_nods_in_ib_elem + mesh%area(:,iceberg_node)  ! increase tot_area_nods_in_ib_elem
+
+            if (iceberg_node<=mydim_nod2d) then
+                ib_nods_in_ib_elem(i)           = iceberg_node
+                num_ib_nods_in_ib_elem          = num_ib_nods_in_ib_elem + 1
+                tot_area_nods_in_ib_elem        = tot_area_nods_in_ib_elem + mesh%area(:,iceberg_node)
             else
-                ib_nods_in_ib_elem(i)           = 0                             ! ... otherwise, add zero to list ib_nods_in_ib_elem
+                ib_nods_in_ib_elem(i)           = 0
             end if
         end do
-        ! at the then of the loop, ib_nods_in_ib_elem holds the nodes of the iceberg elemen - why so complicated?
-        ! ... num_ib_nods_in_ib_elem holds the number of nodes assigned to PE
-        ! ... tot_area_nods_in_ib_elem holds the total area of the (up to three) nodes assigned to iceberg_elem
-        ! -------------------------------------------------------------------
 
-        ! loop over (up to three) nodes of iceberg_elem
         do i=1, 3
             iceberg_node=ib_nods_in_ib_elem(i)
+            ! Skip nodes with no ocean column; do NOT skip valid cavity nodes
+            if (iceberg_node <= 0 .or. ulevels_nod2d(iceberg_node) == 0) cycle
 
-            ! check if iceberg_nod is in cavity, cycle of .true.
-            if ((ulevels_nod2d(iceberg_node) == 0 ) .or. (use_cavity .and. ulevels_nod2d(iceberg_node) > 1)) cycle
-
-            ! if iceberg_node in PE, convert freshwater flux to flux density by dividing with tot_area_nods_in_ib_elem ...
-            ! ... of upper most level. The total iceberg flux is distributed among up to three nodes and only applied to the ...
-            ! ... surface layer.
             if (iceberg_node>0) then
-                ibfwbv(iceberg_node) = ibfwbv(iceberg_node) - fwbv_flux_ib(ib) / tot_area_nods_in_ib_elem(1)
-                ibfwb(iceberg_node) = ibfwb(iceberg_node) - fwb_flux_ib(ib) / tot_area_nods_in_ib_elem(1)
-                ibfwl(iceberg_node) = ibfwl(iceberg_node) - fwl_flux_ib(ib) / tot_area_nods_in_ib_elem(1)
-                ibfwe(iceberg_node) = ibfwe(iceberg_node) - fwe_flux_ib(ib) / tot_area_nods_in_ib_elem(1)
+                ! Guard: tot_area at level 1 is zero when all element nodes are cavity nodes
+                ! (mesh%area(1,cavity_node)==0), which would give division by zero / NaN.
+                if (tot_area_nods_in_ib_elem(1) > 0.0) then
+                    ibfwbv(iceberg_node) = ibfwbv(iceberg_node) - fwbv_flux_ib(ib) / tot_area_nods_in_ib_elem(1)
+                    ibfwb(iceberg_node) = ibfwb(iceberg_node) - fwb_flux_ib(ib) / tot_area_nods_in_ib_elem(1)
+                    ibfwl(iceberg_node) = ibfwl(iceberg_node) - fwl_flux_ib(ib) / tot_area_nods_in_ib_elem(1)
+                    ibfwe(iceberg_node) = ibfwe(iceberg_node) - fwe_flux_ib(ib) / tot_area_nods_in_ib_elem(1)
+                end if
+                !ibhf(iceberg_node) = ibhf(iceberg_node) - hfb_flux_ib(ib) / tot_area_nods_in_ib_elem(1)
+
+                ! Guard against idx_d=0 (can happen if nlevels_nod2D is 0 or loop didn't execute)
+                if (idx_d(i) <= 0) cycle
 
-                ! loop from surface to level below or at iceberg base
                 do j=1,idx_d(i)
-                    lev_up  = mesh%zbar_3d_n(j, iceberg_node)           ! upper level
-                    if( j==nlevels_nod2D(iceberg_node) ) then           ! if bottom level is reached ...
-                        lev_low = mesh%zbar_n_bot(iceberg_node)         ! ... lower level is set to bottom topography
+                    lev_up  = mesh%zbar_3d_n(j, iceberg_node)
+                    if( j==nlevels_nod2D(iceberg_node) ) then
+                        lev_low = mesh%zbar_n_bot(iceberg_node)
                     else
-                        lev_low = mesh%zbar_3d_n(j+1, iceberg_node)     ! ... otherwise, lower level is one level below upper one
+                        lev_low = mesh%zbar_3d_n(j+1, iceberg_node)
                     end if
-                    dz = abs( lev_low - lev_up )                        ! distance between lower and upper level
-                    
-                    ! check if lower level is below iceberg base and upper level is above iceberg base ...
+                    dz = abs( lev_low - lev_up )
                     if( (abs(lev_low)>=abs(depth_ib)) .and. (abs(lev_up)<abs(depth_ib)) ) then 
-                        dz = abs(abs(lev_up) - abs(depth_ib))           ! ... if so, distance dz is calculated as difference ...
-                                                                        ! ... between upper level and iceberg base
+                        dz = abs(abs(lev_up) - abs(depth_ib))
                     end if              
-                  
-                    ! check if iceberg is not yet melted away ...
-                    if( abs(depth_ib) > 0.0 ) then
-
-                        ! ... if so, heat flxues due to iceberg melting applied to layer interface (level) j ...
-                        ! ... is calculated nby weighting the lateral heat fluxes (hfl & hfbv) with the ratio dz/depth_ib ...
-                        ! ... and divided by the tot_area_nods_in_ib_elem of level j.
+                   
+                    ! Also guard tot_area(j): at levels where only some nodes extend,
+                    ! the remaining nodes contribute zero area, making tot_area very small.
+                    if( abs(depth_ib) > 0.0 .and. tot_area_nods_in_ib_elem(j) > 0.0 ) then
                         ibhf_n(j,iceberg_node) = ibhf_n(j,iceberg_node) & 
-                                                    - ((hfbv_flux_ib(ib,j)+hfl_flux_ib(ib,j)) * (dz / abs(depth_ib))) & 
+                                                    - (hfbv_flux_ib(ib,j)+hfl_flux_ib(ib,j)) & 
                                                     / tot_area_nods_in_ib_elem(j)
                     end if
                 end do
                 
-                ! if idx_d is not only upper most level ...
                 if( idx_d(i) > 1 ) then
-
-                    ! ... assign 50% of basal heat flux to level idx_d, i.e. the level below iceberg base ...
-                    ibhf_n(idx_d(i),iceberg_node) = ibhf_n(idx_d(i),iceberg_node) - 0.5 * hfb_flux_ib(ib) / tot_area_nods_in_ib_elem(idx_d(i))
-                    ! ... and 50% to idx_d-1, i.e. the level above iceberg base
-                    ibhf_n(idx_d(i)-1,iceberg_node) = ibhf_n(idx_d(i)-1,iceberg_node) - 0.5 * hfb_flux_ib(ib) / tot_area_nods_in_ib_elem(idx_d(i)-1)
-                else
-                    
-                    ! ... otherwise (should not happen), assign 100% of basal heat flux to level idx_n, i.e. the surface level
-                    ibhf_n(idx_d(i),iceberg_node) = ibhf_n(idx_d(i),iceberg_node) - hfb_flux_ib(ib) / tot_area_nods_in_ib_elem(idx_d(i))
+                    if (tot_area_nods_in_ib_elem(idx_d(i)) > 0.0) &
+                        ibhf_n(idx_d(i),iceberg_node) = ibhf_n(idx_d(i),iceberg_node) - 0.5 * hfb_flux_ib(ib) / tot_area_nods_in_ib_elem(idx_d(i))
+                    if (tot_area_nods_in_ib_elem(idx_d(i)-1) > 0.0) &
+                        ibhf_n(idx_d(i)-1,iceberg_node) = ibhf_n(idx_d(i)-1,iceberg_node) - 0.5 * hfb_flux_ib(ib) / tot_area_nods_in_ib_elem(idx_d(i)-1)
+                else if( idx_d(i) == 1 ) then
+                    if (tot_area_nods_in_ib_elem(idx_d(i)) > 0.0) &
+                        ibhf_n(idx_d(i),iceberg_node) = ibhf_n(idx_d(i),iceberg_node) - hfb_flux_ib(ib) / tot_area_nods_in_ib_elem(idx_d(i))
                 end if
-               
-                ! if iceberg not yet melted away ...
-                if( height_ib_single .ne. 0.0 ) then
-                    
-                    ! ... assign heat flux due to wave erosion to surface level
-                    ibhf_n(1,iceberg_node) = ibhf_n(1,iceberg_node) - hfe_flux_ib(ib) & 
+                
+                ! Wave erosion heat applied at level 1 only on open-ocean nodes (ulevels==1).
+                ! Cavity nodes (ulevels>1) skip nz=1 in oce_ale_tracer, so setting ibhf_n(1)
+                ! for them wastes the flux and overcounts the heat on the open-ocean node.
+                if( height_ib_single .ne. 0.0 .and. tot_area_nods_in_ib_elem(1) > 0.0 &
+                        .and. ulevels_nod2d(iceberg_node) == 1 ) then
+                    ibhf_n(1,iceberg_node) = ibhf_n(1,iceberg_node) - hfe_flux_ib(ib) &
+                            !* abs(height_ib_single) &
+                            !* ((abs(height_ib_single)-abs(depth_ib))/abs(height_ib_single)) &
                             / tot_area_nods_in_ib_elem(1)
                 end if
             end if
@@ -190,7 +179,7 @@ subroutine icb2fesom(mesh, partit, ice)
     use o_param
 
 ! kh 18.03.21 specification of structure used
-    use o_arrays, only: water_flux, heat_flux, wiso_flux_oce
+    use o_arrays, only: water_flux, heat_flux, wiso_flux_oce, Tclim_ib, is_nonlinfs
     use MOD_MESH
     use g_config
     use MOD_PARTIT
@@ -211,7 +200,6 @@ subroutine icb2fesom(mesh, partit, ice)
     do n=1, myDim_nod2d+eDim_nod2D
         if (ulevels_nod2d(n) > 1) cycle
         if (.not.turn_off_fw) then
-                ! LA add heat fluxes here...
                 water_flux(n)  = water_flux(n) - (ibfwb(n)+ibfwl(n)+ibfwe(n)+ibfwbv(n)) !* steps_per_ib_step
         end if
     end do
@@ -225,6 +213,30 @@ subroutine icb2fesom(mesh, partit, ice)
     end do
 !$OMP END PARALLEL DO
   end if
+
+!---enthalpy-correction-begin
+! The ALE kinematic boundary condition in oce_ale_tracer.F90:
+!   bc_surface = -dt*(heat_flux(n)/vcpw + sval*water_flux(n)*is_nonlinfs)
+! uses sval=T_ocean for the entire water_flux, including the iceberg FW
+! contribution. Iceberg meltwater enters at T_melt~=0 degC, not T_ocean.
+! When T_ocean < 0 degC the wrong sign of the kinematic FW term creates
+! spurious cooling: d(T*h)/dt includes T_ocean*dh instead of T_melt*dh=0.
+!
+! Fix: add vcpw*T_ocean*fw_iceberg to heat_flux. This cancels the erroneous
+! T_ocean*fw_ib term in bc_surface, leaving a net contribution of zero
+! (correct for T_melt=0 degC). The sign ensures warming for T_ocean<0 and
+! cooling for T_ocean>0, matching the physics of mixing with 0-degC water.
+  if (.not. (turn_off_hf .or. turn_off_fw)) then
+!$OMP PARALLEL DO
+    do n=1, myDim_nod2d+eDim_nod2D
+        if (use_cavity .and. ulevels_nod2d(n) > 1) cycle
+        heat_flux(n) = heat_flux(n) + vcpw * is_nonlinfs * Tclim_ib(1, n) &
+                       * (ibfwb(n) + ibfwl(n) + ibfwe(n) + ibfwbv(n))
+    end do
+!$OMP END PARALLEL DO
+  end if
+!---enthalpy-correction-end
+
 !---wiso-code-begin
     if(lwiso) then
       do n=1, myDim_nod2D+eDim_nod2D
diff --git a/src/icb_dyn.F90 b/src/icb_dyn.F90
index 7fa5b0508..334bd3f9d 100644
--- a/src/icb_dyn.F90
+++ b/src/icb_dyn.F90
@@ -76,6 +76,7 @@ subroutine iceberg_dyn(mesh, partit, ice, dynamics, ib, new_u_ib, new_v_ib, u_ib
  real			:: T_ave_ib, S_ave_ib, T_keel_ib, S_keel_ib
  character,  intent(IN)	:: file3*80
  real, intent(IN)	:: rho_icb
+ integer                :: max_nlev_alloc
  
 type(t_ice)   , intent(inout), target :: ice
 type(t_mesh), intent(in) , target :: mesh
@@ -87,19 +88,26 @@ subroutine iceberg_dyn(mesh, partit, ice, dynamics, ib, new_u_ib, new_v_ib, u_ib
 #include "associate_mesh_ass.h"
   
 n2=elem2D_nodes(1,iceberg_elem)
-allocate(arr_uo_dz(3,nlevels_nod2D(n2)))
-allocate(arr_vo_dz(3,nlevels_nod2D(n2)))
-allocate(arr_T_dz(3,nlevels_nod2D(n2)))
-allocate(arr_S_dz(3,nlevels_nod2D(n2)))
+! Use max nlevels across all 3 element nodes for allocation (cavity-safe)
+max_nlev_alloc = nlevels_nod2D(elem2D_nodes(1,iceberg_elem))
+do n=2,3
+  max_nlev_alloc = max(max_nlev_alloc, nlevels_nod2D(elem2D_nodes(n,iceberg_elem)))
+end do
+
+n2=elem2D_nodes(1,iceberg_elem)
+allocate(arr_uo_dz(3,max_nlev_alloc))
+allocate(arr_vo_dz(3,max_nlev_alloc))
+allocate(arr_T_dz(3,max_nlev_alloc))
+allocate(arr_S_dz(3,max_nlev_alloc))
 arr_uo_dz = 0.0
 arr_vo_dz = 0.0
 arr_T_dz = 0.0
 arr_S_dz = 0.0
 
-allocate(arr_uo_ib(nlevels_nod2D(n2)))
-allocate(arr_vo_ib(nlevels_nod2D(n2)))
-allocate(arr_T_ave_ib(nlevels_nod2D(n2)))
-allocate(arr_S_ave_ib(nlevels_nod2D(n2)))
+allocate(arr_uo_ib(max_nlev_alloc))
+allocate(arr_vo_ib(max_nlev_alloc))
+allocate(arr_T_ave_ib(max_nlev_alloc))
+allocate(arr_S_ave_ib(max_nlev_alloc))
 arr_uo_ib = 0.0
 arr_vo_ib = 0.0
 arr_T_ave_ib = 0.0
@@ -355,9 +363,15 @@ subroutine iceberg_acceleration(mesh, partit, dynamics, ib, au_ib, av_ib, Ao, Aa
  vel_atm = sqrt(ua_ib**2 + va_ib**2)
  wave_amplitude = 0.5 * 0.02025 * vel_atm**2
  
- !assume that waves have same direction as the winds
- direction_u = ua_ib / vel_atm
- direction_v = va_ib / vel_atm
+ !assume that waves have same direction as the winds;
+ !guard against zero wind speed to avoid division by zero
+ if(vel_atm > 0.) then
+   direction_u = ua_ib / vel_atm
+   direction_v = va_ib / vel_atm
+ else
+   direction_u = 0.
+   direction_v = 0.
+ end if
   
  !absolute values of relative velocities 
  abs_omib = sqrt( (uo_ib - u_ib)**2 + (vo_ib - v_ib)**2 )
@@ -610,6 +624,10 @@ end subroutine compute_areas
 
 
 subroutine iceberg_levelwise_andkeel(mesh, partit, dynamics, uo_dz,vo_dz, uo_keel,vo_keel, T_dz,S_dz, T_keel,S_keel, depth_ib,iceberg_elem, ib, ib_n_lvls)
+  !--------------------------------------------------------------------------
+  ! Per-level velocities and T/S at iceberg location, plus keel values.
+  ! Uses Z_3d_n_ib (mid-level depths) for consistent index mapping to UV_ib.
+  !--------------------------------------------------------------------------
   USE MOD_MESH
   use o_param
   use MOD_PARTIT
@@ -639,9 +657,9 @@ subroutine iceberg_levelwise_andkeel(mesh, partit, dynamics, uo_dz,vo_dz, uo_kee
   REAL, dimension(:,:,:), pointer :: UV_ib
 
   real           :: lev_up, lev_low
-  integer        :: m, k, n2, n_up, n_low, cavity_count, max_node_level_count
+  integer        :: m, k, n2, max_node_level_count, safe_lev
   ! depth over which is integrated (layer and sum)
-  real           :: dz, ufkeel1, ufkeel2, Temkeel, Salkeel, ldepth_up, ldepth_low, dz_depth
+  real           :: dz, ufkeel1, ufkeel2, Temkeel, Salkeel
 
 type(t_mesh), intent(in) , target :: mesh
 type(t_dyn), intent(in) , target :: dynamics
@@ -652,7 +670,6 @@ subroutine iceberg_levelwise_andkeel(mesh, partit, dynamics, uo_dz,vo_dz, uo_kee
 #include "associate_mesh_ass.h"
 
   UV_IB     => dynamics%uv_ib(:,:,:)
-  cavity_count=0
 
   do m=1,3
     if(m==1) then
@@ -685,89 +702,71 @@ subroutine iceberg_levelwise_andkeel(mesh, partit, dynamics, uo_dz,vo_dz, uo_kee
    !for each 2D node of the iceberg element..
    n2=elem2D_nodes(m,iceberg_elem)
 
-   ! LOOP: consider all neighboring pairs (n_up,n_low) of 3D nodes
-   ! below n2..
-   !innerloop: do k=1, nl+1
-   innerloop: do k=1, nlevels_nod2D(n2)
-    lev_up  = mesh%zbar_3d_n(k, n2)
-    !lev_up  = mesh%Z_3d_n(k, n2)
-    ldepth_up = mesh%Z_3d_n(k, n2)
-
-    if( k==nlevels_nod2D(n2) ) then
-        lev_low = mesh%zbar_n_bot(n2)
-        ldepth_low = mesh%zbar_n_bot(n2)
+   ! LOOP over mid-levels: Z_3d_n_ib(k) gives depth where UV_ib(k) lives
+   innerloop: do k=1, nlevels_nod2d(n2)
+    ! Skip levels inside the ice shelf for cavity nodes.  Tclim_ib/UV_ib at
+    ! k < ulevels_nod2d are never updated by the tracer solver; they hold stale
+    ! initial values (frozen at model start) that drive extreme spurious melt.
+    if (use_cavity .AND. mesh%cavity_depth(n2) /= 0.0 &
+        .AND. k < ulevels_nod2d(n2)) cycle
+
+    if( k==1 ) then
+        lev_up = 0.0
+    else if (use_cavity .AND. mesh%cavity_depth(n2) /= 0.0 &
+             .AND. k == ulevels_nod2d(n2)) then
+        ! First valid ocean level for this cavity node: the top boundary is the
+        ! ice shelf base, i.e. zbar_3d_n(ulevels_nod2d, n2) = zbar_n_srf(n2).
+        ! Using Z_3d_n_ib(k-1) here would give a dummy ice-shelf mid-depth.
+        lev_up = mesh%zbar_3d_n(k, n2)
+    else
+        lev_up = mesh%Z_3d_n_ib(k-1, n2)
+    end if
+    if( k <= nlevels_nod2D(n2)-1 ) then
+        lev_low = mesh%Z_3d_n_ib(k, n2)
     else
-        lev_low = mesh%zbar_3d_n(k+1, n2)
-        ldepth_low = mesh%Z_3d_n(k+1, n2)
-        !lev_low = mesh%Z_3d_n(k+1, n2)
+        lev_low = mesh%zbar_n_bot(n2)          ! bottom boundary
     end if
 
-    
-    !if( k==1 ) then
-    !    lev_up = 0.0
-    !else
-    !    lev_up  = mesh%Z_3d_n_ib(k-1, n2)
-    !    !lev_up  = mesh%Z_3d_n_ib(k-1, n2)
-    !end if
-
-    !if( k==nlevels_nod2D(n2) ) then
-    !    lev_low = mesh%zbar_n_bot(n2)
-    !else
-    !    lev_low = mesh%Z_3d_n_ib(k, n2)
-    !end if
+    if (lev_up==lev_low) then
+      arr_ib_n_lvls(m) = k
+      exit innerloop
+    end if
     dz = abs( lev_low - lev_up )
-    dz_depth = abs( ldepth_low - ldepth_up )
-    
-    !if( abs(lev_up)>=abs(depth_ib) ) then
-    !    ! ...icb bottom above lev_up --> no further integration
-    !end if
-    
-    !if( (abs(coord_nod3D(3, n_low))>abs(depth_ib)) .AND. (abs(coord_nod3D(3, n_up))>abs(depth_ib)) ) then
-    ! write(*,*) 'INFO, k:',k,'z_up:',coord_nod3D(3, n_up),'z_lo:',coord_nod3D(3, n_low),'depth:',depth_ib,'cavity:',(cavity_flag_nod2d(elem2D_nodes(m,iceberg_elem))==1)
-    !end if
-
-    ! if cavity node ..
-if (use_cavity .AND. mesh%cavity_depth(elem2D_nodes(m,iceberg_elem)) /= 0.0 .AND. abs(depth_ib) < abs(lev_up)) then
-    ! LA: Never go here for k=1, because abs(depth_ib)>=0.0 for all icebergs
-      
-      uo_dz(m,k)=UV_ib(1,k-1,n2)*abs(depth_ib)
-      vo_dz(m,k)=UV_ib(2,k-1,n2)*abs(depth_ib)
-      uo_keel(m)=UV_ib(1,k-1,n2)
-      vo_keel(m)=UV_ib(2,k-1,n2)
-
-      T_dz(m,k)=Tclim_ib(k-1,n2)*abs(depth_ib)
-      S_dz(m,k)=Sclim_ib(k-1,n2)*abs(depth_ib)
-      T_keel(m)=Tclim_ib(k-1,n2)
-      S_keel(m)=Sclim_ib(k-1,n2) ! check those choices with RT: OK
 
+    ! Cavity check: use zbar_3d_n for geometric comparison
+    if (use_cavity .AND. mesh%cavity_depth(elem2D_nodes(m,iceberg_elem)) /= 0.0 &
+        .AND. abs(depth_ib) < abs(mesh%zbar_3d_n(k, n2))) then
+      ! Iceberg draft is above the ocean surface at this cavity node.
+      ! Clamp to shallowest valid ocean level (ulevels_nod2d) to avoid
+      ! reading ice-shelf level T/S which are never initialised by the
+      ! tracer solver and may be garbage / NaN.
+      safe_lev = max(k-1, ulevels_nod2d(n2))
+      uo_keel(m)=UV_ib(1,safe_lev,n2)
+      vo_keel(m)=UV_ib(2,safe_lev,n2)
+      T_keel(m)=Tclim_ib(safe_lev,n2)
+      S_keel(m)=Sclim_ib(safe_lev,n2)
+      uo_dz(m,k)=UV_ib(1,safe_lev,n2)
+      vo_dz(m,k)=UV_ib(2,safe_lev,n2)
+      T_dz(m,k)=Tclim_ib(safe_lev,n2)
+      S_dz(m,k)=Sclim_ib(safe_lev,n2)
       exit innerloop
 
-    ! if the lowest z coord is below the iceberg draft, exit
-    !else if( abs(coord_nod3D(3, n_low))>=abs(depth_ib) .AND. abs(coord_nod3D(3, n_up))<=abs(depth_ib) ) then
-    
-    !****************************************************************
-    ! LA 23.11.21 case if depth_ib<lev_up
-else !# comp cav flag
-    if( abs(lev_low)>=abs(depth_ib) ) then !.AND. (abs(lev_up)<=abs(depth_ib)) ) then
-      if( abs(lev_up)<abs(depth_ib) ) then
-        dz = abs ( lev_up - depth_ib )
-      else
-        ! LA: Never go here, when starting with k=1
-        dz = abs(depth_ib)
-      end if
-    !****************************************************************
+    ! Keel layer: mid-level k is at or below the iceberg draft
+    else if( abs(lev_low)>=abs(depth_ib) ) then
       if( k==1 ) then
-
-        ufkeel1 = UV_ib(1,k,n2)
-        ufkeel2 = UV_ib(2,k,n2)
-        Temkeel = Tclim_ib(k,n2)
-        Salkeel = Sclim_ib(k,n2)
-      else if( k.eq.nlevels_nod2D(n2) ) then
+        ! Draft within first half-layer: use first mid-level value
+        ufkeel1 = UV_ib(1,1,n2)
+        ufkeel2 = UV_ib(2,1,n2)
+        Temkeel = Tclim_ib(1,n2)
+        Salkeel = Sclim_ib(1,n2)
+      else if( k > nlevels_nod2D(n2)-1 ) then
+        ! Last half-layer to bottom: no UV_ib(k), use piecewise constant
         ufkeel1 = UV_ib(1,k-1,n2)
         ufkeel2 = UV_ib(2,k-1,n2)
         Temkeel = Tclim_ib(k-1,n2)
         Salkeel = Sclim_ib(k-1,n2)
       else
+        ! Interpolate keel values between mid-levels k-1 and k
         ufkeel1 = interpol1D(abs(lev_up),UV_ib(1,k-1,n2),abs(lev_low),UV_ib(1,k,n2),abs(depth_ib))
         ufkeel2 = interpol1D(abs(lev_up),UV_ib(2,k-1,n2),abs(lev_low),UV_ib(2,k,n2),abs(depth_ib))
         Temkeel = interpol1D(abs(lev_up),Tclim_ib(k-1,n2),abs(lev_low),Tclim_ib(k,n2),abs(depth_ib))
@@ -786,53 +785,32 @@ subroutine iceberg_levelwise_andkeel(mesh, partit, dynamics, uo_dz,vo_dz, uo_kee
 
       arr_ib_n_lvls(m) = k
       exit innerloop
-    
-    !****************************************************************
-    ! LA 23.11.21 case if lev_low==0
-    else if(lev_low==lev_up) then
-      arr_ib_n_lvls(m) = k
-      exit innerloop
-    !****************************************************************
-    
-    else
-      if( k==1 ) then
 
-        uo_dz(m,k)=UV_ib(1,k,n2)
-        vo_dz(m,k)=UV_ib(2,k,n2)
-        T_dz(m,k)=Tclim_ib(k,n2)
-        S_dz(m,k)=Sclim_ib(k,n2)
-      elseif (k.eq.nlevels_nod2D(n2)) then  ! LA 2023-08-31
-
-        ! .. and sum up the layer-integrated velocities ..
-  ! kh 08.03.21 use UV_ib buffered values here      
+    ! Regular layer: iceberg extends deeper
+    else
+      if( k > nlevels_nod2D(n2)-1 ) then
+        ! Last half-layer to bottom: no UV_ib(k), use piecewise constant
         uo_dz(m,k)=UV_ib(1,k-1,n2)
         vo_dz(m,k)=UV_ib(2,k-1,n2)
         T_dz(m,k)=Tclim_ib(k-1,n2)
         S_dz(m,k)=Sclim_ib(k-1,n2)
-      else  
-
-        uo_dz(m,k)=0.5*(UV_ib(1,k-1,n2)+UV_ib(1,k,n2))
-        vo_dz(m,k)=0.5*(UV_ib(2,k-1,n2)+UV_ib(2,k,n2))
-        T_dz(m,k)=0.5*(Tclim_ib(k-1,n2)+ Tclim_ib(k,n2))
-        S_dz(m,k)=0.5*(Sclim_ib(k-1,n2)+ Sclim_ib(k,n2))
-      end if
-      
-      if (k.eq.nlevels_nod2D(n2)) then  ! LA 2023-08-31
-
         uo_keel(m)=UV_ib(1,k-1,n2)
         vo_keel(m)=UV_ib(2,k-1,n2)
-
         T_keel(m)=Tclim_ib(k-1,n2)
         S_keel(m)=Sclim_ib(k-1,n2)
-      else  
+      else
+        ! Store per-level value (UV_ib(k) lives at Z_3d_n_ib(k))
+        uo_dz(m,k)=UV_ib(1,k,n2)
+        vo_dz(m,k)=UV_ib(2,k,n2)
+        T_dz(m,k)=Tclim_ib(k,n2)
+        S_dz(m,k)=Sclim_ib(k,n2)
+        ! Update keel to deepest level so far
         uo_keel(m)=UV_ib(1,k,n2)
         vo_keel(m)=UV_ib(2,k,n2)
-
         T_keel(m)=Tclim_ib(k,n2)
         S_keel(m)=Sclim_ib(k,n2)
       end if
     end if
-end if !cavity
 
    end do innerloop
  end do nodeloop !loop over all nodes of iceberg element
@@ -854,6 +832,12 @@ real function interpol1D(x0,f0,x1,f1,x)
    real, intent(IN) :: x0,f0,x1,f1,x
    real :: frac
    
+   ! Guard against near-zero denominator (thin ALE layer).
+   ! An exact-zero check is insufficient; 0*Inf = NaN when x≈x0 too.
+   if (abs(x1-x0) < 1.0e-6) then
+     interpol1D = 0.5*(f0+f1)
+     return
+   end if
    frac = (f1 - f0)/(x1 - x0)
    interpol1D = f0 + frac * (x - x0)
   	
@@ -865,6 +849,10 @@ end subroutine iceberg_levelwise_andkeel
 
 
 subroutine iceberg_average_andkeel(mesh, partit, dynamics, uo_dz,vo_dz, uo_keel,vo_keel, T_dz,S_dz, T_keel,S_keel, depth_ib,iceberg_elem, ib)
+  !--------------------------------------------------------------------------
+  ! Depth-averaged velocities and T/S from surface to iceberg draft, plus
+  ! keel values.  Uses Z_3d_n_ib with trapezoidal integration (as iceberg_avvelo).
+  !--------------------------------------------------------------------------
   USE MOD_MESH
   use o_param
   use MOD_PARTIT
@@ -891,7 +879,7 @@ subroutine iceberg_average_andkeel(mesh, partit, dynamics, uo_dz,vo_dz, uo_keel,
   REAL, dimension(:,:,:), pointer :: UV_ib
 
   real           :: lev_up, lev_low
-  integer        :: m, k, n2, n_up, n_low, cavity_count
+  integer        :: m, k, n2, safe_lev
   ! depth over which is integrated (layer and sum)
   real           :: dz, ufkeel1, ufkeel2, Temkeel, Salkeel
 
@@ -904,7 +892,6 @@ subroutine iceberg_average_andkeel(mesh, partit, dynamics, uo_dz,vo_dz, uo_keel,
 #include "associate_mesh_ass.h"
 
   UV_IB     => dynamics%uv_ib(:,:,:)
-  cavity_count=0
 
   !LOOP: over all nodes of the iceberg element
   nodeloop: do m=1, 3
@@ -920,97 +907,90 @@ subroutine iceberg_average_andkeel(mesh, partit, dynamics, uo_dz,vo_dz, uo_keel,
    T_keel(m)=0.0
    S_keel(m)=0.0
 
-   ! LOOP: consider all neighboring pairs (n_up,n_low) of 3D nodes
-   ! below n2..
-   !innerloop: do k=1, nl+1
-   innerloop: do k=1, nlevels_nod2D(n2)
-    lev_up  = mesh%zbar_3d_n(k, n2)
+   ! LOOP over mid-levels: Z_3d_n_ib(k) gives depth where UV_ib(k) lives
+   innerloop: do k=1, nlevels_nod2d(n2)
+    ! Skip levels inside the ice shelf for cavity nodes.  Tclim_ib/UV_ib at
+    ! k < ulevels_nod2d are never updated by the tracer solver; they hold stale
+    ! initial values (frozen at model start) that drive extreme spurious melt.
+    if (use_cavity .AND. mesh%cavity_depth(n2) /= 0.0 &
+        .AND. k < ulevels_nod2d(n2)) cycle
 
-    if( k==nlevels_nod2D(n2) ) then
-        lev_low = mesh%zbar_n_bot(n2)
+    if( k==1 ) then
+        lev_up = 0.0
+    else if (use_cavity .AND. mesh%cavity_depth(n2) /= 0.0 &
+             .AND. k == ulevels_nod2d(n2)) then
+        ! First valid ocean level for this cavity node: the top boundary is the
+        ! ice shelf base, i.e. zbar_3d_n(ulevels_nod2d, n2) = zbar_n_srf(n2).
+        ! Using Z_3d_n_ib(k-1) here would give a dummy ice-shelf mid-depth.
+        lev_up = mesh%zbar_3d_n(k, n2)
     else
-        lev_low = mesh%zbar_3d_n(k+1, n2)
+        lev_up = mesh%Z_3d_n_ib(k-1, n2)
     end if
 
-    !if( k==1 ) then
-    !    lev_up = 0.0
-    !else
-    !    lev_up  = mesh%Z_3d_n_ib(k-1, n2)
-    !    !lev_up  = mesh%Z_3d_n_ib(k-1, n2)
-    !end if
-
-    !if( k==nlevels_nod2D(n2) ) then
-    !    lev_low = mesh%zbar_n_bot(n2)
-    !else
-    !    lev_low = mesh%Z_3d_n_ib(k, n2)
-    !end if
-    dz = abs( lev_low - lev_up )
-    
-    !if( abs(lev_up)>=abs(depth_ib) ) then
-    !    ! ...icb bottom above lev_up --> no further integration
-    !end if
-    
-    !if( (abs(coord_nod3D(3, n_low))>abs(depth_ib)) .AND. (abs(coord_nod3D(3, n_up))>abs(depth_ib)) ) then
-    ! write(*,*) 'INFO, k:',k,'z_up:',coord_nod3D(3, n_up),'z_lo:',coord_nod3D(3, n_low),'depth:',depth_ib,'cavity:',(mesh%cavity_flag_n(elem2D_nodes(m,iceberg_elem))==1)
-    !end if
-
-if (use_cavity .AND. mesh%cavity_depth(elem2D_nodes(m,iceberg_elem)) /= 0.0 .AND. abs(depth_ib) < abs(lev_up)) then
-    ! LA: Never go here for k=1, because abs(depth_ib)>=0.0 for all icebergs
-      
-      uo_dz(m)=UV_ib(1,k-1,n2)*abs(depth_ib)
-      vo_dz(m)=UV_ib(2,k-1,n2)*abs(depth_ib)
-      uo_keel(m)=UV_ib(1,k-1,n2)
-      vo_keel(m)=UV_ib(2,k-1,n2)
+    if( k <= nlevels_nod2D(n2)-1 ) then
+        lev_low = mesh%Z_3d_n_ib(k, n2)
+    else
+        lev_low = mesh%zbar_n_bot(n2)          ! bottom boundary
+    end if
 
-      T_dz(m)=Tclim_ib(k-1,n2)*abs(depth_ib)
-      S_dz(m)=Sclim_ib(k-1,n2)*abs(depth_ib)
-      T_keel(m)=Tclim_ib(k-1,n2)
-      S_keel(m)=Sclim_ib(k-1,n2) ! check those choices with RT: OK
+    if (lev_up==lev_low) then
+      exit innerloop
+    end if
+    dz = abs( lev_low - lev_up )
 
+    ! Cavity check: use zbar_3d_n for geometric comparison
+    if (use_cavity .AND. mesh%cavity_depth(elem2D_nodes(m,iceberg_elem)) /= 0.0 &
+        .AND. abs(depth_ib) < abs(mesh%zbar_3d_n(k, n2))) then
+      ! Iceberg draft is above the ocean surface at this cavity node.
+      ! Clamp to shallowest valid ocean level (ulevels_nod2d) to avoid
+      ! reading ice-shelf level T/S which are never initialised by the
+      ! tracer solver and may be garbage / NaN.
+      safe_lev = max(k-1, ulevels_nod2d(n2))
+      uo_dz(m)=UV_ib(1,safe_lev,n2)*abs(depth_ib)
+      vo_dz(m)=UV_ib(2,safe_lev,n2)*abs(depth_ib)
+      uo_keel(m)=UV_ib(1,safe_lev,n2)
+      vo_keel(m)=UV_ib(2,safe_lev,n2)
+      T_dz(m)=Tclim_ib(safe_lev,n2)*abs(depth_ib)
+      S_dz(m)=Sclim_ib(safe_lev,n2)*abs(depth_ib)
+      T_keel(m)=Tclim_ib(safe_lev,n2)
+      S_keel(m)=Sclim_ib(safe_lev,n2)
       exit innerloop
 
-    ! if the lowest z coord is below the iceberg draft, exit
-    !else if( abs(coord_nod3D(3, n_low))>=abs(depth_ib) .AND. abs(coord_nod3D(3, n_up))<=abs(depth_ib) ) then
-    
-    !****************************************************************
-    ! LA 23.11.21 case if depth_ib<lev_up
-else !# comp cav flag
-    if( abs(lev_low)>=abs(depth_ib) ) then !.AND. (abs(lev_up)<=abs(depth_ib)) ) then
-      if( abs(lev_up)<abs(depth_ib) ) then
-        dz = abs ( lev_up - depth_ib )
-      else
-        ! LA: Never go here, when starting with k=1
-        dz = abs(depth_ib)
-      end if
-    !****************************************************************
+    ! Keel layer: mid-level k is at or below the iceberg draft
+    else if( abs(lev_low)>=abs(depth_ib) ) then
+      dz = abs( lev_up - depth_ib )
 
       if( k==1 ) then
-        ufkeel1 = UV_ib(1,k,n2)
-        ufkeel2 = UV_ib(2,k,n2)
-        Temkeel = Tclim_ib(k,n2)
-        Salkeel = Sclim_ib(k,n2)
+        ! Draft within first half-layer: piecewise constant
+        ufkeel1 = UV_ib(1,1,n2)
+        ufkeel2 = UV_ib(2,1,n2)
+        Temkeel = Tclim_ib(1,n2)
+        Salkeel = Sclim_ib(1,n2)
         uo_dz(m)=ufkeel1*dz 
         vo_dz(m)=ufkeel2*dz
         T_dz(m)=Temkeel*dz
         S_dz(m)=Salkeel*dz
-      else if( k.eq.nlevels_nod2D(n2) ) then
+      else if( k > nlevels_nod2D(n2)-1 ) then
+        ! Last half-layer to bottom: no UV_ib(k), use piecewise constant
         ufkeel1 = UV_ib(1,k-1,n2)
         ufkeel2 = UV_ib(2,k-1,n2)
         Temkeel = Tclim_ib(k-1,n2)
         Salkeel = Sclim_ib(k-1,n2)
-        uo_dz(m)=uo_dz(m)+ 0.5*(UV_ib(1,k-1,n2) + ufkeel1)*dz 
-        vo_dz(m)=vo_dz(m)+ 0.5*(UV_ib(2,k-1,n2) + ufkeel2)*dz
-        T_dz(m)=T_dz(m)+ 0.5*(Tclim_ib(k-1,n2)+ Temkeel)*dz
-        S_dz(m)=S_dz(m)+ 0.5*(Sclim_ib(k-1,n2)+ Salkeel)*dz
+        uo_dz(m)=uo_dz(m)+ ufkeel1*dz
+        vo_dz(m)=vo_dz(m)+ ufkeel2*dz
+        T_dz(m)=T_dz(m)+ Temkeel*dz
+        S_dz(m)=S_dz(m)+ Salkeel*dz
       else
+        ! Interpolate keel values between mid-levels k-1 and k
         ufkeel1 = interpol1D(abs(lev_up),UV_ib(1,k-1,n2),abs(lev_low),UV_ib(1,k,n2),abs(depth_ib))
         ufkeel2 = interpol1D(abs(lev_up),UV_ib(2,k-1,n2),abs(lev_low),UV_ib(2,k,n2),abs(depth_ib))
         Temkeel = interpol1D(abs(lev_up),Tclim_ib(k-1,n2),abs(lev_low),Tclim_ib(k,n2),abs(depth_ib))
         Salkeel = interpol1D(abs(lev_up),Sclim_ib(k-1,n2),abs(lev_low),Sclim_ib(k,n2),abs(depth_ib))
-        uo_dz(m)=uo_dz(m)+ 0.5*(UV_ib(1,k-1,n2) + ufkeel1)*dz 
+        ! Trapezoidal integration for partial layer to keel
+        uo_dz(m)=uo_dz(m)+ 0.5*(UV_ib(1,k-1,n2) + ufkeel1)*dz
         vo_dz(m)=vo_dz(m)+ 0.5*(UV_ib(2,k-1,n2) + ufkeel2)*dz
-        T_dz(m)=T_dz(m)+ 0.5*(Tclim_ib(k-1,n2)+ Temkeel)*dz
-        S_dz(m)=S_dz(m)+ 0.5*(Sclim_ib(k-1,n2)+ Salkeel)*dz
+        T_dz(m)=T_dz(m)+ 0.5*(Tclim_ib(k-1,n2) + Temkeel)*dz
+        S_dz(m)=S_dz(m)+ 0.5*(Sclim_ib(k-1,n2) + Salkeel)*dz
       end if
 
       uo_keel(m)=ufkeel1
@@ -1018,74 +998,52 @@ subroutine iceberg_average_andkeel(mesh, partit, dynamics, uo_dz,vo_dz, uo_keel,
 
       T_keel(m)=Temkeel
       S_keel(m)=Salkeel
-      
-      if(S_dz(m)/abs(depth_ib)>70.) then
-       write(*,*) 'innerloop, dz:',dz,', depth:',depth_ib,',S_dz(m):',S_dz(m),"m:",m,", k:",k,", Tclim_ib(k-1,n2):",Tclim_ib(k-1,n2),", Tclim_ib(k,n2):", Tclim_ib(k,n2),", Salkeel:",Salkeel,", lev_low:",lev_low,", lev_up:",lev_up
-      end if
-      
-      if(T_dz(m)/abs(depth_ib)>70.) then
-       write(*,*) 'innerloop, dz:',dz,', depth:',depth_ib,',T_dz(m):',T_dz(m),"m:",m,", k:",k,", Sclim_ib(k-1,n2):",Sclim_ib(k-1,n2),", Sclim_ib(k,n2):", Sclim_ib(k,n2),",Temkeel:",Temkeel,", lev_low:",lev_low,", lev_up:",lev_up
-      end if
-
       exit innerloop
-    
-    !****************************************************************
-    ! LA 23.11.21 case if lev_low==0
-    else if(lev_low==lev_up) then
-      exit innerloop
-    !****************************************************************
-    
-    else	
+
+    ! Regular layer: iceberg extends deeper
+    else
       if( k==1 ) then
-          cycle
+        ! First half-layer (surface to first mid-level): piecewise constant
+        uo_dz(m)=uo_dz(m)+ UV_ib(1,1,n2)*dz
+        vo_dz(m)=vo_dz(m)+ UV_ib(2,1,n2)*dz
+        T_dz(m)=T_dz(m)+ Tclim_ib(1,n2)*dz
+        S_dz(m)=S_dz(m)+ Sclim_ib(1,n2)*dz
+        cycle
       end if
 
-
-      if (k.eq.nlevels_nod2D(n2)) then  ! LA 2023-08-31
-        ! .. and sum up the layer-integrated velocities ..
-  ! kh 08.03.21 use UV_ib buffered values here      
+      if( k > nlevels_nod2D(n2)-1 ) then
+        ! Last half-layer to bottom: no UV_ib(k), use piecewise constant
         uo_dz(m)=uo_dz(m)+ UV_ib(1,k-1,n2)*dz
         vo_dz(m)=vo_dz(m)+ UV_ib(2,k-1,n2)*dz
         T_dz(m)=T_dz(m)+ Tclim_ib(k-1,n2)*dz
         S_dz(m)=S_dz(m)+ Sclim_ib(k-1,n2)*dz
-      else  
-        uo_dz(m)=uo_dz(m)+ 0.5*(UV_ib(1,k-1,n2)+UV_ib(1,k,n2))*dz
-        vo_dz(m)=vo_dz(m)+ 0.5*(UV_ib(2,k-1,n2)+UV_ib(2,k,n2))*dz
-        T_dz(m)=T_dz(m)+ 0.5*(Tclim_ib(k-1,n2)+ Tclim_ib(k,n2))*dz
-        S_dz(m)=S_dz(m)+ 0.5*(Sclim_ib(k-1,n2)+ Sclim_ib(k,n2))*dz
-      end if
-
-   
-      if(S_dz(m)/abs(depth_ib)>70.) then
-       write(*,*) 'innerloop, dz:',dz,', depth:',depth_ib,',S_dz(m):',S_dz(m),"m:",m,", k:",k,", Tclim_ib(k-1,n2):",Tclim_ib(k-1,n2),", Tclim_ib(k,n2):", Tclim_ib(k,n2),", lev_low:",lev_low,", lev_up:",lev_up
-      end if
-      
-      if(T_dz(m)/abs(depth_ib)>70.) then
-       write(*,*) 'innerloop, dz:',dz,', depth:',depth_ib,',T_dz(m):',T_dz(m),"m:",m,", k:",k,", Sclim_ib(k-1,n2):",Sclim_ib(k-1,n2),", Sclim_ib(k,n2):", Sclim_ib(k,n2),", lev_low:",lev_low,", lev_up:",lev_up
-      end if
-
-      if (k.eq.nlevels_nod2D(n2)) then  ! LA 2023-08-31
         uo_keel(m)=UV_ib(1,k-1,n2)
         vo_keel(m)=UV_ib(2,k-1,n2)
-
         T_keel(m)=Tclim_ib(k-1,n2)
         S_keel(m)=Sclim_ib(k-1,n2)
-      else  
+      else
+        ! Trapezoidal integration between consecutive mid-levels
+        uo_dz(m)=uo_dz(m)+ 0.5*(UV_ib(1,k-1,n2)+UV_ib(1,k,n2))*dz
+        vo_dz(m)=vo_dz(m)+ 0.5*(UV_ib(2,k-1,n2)+UV_ib(2,k,n2))*dz
+        T_dz(m)=T_dz(m)+ 0.5*(Tclim_ib(k-1,n2)+Tclim_ib(k,n2))*dz
+        S_dz(m)=S_dz(m)+ 0.5*(Sclim_ib(k-1,n2)+Sclim_ib(k,n2))*dz
+        ! Update keel to deepest level so far
         uo_keel(m)=UV_ib(1,k,n2)
         vo_keel(m)=UV_ib(2,k,n2)
-
         T_keel(m)=Tclim_ib(k,n2)
         S_keel(m)=Sclim_ib(k,n2)
       end if
     end if
-end if !cavity
+ 
    end do innerloop
 
    ! divide by depth over which was integrated
-   uo_dz(m)=uo_dz(m)/abs(depth_ib)
-   vo_dz(m)=vo_dz(m)/abs(depth_ib)
-   T_dz(m)=T_dz(m)/abs(depth_ib)
-   S_dz(m)=S_dz(m)/abs(depth_ib)
+   if (abs(depth_ib) > 0.0) then
+     uo_dz(m)=uo_dz(m)/abs(depth_ib)
+     vo_dz(m)=vo_dz(m)/abs(depth_ib)
+     T_dz(m)=T_dz(m)/abs(depth_ib)
+     S_dz(m)=S_dz(m)/abs(depth_ib)
+   end if
 
  end do nodeloop !loop over all nodes of iceberg element
        
@@ -1096,6 +1054,12 @@ real function interpol1D(x0,f0,x1,f1,x)
    real, intent(IN) :: x0,f0,x1,f1,x
    real :: frac
    
+   ! Guard against near-zero denominator (thin ALE layer).
+   ! An exact-zero check is insufficient; 0*Inf = NaN when x≈x0 too.
+   if (abs(x1-x0) < 1.0e-6) then
+     interpol1D = 0.5*(f0+f1)
+     return
+   end if
    frac = (f1 - f0)/(x1 - x0)
    interpol1D = f0 + frac * (x - x0)
   	
@@ -1193,10 +1157,11 @@ subroutine iceberg_avvelo(mesh, partit, dynamics, uo_dz,vo_dz,depth_ib,iceberg_e
 	 
     else	
       if( k==1 ) then
+        uo_dz(m)=uo_dz(m)+ UV_ib(1,k,n2)*dz
+        vo_dz(m)=vo_dz(m)+ UV_ib(2,k,n2)*dz
         cycle
       end if
       ! .. and sum up the layer-integrated velocities:
-! kh 08.03.21 use UV_ib buffered values here
       uo_dz(m)=uo_dz(m)+0.5*(UV_ib(1,k-1,n2)+UV_ib(1,k,n2))*dz
       vo_dz(m)=vo_dz(m)+0.5*(UV_ib(2,k-1,n2)+UV_ib(2,k,n2))*dz
 	 
diff --git a/src/icb_elem.F90 b/src/icb_elem.F90
index c87f6eb9c..c0a0c8421 100644
--- a/src/icb_elem.F90
+++ b/src/icb_elem.F90
@@ -111,25 +111,23 @@ subroutine nodal_average(mesh, partit, dynamics, local_idx, gradientx, gradienty
   return !not my node
  end if 
 
- do node=1,myDim_nod2D
-   do k=1, nod_in_elem2D_num(node) 
-     elem  = nod_in_elem2D(k,node) 
-
-     !do idx_elem = 1,  nod_in_elem2D(local_idx)%nmb
-     !elem = nod_in_elem2D(local_idx)%addresses(idx_elem)
-     !area_ = voltriangle(elem)
-     area_ = elem_area(elem)
-     patch= patch + area_
-     
-     !gradientx = gradientx + area * sum( ssh(elem2D_nodes(:,elem)) * bafux_2D(:,elem) )
-     !gradienty = gradienty + area * sum( ssh(elem2D_nodes(:,elem)) * bafuy_2D(:,elem) )
+ do k=1, nod_in_elem2D_num(local_idx)
+   elem  = nod_in_elem2D(k,local_idx)
 
-!LA 2023-03-07
+   !do idx_elem = 1,  nod_in_elem2D(local_idx)%nmb
+   !elem = nod_in_elem2D(local_idx)%addresses(idx_elem)
+   !area_ = voltriangle(elem)
+   area_ = elem_area(elem)
+   patch= patch + area_
+
+   !gradientx = gradientx + area * sum( ssh(elem2D_nodes(:,elem)) * bafux_2D(:,elem) )
+   !gradienty = gradienty + area * sum( ssh(elem2D_nodes(:,elem)) * bafuy_2D(:,elem) )
+
+! 2023-03-07
 eta_n_ib => dynamics%eta_n_ib(:)
-! kh 18.03.21 use eta_n_ib buffered values here
-     gradientx = gradientx + area_ * sum( eta_n_ib(elem2D_nodes(:,elem)) * gradient_sca(1:3, elem)) 
-     gradienty = gradienty + area_ * sum( eta_n_ib(elem2D_nodes(:,elem)) * gradient_sca(4:6, elem)) 
-   end do
+!h 18.03.21 use eta_n_ib buffered values here
+   gradientx = gradientx + area_ * sum( eta_n_ib(elem2D_nodes(:,elem)) * gradient_sca(1:3, elem))
+   gradienty = gradienty + area_ * sum( eta_n_ib(elem2D_nodes(:,elem)) * gradient_sca(4:6, elem))
  end do
  
  gradientx = gradientx / patch
@@ -399,13 +397,12 @@ subroutine iceberg_elem4all(mesh, partit, elem, lon_deg, lat_deg)
   if(i_have_element) then
    i_have_element= elem2D_nodes(1,elem) <= myDim_nod2D !1 PE still .true.
    if (use_cavity) then 
-      !reject_tmp = all( (mesh%cavity_depth(elem2D_nodes(:,elem))/=0.0) .OR. (mesh%bc_index_nod2D(elem2D_nodes(:,elem))==0.0) )
-      reject_tmp = any(mesh%cavity_depth(elem2D_nodes(:,elem))/=0.0) .OR. all(mesh%bc_index_nod2D(elem2D_nodes(:,elem))==0.0)
+      reject_tmp = all( (mesh%cavity_depth(elem2D_nodes(:,elem))/=0.0) .OR. (mesh%bc_index_nod2D(elem2D_nodes(:,elem))==0.0) )
       if(reject_tmp) then
       !if( reject_elem(mesh, partit, elem) ) then
        elem=0 !reject element
        i_have_element=.false.
-       write(*,*) 'elem4all: iceberg found in shelf region: elem = 0'
+       !write(*,*) 'elem4all: iceberg found in shelf region: elem = 0'
       else 
        elem=myList_elem2D(elem) !global now
       end if 
@@ -462,12 +459,17 @@ subroutine find_new_iceberg_elem(mesh, partit, old_iceberg_elem, pt, left_mype)
    old_iceberg_elem=elem_containing_n2
    if (use_cavity) then 
       !if( reject_elem(mesh, partit, old_iceberg_elem) ) then
-      !reject_tmp = all( (mesh%cavity_depth(elem2D_nodes(:,ibelem_tmp))/=0.0) .OR. (mesh%bc_index_nod2D(elem2D_nodes(:,ibelem_tmp))==0.0) )
-      reject_tmp = any(mesh%cavity_depth(elem2D_nodes(:,ibelem_tmp))/=0.0) .OR. all(mesh%bc_index_nod2D(elem2D_nodes(:,ibelem_tmp))==0.0)
+      reject_tmp = all( (mesh%cavity_depth(elem2D_nodes(:,elem_containing_n2))/=0.0) .OR. (mesh%bc_index_nod2D(elem2D_nodes(:,elem_containing_n2))==0.0) )
       if(reject_tmp) then
          left_mype=1.0
-         write(*,*) 'iceberg found in shelf region: left_mype = 1'
+         !write(*,*) 'iceberg found in shelf region: left_mype = 1'
          old_iceberg_elem=ibelem_tmp
+         ! Do NOT return here.  At a cavity boundary the trajectory endpoint
+         ! can lie inside both the cavity element and an adjacent non-cavity
+         ! neighbour (shared edge).  If the cavity element is encountered
+         ! first in the search order, returning immediately would miss the
+         ! valid non-cavity alternative and trigger a spurious left_mype=1.
+         cycle
       end if
    endif
    RETURN 
diff --git a/src/icb_modules.F90 b/src/icb_modules.F90
index ee309ee6a..6b3787792 100644
--- a/src/icb_modules.F90
+++ b/src/icb_modules.F90
@@ -141,8 +141,7 @@ logical function reject_elem(mesh, partit, elem)
 if (use_cavity) then
 ! kh 09.08.21 change index_nod2d -> bc_index_nod2d?
  if (.not. use_cavityonelem) then
-   reject_elem = any(mesh%cavity_depth(mesh%elem2D_nodes(:,elem))/=0.0) .OR. all(mesh%bc_index_nod2D(mesh%elem2D_nodes(:,elem))==0.0) 
-   !reject_elem = all( (mesh%cavity_depth(mesh%elem2D_nodes(:,elem))/=0.0) .OR. (mesh%bc_index_nod2D(mesh%elem2D_nodes(:,elem))==0.0) )
+   reject_elem = all( (mesh%cavity_depth(mesh%elem2D_nodes(:,elem))/=0.0) .OR. (mesh%bc_index_nod2D(mesh%elem2D_nodes(:,elem))==0.0) )
  !else
  end if
 else
diff --git a/src/icb_step.F90 b/src/icb_step.F90
index 62521a85d..4881dc090 100644
--- a/src/icb_step.F90
+++ b/src/icb_step.F90
@@ -56,14 +56,15 @@ subroutine iceberg_calculation(ice, mesh, partit, dynamics, istep)
  logical	:: firstcall=.true. 					!=
  logical	:: lastsubstep  					!=
 
- real		:: arr_from_block(15)					!=
+ real		:: arr_from_block(16)					!=
  integer	:: elem_from_block					!=  
  real		:: vl_from_block(4)					!=	
- real,dimension(15*ib_num):: arr_block_red				!=
+ real,dimension(16*ib_num):: arr_block_red				!=
  integer,dimension(ib_num):: elem_block_red				!=
  integer,dimension(ib_num):: pe_block_red				!=
  integer    :: n
  real, dimension(4*ib_num):: vl_block_red				!=
+ !integer,dimension(ib_num):: grounded_int, grounded_int_red		!=
 
 type(t_ice), intent(inout), target :: ice
 type(t_mesh), intent(in) , target :: mesh
@@ -143,10 +144,11 @@ subroutine iceberg_calculation(ice, mesh, partit, dynamics, istep)
  arr_block_red = 0.0
  elem_block_red= 0
  pe_block_red= 0
+
  vl_block_red = 0.0
 
 !$omp critical 
- call MPI_IAllREDUCE(arr_block, arr_block_red, 15*ib_num, MPI_DOUBLE_PRECISION, MPI_SUM, partit%MPI_COMM_FESOM_IB, req, partit%MPIERR_IB)
+ call MPI_IAllREDUCE(arr_block, arr_block_red, 16*ib_num, MPI_DOUBLE_PRECISION, MPI_SUM, partit%MPI_COMM_FESOM_IB, req, partit%MPIERR_IB)
 !$omp end critical
 
  completed = .false.
@@ -213,7 +215,7 @@ subroutine iceberg_calculation(ice, mesh, partit, dynamics, istep)
   !get the smaller array arr(15) and 
   !the iceberg element for iceberg ib
   !as before
-  arr_from_block = arr_block_red( (ib-1)*15+1 : ib*15)
+  arr_from_block = arr_block_red( (ib-1)*16+1 : ib*16)
   elem_from_block= elem_block_red(ib)
   !averaged volume losses
   vl_from_block = vl_block_red( (ib-1)*4+1 : ib*4)
@@ -294,7 +296,7 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib_single,lengt
  use g_rotate_grid	!for subroutine g2r, logfile_outfreq					!=
  use g_config, only: steps_per_ib_step
  !=
-use iceberg_params, only: length_ib, width_ib, scaling, elem_block, elem_area_glob !smallestvol_icb, arr_block, elem_block, l_geo_out, icb_outfreq, l_allowgrounding, draft_scale, reject_elem, melted, grounded, scaling !, length_ib, width_ib, scaling
+use iceberg_params, only: length_ib, width_ib, scaling, elem_block, elem_area_glob !, smallestvol_icb, arr_block, elem_block, l_geo_out, icb_outfreq, l_allowgrounding, draft_scale, reject_elem, melted, grounded !, length_ib, width_ib, scaling
 !#else
 ! use iceberg_params, only: smallestvol_icb, arr_block, elem_block, l_geo_out, icb_outfreq, l_allowgrounding, draft_scale, melted, grounded, scaling !, length_ib, width_ib, scaling
 !#endif
@@ -320,8 +322,9 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib_single,lengt
  
  integer, dimension(:), save, allocatable :: local_idx_of
  real      			:: depth_ib, volume_ib, mass_ib
+ real      			:: depth_ib_premelt, height_ib_premelt
  real				:: lon_rad, lat_rad, new_u_ib, new_v_ib
- real	   			:: old_lon,old_lat, frozen_in, P_ib, conci_ib
+ real	   			:: old_lon,old_lat, frozen_in, P_ib, conci_ib, grounded_ib
  real				:: lon_rad_out, lat_rad_out  !for unrotated output
  real				:: lon_deg_out, lat_deg_out  !for unrotated output
  integer   			:: i, iceberg_node  
@@ -339,7 +342,7 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib_single,lengt
  logical			:: l_output
  							
  !MPI											
- real	   			:: arr(15)						!=
+ real	   			:: arr(16)						!=
  logical   			:: i_have_element					!=
  real	   			:: left_mype						!=
  integer   			:: old_element						!=
@@ -412,12 +415,11 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib_single,lengt
    
    
    if (use_cavity) then
-      reject_tmp = any(mesh%cavity_depth(mesh%elem2D_nodes(:,iceberg_elem))/=0.0) .OR. all(mesh%bc_index_nod2D(mesh%elem2D_nodes(:,iceberg_elem))==0.0) 
-      !reject_tmp = all( (mesh%cavity_depth(mesh%elem2D_nodes(:,iceberg_elem))/=0.0) .OR. (mesh%bc_index_nod2D(mesh%elem2D_nodes(:,iceberg_elem))==0.0) )
+      reject_tmp = all( (mesh%cavity_depth(mesh%elem2D_nodes(:,iceberg_elem))/=0.0) .OR. (mesh%bc_index_nod2D(mesh%elem2D_nodes(:,iceberg_elem))==0.0) )
       if(reject_tmp) then
-       write(*,*) " * set IB elem ",iceberg_elem,"to zero for IB=",ib
-       write(*,*) " cavity: ",all((mesh%cavity_depth(mesh%elem2D_nodes(:,iceberg_elem))/=0.0))
-       write(*,*) " boundary: ", all(mesh%bc_index_nod2D(mesh%elem2D_nodes(:,iceberg_elem))==0)
+!       write(*,*) " * set IB elem ",iceberg_elem,"to zero for IB=",ib
+!       write(*,*) " cavity: ",all((mesh%cavity_depth(mesh%elem2D_nodes(:,iceberg_elem))/=0.0))
+!       write(*,*) " boundary: ", all(mesh%bc_index_nod2D(mesh%elem2D_nodes(:,iceberg_elem))==1)
        iceberg_elem=0 !reject element
        i_have_element=.false.
       else 
@@ -430,9 +432,9 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib_single,lengt
   call com_integer(partit, i_have_element,iceberg_elem)
  
   if(iceberg_elem .EQ. 0) then
-    write(*,*) 'IB ',ib,' rot. coords:', lon_deg, lat_deg !,lon_rad, lat_rad
-    call par_ex (partit%MPI_COMM_FESOM, partit%mype)
-    stop 'ICEBERG OUTSIDE MODEL DOMAIN OR IN ICE SHELF REGION'
+        write(*,*) 'IB ',ib,' rot. coords:', lon_deg, lat_deg !,lon_rad, lat_rad
+   	call par_ex (partit%MPI_COMM_FESOM, partit%mype)
+   	stop 'ICEBERG OUTSIDE MODEL DOMAIN OR IN ICE SHELF REGION'
   end if
   
   ! initialize the iceberg velocity
@@ -463,6 +465,7 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib_single,lengt
  
  arr=0.
  frozen_in = 0.
+ grounded_ib = 0.
  i_have_element=.false.
  !if the first node belongs to this processor.. (just one processor enters here!)
  !if( local_idx_of(iceberg_elem) > 0 .and. elem2D_nodes(1,local_idx_of(iceberg_elem)) <= myDim_nod2D ) then
@@ -477,6 +480,9 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib_single,lengt
   
   !===========================DYNAMICS===============================
   
+  ! Save pre-melt geometry for heat flux distribution (bug 4)
+  depth_ib_premelt  = depth_ib
+  height_ib_premelt = height_ib_single
 
   call iceberg_dyn(mesh, partit, ice, dynamics, ib, new_u_ib, new_v_ib, u_ib, v_ib, lon_rad,lat_rad, depth_ib, &
                    height_ib_single, length_ib_single, width_ib_single, local_idx_of(iceberg_elem), &
@@ -510,13 +516,14 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib_single,lengt
     old_lat = lat_rad
  
 ! kh 16.03.21 (asynchronous) iceberg calculation starts with the content in common arrays at istep and will merge its results at istep_end_synced
-    grounded(ib) = .true.
+    grounded_ib = 1.
     !if (mod(istep_end_synced,logfile_outfreq)==0) then
-    if (lverbose_icb) write(*,*) 'iceberg ib ', ib, 'is grounded'
+    if (lverbose_icb) write(*,*) 'iceberg ib ', ib, 'is grounded; ib_depth=',draft_scale(ib)*abs(depth_ib),'; Zdepth=',Zdepth
     !end if
  	
  else 
   !===================...ELSE CALCULATE TRAJECTORY====================
+    grounded_ib = 0.
 
  
  t0=MPI_Wtime()
@@ -535,8 +542,7 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib_single,lengt
    lon_rad = old_lon
    lat_rad = old_lat
  t4=MPI_Wtime()
-   ! LA: Does this work corretly?
-   call parallel2coast(mesh,partit, new_u_ib, new_v_ib, lon_rad,lat_rad, local_idx_of(iceberg_elem), ib, iceberg_elem)
+   call parallel2coast(mesh,partit, new_u_ib, new_v_ib, lon_rad,lat_rad, local_idx_of(iceberg_elem))
  t5=MPI_Wtime()
    call trajectory( lon_rad,lat_rad, new_u_ib,new_v_ib, new_u_ib,new_v_ib, &
 		   lon_deg,lat_deg,old_lon,old_lat, dt*REAL(steps_per_ib_step))
@@ -567,22 +573,23 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib_single,lengt
      case(2) 
       area_ib_tot = length_ib_single*width_ib_single*scaling(ib)
     end select
-!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(idx) REDUCTION(+:area_ib_tot)
-!$OMP DO
+!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(idx) REDUCTION(+:area_ib_tot)
     do idx = 1, size(elem_block)
         if (elem_block(idx) == iceberg_elem) then
             area_ib_tot = area_ib_tot + length_ib(idx) * width_ib(idx) * scaling(idx)
         end if
     end do
-!$OMP END DO
-!$OMP END PARALLEL
+!$OMP END PARALLEL DO
   !-----------------------------
 
     if ((area_ib_tot > elem_area_glob(iceberg_elem)) .and. (old_element.ne.0)) then ! .and. (left_mype == 0)) then 
         if( lverbose_icb) then
             write(*,*) " *******************************************************"
             write(*,*) " * set iceberg ", ib, " back to elem ", old_element, " from elem ", iceberg_elem
-            write(*,*) " * area_ib_tot = ", area_ib_tot, "; elem_area = ", elem_area(local_idx_of(iceberg_elem))
+            ! elem_area is indexed by local element number.  When the new
+            ! element belongs to a different PE, local_idx_of() returns 0
+            ! and elem_area(0) is an out-of-bounds read.  Use the global array.
+            write(*,*) " * area_ib_tot = ", area_ib_tot, "; elem_area = ", elem_area_glob(iceberg_elem)
         end if
         i_have_element = .true.
         left_mype = 0.0
@@ -599,15 +606,15 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib_single,lengt
  
   !values for communication
   arr= (/ height_ib_single,length_ib_single,width_ib_single, u_ib,v_ib, lon_rad,lat_rad, &
-          left_mype, old_lon,old_lat, frozen_in, dudt, dvdt, P_ib, conci_ib/) 
+          left_mype, old_lon,old_lat, frozen_in, dudt, dvdt, P_ib, conci_ib, grounded_ib /) 
 
   !save in larger array	  
-  arr_block((ib-1)*15+1 : ib*15)=arr
+  arr_block((ib-1)*16+1 : ib*16)=arr
   elem_block(ib)=iceberg_elem
   pe_block(ib)=mype
 
   volume_ib=height_ib_single*length_ib_single*width_ib_single
-  call prepare_icb2fesom(mesh,partit,ib,i_have_element,local_idx_of(iceberg_elem),depth_ib,height_ib_single)
+  call prepare_icb2fesom(mesh,partit,ib,i_have_element,local_idx_of(iceberg_elem),depth_ib_premelt,height_ib_premelt)
  end if !processor has element?
 end if !... and first node belongs to processor?
 
@@ -642,14 +649,14 @@ subroutine iceberg_step2(mesh, partit,arr, elem_from_block, ib, height_ib_single
  use g_comm_auto
 !=
 use g_comm
-use iceberg_params, only: length_ib, width_ib, scaling !smallestvol_icb, arr_block, elem_block, l_geo_out, icb_outfreq, l_allowgrounding, draft_scale, reject_elem, melted, grounded, scaling !, length_ib, width_ib, scaling
+use iceberg_params, only: length_ib, width_ib, scaling !, smallestvol_icb, arr_block, elem_block, l_geo_out, icb_outfreq, l_allowgrounding, draft_scale, reject_elem, melted, grounded !, length_ib, width_ib, scaling
 !#else
 ! use iceberg_params, only: smallestvol_icb, buoy_props, bvl_mean, lvlv_mean, lvle_mean, lvlb_mean, ascii_out, l_geo_out, icb_outfreq, l_allowgrounding, draft_scale, elem_block
 !#endif
  												!=
  implicit none											!=
  
- real, 	  intent(in)	:: arr(15)
+ real, 	  intent(in)	:: arr(16)
  integer, intent(in)	:: elem_from_block
  integer, intent(in)	:: ib
  real,    intent(inout)	:: height_ib_single, length_ib_single, width_ib_single
@@ -749,6 +756,7 @@ subroutine iceberg_step2(mesh, partit,arr, elem_from_block, ib, height_ib_single
  dvdt	  = arr(13)
  P_ib	  = arr(14)
  conci_ib = arr(15) 
+ grounded(ib) = (arr(16) > 0.5)
 
  !**** check if iceberg melted in step 1 ****! 
  volume_ib = height_ib_single * length_ib_single * width_ib_single ! * rho_icb
@@ -790,16 +798,14 @@ subroutine iceberg_step2(mesh, partit,arr, elem_from_block, ib, height_ib_single
       case(2) 
        area_ib_tot = length_ib_single*width_ib_single*scaling(ib)
      end select
-!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(idx) REDUCTION(+:area_ib_tot)
-!$OMP DO
+!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(idx) REDUCTION(+:area_ib_tot)
      do idx = 1, size(elem_block_red)
          if (elem_block_red(idx) == iceberg_elem) then
      !        write(*,*) " * Found element ",elem_block_red(idx), " for ib ",idx, ", elem_area=",elem_area_tmp
              area_ib_tot = area_ib_tot + length_ib(idx) * width_ib(idx) * scaling(idx)
          end if
      end do
-!$OMP END DO
-!$OMP END PARALLEL
+!$OMP END PARALLEL DO
      if((area_ib_tot > elem_area_tmp) .and. (elem_area_tmp > 0.0) .and. (old_element.ne.0)) then
          if(mype==pe_block_red(ib) .and. lverbose_icb) then
             write(*,*) " *******************************************************"
@@ -816,9 +822,10 @@ subroutine iceberg_step2(mesh, partit,arr, elem_from_block, ib, height_ib_single
          u_ib    = 0.
          v_ib    = 0.
 
-         i_have_element= (local_idx_of(iceberg_elem) .ne. 0) 
+         iceberg_elem = old_element
+         i_have_element= (local_idx_of(old_element) .ne. 0)
          if(i_have_element) then
-            i_have_element= mesh%elem2D_nodes(1,local_idx_of(iceberg_elem)) <= partit%myDim_nod2D !1 PE still .true.
+            i_have_element= mesh%elem2D_nodes(1,local_idx_of(old_element)) <= partit%myDim_nod2D !1 PE still .true.
          end if
      end if
    else 
@@ -888,6 +895,7 @@ subroutine iceberg_step2(mesh, partit,arr, elem_from_block, ib, height_ib_single
   buoy_props(ib,11) = length_ib_single
   buoy_props(ib,12) = width_ib_single
   buoy_props(ib,13) = iceberg_elem
+  buoy_props(ib,14) = arr(16)
 
 ! end if
 
@@ -901,7 +909,7 @@ end subroutine iceberg_step2
 !****************************************************************************************************************************
 
 subroutine initialize_velo(mesh,partit,dynamics, i_have_element, ib, u_ib, v_ib, lon_rad, lat_rad, depth_ib, localelem)			
- 
+
  use g_rotate_grid,  only: vector_g2r
 ! use iceberg_params, only: l_initial, l_iniuser, ini_u, ini_v
 implicit none
@@ -1036,16 +1044,14 @@ end subroutine depth_bathy
 !****************************************************************************************************************************
 !****************************************************************************************************************************
 
-subroutine parallel2coast(mesh, partit,u, v, lon,lat, elem, ib, elem_global)
+subroutine parallel2coast(mesh, partit,u, v, lon,lat, elem)
  use iceberg_params, only: coastal_nodes
- use g_config, only: lverbose_icb
- 
  implicit none
  
  real, intent(inout) 	:: u, v 	!velocity
  real, intent(in)	:: lon, lat 	!radiant
  integer, intent(in)	:: elem
- integer, intent(in)    :: ib, elem_global
+ 
  integer :: fld_tmp
  integer, dimension(3) :: n
  integer :: node, m, i
@@ -1058,85 +1064,47 @@ subroutine parallel2coast(mesh, partit,u, v, lon,lat, elem, ib, elem_global)
 #include "associate_mesh_def.h"
 #include "associate_part_ass.h"
 #include "associate_mesh_ass.h"
- 
+  
   if( use_cavity ) then
     !fld_tmp = coastal_nodes(mesh, partit, elem)
     fld_tmp =  count( (mesh%cavity_depth(mesh%elem2D_nodes(:,elem))/=0.0) .OR. (mesh%bc_index_nod2D(mesh%elem2D_nodes(:,elem))==0.0) )
   else
     fld_tmp =  count( (mesh%bc_index_nod2D(mesh%elem2D_nodes(:,elem))==0.0) )
   end if
- 
-  if( lverbose_icb ) then
-    write(*,*) " - running parallel2coast for  IB ",ib," on elem ",elem_global,", fld_tmp=",fld_tmp
-    write(*,*) " - u_old, v_old =",u,v
-  end if
-
+  
   SELECT CASE ( fld_tmp ) !num of coastal points
    CASE (0) !...coastal points: do nothing
     return
     
    CASE (1) !...coastal point
    n = 0
-   i = 1
+   i = 2
    velocity = [ u, v ]
+    ! Classify element nodes: n(1) = barrier node, n(2),n(3) = non-barrier nodes
     do m = 1, 3
       node = mesh%elem2D_nodes(m,elem)
-      !write(*,*) 'index ', m, ':', index_nod2D(node)
       if( use_cavity ) then
-        !if( mesh%bc_index_nod2D(node)==1 .OR. cavity_flag_nod2d(node)==1 ) then
         if( mesh%bc_index_nod2D(node)==0.0 .OR.  (mesh%cavity_depth(node)/=0.0) ) then
-         n(i) = node
-         exit
-        end if 
-      else
-        if( mesh%bc_index_nod2D(node)==0.0 ) then
-         n(i) = node
-         exit
-        end if 
-      end if
-    end do 
-    
-   !write(*,*) 'one coastal node ', n(1)
-  
-  !LA comment for testing
-   i = 2 
-   !if ( n(1) <= myDim_nod2D ) then	!all neighbours known
-    do m = 1, 3 !nghbr_nod2D(n(1))%nmb
-      node = mesh%elem2D_nodes(m,elem)
-      if( use_cavity ) then
-        if ( (node /= n(1)) .and. ( (mesh%bc_index_nod2D(node)==0.0) .OR. (mesh%cavity_depth(node)/=0.0) ) ) then   
-            n(i) = node
-            i = i+1
-            if(i==4) exit
+          n(1) = node
+        else
+          n(i) = node
+          i = i + 1
         end if
       else
-        if ( (node /= n(1)) .and. ( (mesh%bc_index_nod2D(node)==0.0))) then
-            n(i) = node
-            i = i+1
-            if(i==4) exit
+        if( mesh%bc_index_nod2D(node)==1 ) then
+          n(1) = node
+        else
+          n(i) = node
+          i = i + 1
         end if
       end if
     end do
-    
-    !write(*,*) 'nodes n(i) ', n
-    
-    d1 = sqrt( (lon - coord_nod2D(1, n(2)))**2 + (lat - coord_nod2D(2, n(2)))**2 )
-    d2 = sqrt( (lon - coord_nod2D(1, n(3)))**2 + (lat - coord_nod2D(2, n(3)))**2 )
-    !write(*,*) 'distances :' , d1, d2
-    !write(*,*) 'velocity vor :' , velocity
-    if (d1 < d2) then
-      call projection(mesh,partit,  velocity, n(2), n(1))
-    else
-      call projection(mesh,partit,  velocity, n(3), n(1))
+
+    ! Project velocity along edge between the two non-barrier nodes
+    ! (the open-ocean face opposite the barrier node)
+    if (n(2) /= 0 .AND. n(3) /= 0) then
+      call projection(mesh,partit,  velocity, n(2), n(3))
     end if
-    !write(*,*) 'velocity nach:', velocity
-    !call projection(velocity, n(3), n(2))
-      
-   !else
-   ! !if coastal point is not first node of element, the coastal point could be in eDim_nod2D,
-   ! !so not all neighbours of this node are known to PE. WHAT SHOULD BE DONE?
-   !end if    
-    
     
    CASE (2) !...coastal points
     n = 0
@@ -1151,7 +1119,7 @@ subroutine parallel2coast(mesh, partit,u, v, lon,lat, elem, ib, elem_global)
          i = i+1
         end if
       else
-        if( mesh%bc_index_nod2D(node)==0.0 ) then
+        if( mesh%bc_index_nod2D(node)==1 ) then
          n(i) = node
          i = i+1
         end if
@@ -1160,16 +1128,18 @@ subroutine parallel2coast(mesh, partit,u, v, lon,lat, elem, ib, elem_global)
     call projection(mesh,partit,  velocity, n(1), n(2))
     
    
-   CASE DEFAULT 
-    return  	!mesh element MUST NOT have 3 coastal points!
+   CASE DEFAULT
+    ! All 3 nodes flagged cavity/boundary.  This element should have been
+    ! rejected by find_new_iceberg_elem / iceberg_elem4all before reaching
+    ! here.  Velocity is not modified; step2's iceberg_elem4all will stop
+    ! the iceberg by reverting it to the previous element with zero velocity.
+    if (mype==0) write(*,*) 'WARNING: parallel2coast called with fully-cavity element', elem
+    return
 
  END SELECT
  
  u = velocity(1)
  v = velocity(2)
-  if( lverbose_icb ) then
-    write(*,*) " - u_new, v_new =",u,v
-  end if
 
 end subroutine parallel2coast
 
@@ -1558,15 +1528,15 @@ subroutine determine_save_count(partit)
   ! computes save_count_buoys and prev_sec_in_year from records in existing netcdf file
   !-----------------------------------------------------------  
   use g_clock
-  use netcdf
 !  use iceberg_params, only : file_icb_netcdf, save_count_buoys, prev_sec_in_year
   !use iceberg_params, only : save_count_buoys, prev_sec_in_year
   implicit none
+
+#include "netcdf.inc" 
   integer		    :: buoy_nrec
   integer                   :: status, ncid
   integer                   :: dimid_rec
   integer                   :: time_varid
-  integer                   :: start_1d(1), count_1d(1)
 type(t_partit), intent(inout), target :: partit
 #include "associate_part_def.h"
 #include "associate_part_ass.h"
@@ -1574,14 +1544,14 @@ subroutine determine_save_count(partit)
  if (mype==0) then
 
   ! open file
-  status = nf90_open(file_icb_netcdf, nf90_nowrite, ncid)
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_open(file_icb_netcdf, nf_nowrite, ncid)
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
   ! inquire time dimension ID and its length
-  status = nf90_inq_dimid(ncid, 'time', dimid_rec)
-  if(status .ne. nf90_noerr) call handle_err(status, partit)
-  status = nf90_inquire_dimension(ncid, dimid_rec, len=buoy_nrec)
-  if(status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_inq_dimid(ncid, 'time', dimid_rec)
+  if(status .ne. nf_noerr) call handle_err(status, partit)
+  status = nf_inq_dimlen(ncid, dimid_rec, buoy_nrec)
+  if(status .ne. nf_noerr) call handle_err(status, partit)
 
   ! the next buoy/iceberg record to be saved
   save_count_buoys=buoy_nrec+1
@@ -1589,16 +1559,16 @@ subroutine determine_save_count(partit)
   write(*,*) 'next record is #',save_count_buoys
 
   ! load sec_in_year up to now in 'prev_sec_in_year', time axis will be continued
-  status=nf90_inq_varid(ncid, 'time', time_varid)
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
-  status=nf90_get_var(ncid, time_varid, prev_sec_in_year) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status=nf_inq_varid(ncid, 'time', time_varid)
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  status=nf_get_vara_double(ncid, time_varid, save_count_buoys-1, 1, prev_sec_in_year) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
   write(*,*) 'seconds passed up to now: ',prev_sec_in_year
 
   !close file
-  status=nf90_close(ncid)
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status=nf_close(ncid)
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
 end if
 
@@ -1615,18 +1585,19 @@ subroutine init_buoy_output(partit)
   !-----------------------------------------------------------  
   use g_clock
   use g_config, only : ib_num
-  use netcdf
 !  use iceberg_params, only : file_icb_netcdf, save_count_buoys !ggf in namelist
   implicit none
-  integer                   :: status,ncid,year_start,month_start,day_start
+
+#include "netcdf.inc" 
+  integer                   :: status,ncid
   integer                   :: dimid_ib, dimid_rec, dimids(2)
   integer                   :: time_varid, iter_varid
   integer                   :: lonrad_id, latrad_id, londeg_id, latdeg_id
   integer                   :: frozen_id, dudt_id, dvdt_id
   integer                   :: uib_id, vib_id
   integer                   :: height_id, length_id, width_id
-  integer                   :: bvl_id, lvlv_id, lvle_id, lvlb_id, felem_id
-  character(100)            :: longname, att_text, description, units
+  integer                   :: bvl_id, lvlv_id, lvle_id, lvlb_id, felem_id, grounded_id
+  character(100)            :: longname, att_text, description
 type(t_partit), intent(inout), target :: partit
 #include "associate_part_def.h"
 #include "associate_part_ass.h"
@@ -1636,21 +1607,21 @@ subroutine init_buoy_output(partit)
  
 
  ! create a file
-  status = nf90_create(file_icb_netcdf, nf90_clobber, ncid)
-  if (status.ne.nf90_noerr)  call handle_err(status, partit)
+  status = nf_create(file_icb_netcdf, nf_clobber, ncid)
+  if (status.ne.nf_noerr)  call handle_err(status, partit)
 
   ! Define the dimensions
-  status = nf90_def_dim(ncid, 'number_tracer', ib_num, dimid_ib)
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
-  status = nf90_def_dim(ncid, 'time', nf90_unlimited, dimid_rec)
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_def_dim(ncid, 'number_tracer', ib_num, dimid_ib)
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  status = nf_def_dim(ncid, 'time', NF_UNLIMITED, dimid_rec)
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
 
   ! Define the time and iteration variables
-  status = nf90_def_var(ncid, 'time', nf90_double, dimids=(/dimid_rec/), varid=time_varid)
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
-  status = nf90_def_var(ncid, 'iter', nf90_int, dimids=(/dimid_rec/), varid=iter_varid)
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_def_var(ncid, 'time', NF_DOUBLE, 1, dimid_rec, time_varid)
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  status = nf_def_var(ncid, 'iter', NF_INT, 1, dimid_rec, iter_varid)
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
   ! Define the netCDF variables for the tracers.
   ! In Fortran, the unlimited dimension must come
@@ -1659,276 +1630,290 @@ subroutine init_buoy_output(partit)
   dimids(2) = dimid_rec
 
 
-  status = nf90_def_var(ncid, 'pos_lon_rad', nf90_double, dimids=dimids, varid=lonrad_id)
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_def_var(ncid, 'pos_lon_rad', NF_DOUBLE, 2, dimids, lonrad_id)
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
-  status = nf90_def_var(ncid, 'pos_lat_rad', nf90_double, dimids=dimids, varid=latrad_id)
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_def_var(ncid, 'pos_lat_rad', NF_DOUBLE, 2, dimids, latrad_id)
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
-  status = nf90_def_var(ncid, 'pos_lon_deg', nf90_double, dimids=dimids, varid=londeg_id)
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_def_var(ncid, 'pos_lon_deg', NF_DOUBLE, 2, dimids, londeg_id)
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
-  status = nf90_def_var(ncid, 'pos_lat_deg', nf90_double, dimids=dimids, varid=latdeg_id)
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_def_var(ncid, 'pos_lat_deg', NF_DOUBLE, 2, dimids, latdeg_id)
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
-  status = nf90_def_var(ncid, 'frozen_in', nf90_double, dimids=dimids, varid=frozen_id)
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_def_var(ncid, 'frozen_in', NF_DOUBLE, 2, dimids, frozen_id)
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
-  status = nf90_def_var(ncid, 'du_dt', nf90_double, dimids=dimids, varid=dudt_id)
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_def_var(ncid, 'du_dt', NF_DOUBLE, 2, dimids, dudt_id)
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
-  status = nf90_def_var(ncid, 'dv_dt', nf90_double, dimids=dimids, varid=dvdt_id)
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_def_var(ncid, 'dv_dt', NF_DOUBLE, 2, dimids, dvdt_id)
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
-  status = nf90_def_var(ncid, 'icb_vel_u', nf90_double, dimids=dimids, varid=uib_id)
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_def_var(ncid, 'icb_vel_u', NF_DOUBLE, 2, dimids, uib_id)
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
-  status = nf90_def_var(ncid, 'icb_vel_v', nf90_double, dimids=dimids, varid=vib_id)
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_def_var(ncid, 'icb_vel_v', NF_DOUBLE, 2, dimids, vib_id)
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
   ! 3 dimensions of the iceberg, comment for buoy case
 
-  status = nf90_def_var(ncid, 'height', nf90_double, dimids=dimids, varid=height_id)
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_def_var(ncid, 'height', NF_DOUBLE, 2, dimids, height_id)
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
-  status = nf90_def_var(ncid, 'length', nf90_double, dimids=dimids, varid=length_id)
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_def_var(ncid, 'length', NF_DOUBLE, 2, dimids, length_id)
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
-  status = nf90_def_var(ncid, 'width', nf90_double, dimids=dimids, varid=width_id)
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_def_var(ncid, 'width', NF_DOUBLE, 2, dimids, width_id)
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
   ! 4 additional iceberg variables (meltrates), comment for buoy case
 
-  status = nf90_def_var(ncid, 'bvl', nf90_double, dimids=dimids, varid=bvl_id)
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_def_var(ncid, 'bvl', NF_DOUBLE, 2, dimids, bvl_id)
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
-  status = nf90_def_var(ncid, 'lvlv', nf90_double, dimids=dimids, varid=lvlv_id)
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_def_var(ncid, 'lvlv', NF_DOUBLE, 2, dimids, lvlv_id)
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
-  status = nf90_def_var(ncid, 'lvle', nf90_double, dimids=dimids, varid=lvle_id)
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_def_var(ncid, 'lvle', NF_DOUBLE, 2, dimids, lvle_id)
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
-  status = nf90_def_var(ncid, 'lvlb', nf90_double, dimids=dimids, varid=lvlb_id)
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_def_var(ncid, 'lvlb', NF_DOUBLE, 2, dimids, lvlb_id)
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
   ! LA: add felem
-  status = nf90_def_var(ncid, 'felem', nf90_double, dimids=dimids, varid=felem_id)
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_def_var(ncid, 'felem', NF_DOUBLE, 2, dimids, felem_id)
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
+  status = nf_def_var(ncid, 'grounded', NF_DOUBLE, 2, dimids, grounded_id)
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  
   ! Assign long_name and units attributes to variables.
   longname='time' ! use NetCDF Climate and Forecast (CF) Metadata Convention
-  status=nf90_put_att(ncid, time_varid, 'long_name', trim(longname))
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
-  units='seconds since model start'
-  status=nf90_put_att(ncid, time_varid, 'units', trim(units))
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_PUT_ATT_TEXT(ncid, time_varid, 'long_name', len_trim(longname), trim(longname)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  write(att_text, '(a14,I4.4,a1,I2.2,a1,I2.2,a9)') 'seconds since ', yearstart, '-', 1, '-', 1, ' 00:00:00'
+  status = nf_PUT_ATT_TEXT(ncid, time_varid, 'units', len_trim(att_text), trim(att_text))
+  if (status .ne. nf_noerr) call handle_err(status, partit)
   if (include_fleapyear) then
       att_text='standard'
   else
       att_text='noleap'
   end if
-  status = nf90_put_att(ncid, time_varid, 'calendar', trim(att_text))
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, time_varid, 'calendar', len_trim(att_text), trim(att_text))
+  if (status .ne. nf_noerr) call handle_err(status, partit)
   longname='iteration_count'
-  status = nf90_put_att(ncid, iter_varid, 'long_name', trim(longname)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_PUT_ATT_TEXT(ncid, iter_varid, 'long_name', len_trim(longname), trim(longname)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
   longname='longitude of buoy/iceberg position in radiant'
-  status = nf90_put_att(ncid, lonrad_id, 'long_name', trim(longname)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
-  status = nf90_put_att(ncid, lonrad_id, 'units', 'radiant')
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_PUT_ATT_TEXT(ncid, lonrad_id, 'long_name', len_trim(longname), trim(longname)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, lonrad_id, 'units', 7, 'radiant')
+  if (status .ne. nf_noerr) call handle_err(status, partit)
   !rotated or not rotated due to setting in iceberg module
   description='(un)rotated according to setting of l_geo_out in iceberg module'
-  status = nf90_put_att(ncid, lonrad_id, 'description', trim(description)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, lonrad_id, 'description', len_trim(description), trim(description)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
   longname='latitude of buoy/iceberg position in radiant'
-  status = nf90_put_att(ncid, latrad_id, 'long_name', trim(longname)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
-  status = nf90_put_att(ncid, latrad_id, 'units', 'radiant')
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_PUT_ATT_TEXT(ncid, latrad_id, 'long_name', len_trim(longname), trim(longname)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, latrad_id, 'units', 7, 'radiant')
+  if (status .ne. nf_noerr) call handle_err(status, partit)
   !rotated or not rotated due to setting in iceberg module
   description='(un)rotated according to setting of l_geo_out in iceberg module'
-  status = nf90_put_att(ncid, latrad_id, 'description', trim(description)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, latrad_id, 'description', len_trim(description), trim(description)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
   longname='longitude of buoy/iceberg position in degree'
-  status = nf90_put_att(ncid, londeg_id, 'long_name', trim(longname)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
-  status = nf90_put_att(ncid, londeg_id, 'units', 'degree')
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_PUT_ATT_TEXT(ncid, londeg_id, 'long_name', len_trim(longname), trim(longname)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, londeg_id, 'units', 12, 'degrees_east')
+  if (status .ne. nf_noerr) call handle_err(status, partit)
   !rotated or not rotated due to setting in iceberg module
   description='(un)rotated according to setting of l_geo_out in iceberg module'
-  status = nf90_put_att(ncid, londeg_id, 'description', trim(description)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, londeg_id, 'description', len_trim(description), trim(description)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
   longname='latitude of buoy/iceberg position in degree'
-  status = nf90_put_att(ncid, latdeg_id, 'long_name', trim(longname)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
-  status = nf90_put_att(ncid, latdeg_id, 'units', 'degree')
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_PUT_ATT_TEXT(ncid, latdeg_id, 'long_name', len_trim(longname), trim(longname)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, latdeg_id, 'units', 13, 'degrees_north')
+  if (status .ne. nf_noerr) call handle_err(status, partit)
   !rotated or not rotated due to setting in iceberg module
   description='(un)rotated according to setting of l_geo_out in iceberg module'
-  status = nf90_put_att(ncid, latdeg_id, 'description', trim(description)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, latdeg_id, 'description', len_trim(description), trim(description)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
   longname='status of buoy/iceberg (frozen in/not frozen in)'
-  status = nf90_put_att(ncid, frozen_id, 'long_name', trim(longname)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_PUT_ATT_TEXT(ncid, frozen_id, 'long_name', len_trim(longname), trim(longname)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
   description='1 = frozen, 0 = not frozen, else partially frozen'
-  status = nf90_put_att(ncid, frozen_id, 'description', trim(description)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
-  status = nf90_put_att(ncid, frozen_id, 'Coordinates', 'pos_lon_deg pos_lat_deg') ! arcGIS 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, frozen_id, 'description', len_trim(description), trim(description)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  status = nf_PUT_ATT_TEXT(ncid, frozen_id, 'Coordinates', 23, 'pos_lon_deg pos_lat_deg') ! arcGIS 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
   longname='du/dt of buoy/iceberg in last time step'
-  status = nf90_put_att(ncid, dudt_id, 'long_name', trim(longname)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
-  status = nf90_put_att(ncid, dudt_id, 'units', 'm s^(-2)')
-  if (status .ne. nf90_noerr) call handle_err(status, partit)  
+  status = nf_PUT_ATT_TEXT(ncid, dudt_id, 'long_name', len_trim(longname), trim(longname)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, dudt_id, 'units', 8, 'm s^(-2)')
+  if (status .ne. nf_noerr) call handle_err(status, partit)  
   !rotated or not rotated due to setting in iceberg module
   description='(un)rotated according to setting of l_geo_out in iceberg module'
-  status = nf90_put_att(ncid, dudt_id, 'description', trim(description)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
-  status = nf90_put_att(ncid, dudt_id, 'Coordinates', 'pos_lon_deg pos_lat_deg') ! arcGIS 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, dudt_id, 'description', len_trim(description), trim(description)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  status = nf_PUT_ATT_TEXT(ncid, dudt_id, 'Coordinates', 23, 'pos_lon_deg pos_lat_deg') ! arcGIS 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
   longname='dv/dt of buoy/iceberg in last time step'
-  status = nf90_put_att(ncid, dvdt_id, 'long_name', trim(longname)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
-  status = nf90_put_att(ncid, dvdt_id, 'units', 'm s^(-2)')
-  if (status .ne. nf90_noerr) call handle_err(status, partit)  
+  status = nf_PUT_ATT_TEXT(ncid, dvdt_id, 'long_name', len_trim(longname), trim(longname)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, dvdt_id, 'units', 8, 'm s^(-2)')
+  if (status .ne. nf_noerr) call handle_err(status, partit)  
   !rotated or not rotated due to setting in iceberg module
   description='(un)rotated according to setting of l_geo_out in iceberg module'
-  status = nf90_put_att(ncid, dvdt_id, 'description', trim(description)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
-  status = nf90_put_att(ncid, dvdt_id, 'Coordinates', 'pos_lon_deg pos_lat_deg') ! arcGIS 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, dvdt_id, 'description', len_trim(description), trim(description)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  status = nf_PUT_ATT_TEXT(ncid, dvdt_id, 'Coordinates', 23, 'pos_lon_deg pos_lat_deg') ! arcGIS 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
   longname='velocity of buoy/iceberg, u component'
-  status = nf90_put_att(ncid, uib_id, 'long_name', trim(longname)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
-  status = nf90_put_att(ncid, uib_id, 'units', 'm s^(-1)')
-  if (status .ne. nf90_noerr) call handle_err(status, partit)  
+  status = nf_PUT_ATT_TEXT(ncid, uib_id, 'long_name', len_trim(longname), trim(longname)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, uib_id, 'units', 8, 'm s^(-1)')
+  if (status .ne. nf_noerr) call handle_err(status, partit)  
   !rotated or not rotated due to setting in iceberg module
   description='(un)rotated according to setting of l_geo_out in iceberg module'
-  status = nf90_put_att(ncid, uib_id, 'description', trim(description)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
-  status = nf90_put_att(ncid, uib_id, 'Coordinates', 'pos_lon_deg pos_lat_deg') ! arcGIS 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, uib_id, 'description', len_trim(description), trim(description)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  status = nf_PUT_ATT_TEXT(ncid, uib_id, 'Coordinates', 23, 'pos_lon_deg pos_lat_deg') ! arcGIS 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
   longname='velocity of buoy/iceberg, v component'
-  status = nf90_put_att(ncid, vib_id, 'long_name', trim(longname)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
-  status = nf90_put_att(ncid, vib_id, 'units', 'm s^(-1)')
-  if (status .ne. nf90_noerr) call handle_err(status, partit)  
+  status = nf_PUT_ATT_TEXT(ncid, vib_id, 'long_name', len_trim(longname), trim(longname)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, vib_id, 'units', 8, 'm s^(-1)')
+  if (status .ne. nf_noerr) call handle_err(status, partit)  
   !rotated or not rotated due to setting in iceberg module
   description='(un)rotated according to setting of l_geo_out in iceberg module'
-  status = nf90_put_att(ncid, vib_id, 'description', trim(description)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
-  status = nf90_put_att(ncid, vib_id, 'Coordinates', 'pos_lon_deg pos_lat_deg') ! arcGIS 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, vib_id, 'description', len_trim(description), trim(description)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  status = nf_PUT_ATT_TEXT(ncid, vib_id, 'Coordinates', 23, 'pos_lon_deg pos_lat_deg') ! arcGIS 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
   ! 3 dimensions of the iceberg, comment for buoy case
 
   longname='height of the iceberg'
-  status = nf90_put_att(ncid, height_id, 'long_name', trim(longname)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
-  status = nf90_put_att(ncid, height_id, 'units', 'm')
-  if (status .ne. nf90_noerr) call handle_err(status, partit)  
+  status = nf_PUT_ATT_TEXT(ncid, height_id, 'long_name', len_trim(longname), trim(longname)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, height_id, 'units', 1, 'm')
+  if (status .ne. nf_noerr) call handle_err(status, partit)  
   description='freeboard + draft'
-  status = nf90_put_att(ncid, height_id, 'description', trim(description)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
-  status = nf90_put_att(ncid, height_id, 'Coordinates', 'pos_lon_deg pos_lat_deg') ! arcGIS 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, height_id, 'description', len_trim(description), trim(description)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  status = nf_PUT_ATT_TEXT(ncid, height_id, 'Coordinates', 23, 'pos_lon_deg pos_lat_deg') ! arcGIS 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
   longname='length of the iceberg'
-  status = nf90_put_att(ncid, length_id, 'long_name', trim(longname)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
-  status = nf90_put_att(ncid, length_id, 'units', 'm')
-  if (status .ne. nf90_noerr) call handle_err(status, partit)  
+  status = nf_PUT_ATT_TEXT(ncid, length_id, 'long_name', len_trim(longname), trim(longname)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, length_id, 'units', 1, 'm')
+  if (status .ne. nf_noerr) call handle_err(status, partit)  
   description='open'
-  status = nf90_put_att(ncid, length_id, 'description', trim(description)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
-  status = nf90_put_att(ncid, length_id, 'Coordinates', 'pos_lon_deg pos_lat_deg') ! arcGIS 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, length_id, 'description', len_trim(description), trim(description)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  status = nf_PUT_ATT_TEXT(ncid, length_id, 'Coordinates', 23, 'pos_lon_deg pos_lat_deg') ! arcGIS 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
   longname='width of the iceberg'
-  status = nf90_put_att(ncid, width_id, 'long_name', trim(longname)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
-  status = nf90_put_att(ncid, width_id, 'units', 'm')
-  if (status .ne. nf90_noerr) call handle_err(status, partit)  
+  status = nf_PUT_ATT_TEXT(ncid, width_id, 'long_name', len_trim(longname), trim(longname)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, width_id, 'units', 1, 'm')
+  if (status .ne. nf_noerr) call handle_err(status, partit)  
   description='open'
-  status = nf90_put_att(ncid, width_id, 'description', trim(description)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
-  status = nf90_put_att(ncid, width_id, 'Coordinates', 'pos_lon_deg pos_lat_deg') ! arcGIS 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, width_id, 'description', len_trim(description), trim(description)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  status = nf_PUT_ATT_TEXT(ncid, width_id, 'Coordinates', 23, 'pos_lon_deg pos_lat_deg') ! arcGIS 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
 
   ! 4 additional iceberg variables (meltrates), comment for buoy case
 
   longname='basal volume loss'
-  status = nf90_put_att(ncid, bvl_id, 'long_name', trim(longname)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
-  status = nf90_put_att(ncid, bvl_id, 'units', 'm^3 (ice) day^(-1)')
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_PUT_ATT_TEXT(ncid, bvl_id, 'long_name', len_trim(longname), trim(longname)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, bvl_id, 'units', 18, 'm^3 (ice) day^(-1)')
+  if (status .ne. nf_noerr) call handle_err(status, partit)
   description='losses are averaged over the preceding output interval'
-  status = nf90_put_att(ncid, bvl_id, 'description', trim(description)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
-  status = nf90_put_att(ncid, bvl_id, 'Coordinates', 'pos_lon_deg pos_lat_deg') ! arcGIS 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, bvl_id, 'description', len_trim(description), trim(description)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  status = nf_PUT_ATT_TEXT(ncid, bvl_id, 'Coordinates', 23, 'pos_lon_deg pos_lat_deg') ! arcGIS 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
   longname='lateral volume loss due to 1) bouyant convection'
-  status = nf90_put_att(ncid, lvlv_id, 'long_name', trim(longname)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
-  status = nf90_put_att(ncid, lvlv_id, 'units', 'm^3 (ice) day^(-1)')
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_PUT_ATT_TEXT(ncid, lvlv_id, 'long_name', len_trim(longname), trim(longname)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, lvlv_id, 'units', 18, 'm^3 (ice) day^(-1)')
+  if (status .ne. nf_noerr) call handle_err(status, partit)
   description='losses are averaged over the preceding output interval'
-  status = nf90_put_att(ncid, lvlv_id, 'description', trim(description)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
-  status = nf90_put_att(ncid, lvlv_id, 'Coordinates', 'pos_lon_deg pos_lat_deg') ! arcGIS 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, lvlv_id, 'description', len_trim(description), trim(description)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  status = nf_PUT_ATT_TEXT(ncid, lvlv_id, 'Coordinates', 23, 'pos_lon_deg pos_lat_deg') ! arcGIS 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
   longname='lateral volume loss due to 2) wave erosion'
-  status = nf90_put_att(ncid, lvle_id, 'long_name', trim(longname)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
-  status = nf90_put_att(ncid, lvle_id, 'units', 'm^3 (ice) day^(-1)')
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_PUT_ATT_TEXT(ncid, lvle_id, 'long_name', len_trim(longname), trim(longname)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, lvle_id, 'units', 18, 'm^3 (ice) day^(-1)')
+  if (status .ne. nf_noerr) call handle_err(status, partit)
   description='losses are averaged over the preceding output interval'
-  status = nf90_put_att(ncid, lvle_id, 'description', trim(description)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
-  status = nf90_put_att(ncid, lvle_id, 'Coordinates', 'pos_lon_deg pos_lat_deg') ! arcGIS 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, lvle_id, 'description', len_trim(description), trim(description)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  status = nf_PUT_ATT_TEXT(ncid, lvle_id, 'Coordinates', 23, 'pos_lon_deg pos_lat_deg') ! arcGIS 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
   longname='lateral volume loss due to 3) "basal" formulation'
-  status = nf90_put_att(ncid, lvlb_id, 'long_name', trim(longname)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
-  status = nf90_put_att(ncid, lvlb_id, 'units', 'm^3 (ice) day^(-1)')
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_PUT_ATT_TEXT(ncid, lvlb_id, 'long_name', len_trim(longname), trim(longname)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, lvlb_id, 'units', 18, 'm^3 (ice) day^(-1)')
+  if (status .ne. nf_noerr) call handle_err(status, partit)
   description='losses are averaged over the preceding output interval'
-  status = nf90_put_att(ncid, lvlb_id, 'description', trim(description)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
-  status = nf90_put_att(ncid, lvlb_id, 'Coordinates', 'pos_lon_deg pos_lat_deg') ! arcGIS 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, lvlb_id, 'description', len_trim(description), trim(description)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  status = nf_PUT_ATT_TEXT(ncid, lvlb_id, 'Coordinates', 23, 'pos_lon_deg pos_lat_deg') ! arcGIS 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
   ! LA: add felem
   longname='fesom element'
-  status = nf90_put_att(ncid, felem_id, 'long_name', trim(longname)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
-  status = nf90_put_att(ncid, felem_id, 'units', '')
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_PUT_ATT_TEXT(ncid, felem_id, 'long_name', len_trim(longname), trim(longname)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, felem_id, 'units', 0, '')
+  if (status .ne. nf_noerr) call handle_err(status, partit)
   description=''
-  status = nf90_put_att(ncid, felem_id, 'description', trim(description)) 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
-  status = nf90_put_att(ncid, felem_id, 'Coordinates', 'pos_lon_deg pos_lat_deg') ! arcGIS 
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
-
-  status = nf90_enddef(ncid)
-  if (status .ne. nf90_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, felem_id, 'description', len_trim(description), trim(description)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  status = nf_PUT_ATT_TEXT(ncid, felem_id, 'Coordinates', 23, 'pos_lon_deg pos_lat_deg') ! arcGIS 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+
+  longname='grounded'
+  status = nf_PUT_ATT_TEXT(ncid, grounded_id, 'long_name', len_trim(longname), trim(longname)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  status = nf_put_att_text(ncid, grounded_id, 'units', 0, '')
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  description=''
+  status = nf_put_att_text(ncid, grounded_id, 'description', len_trim(description), trim(description)) 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  status = nf_PUT_ATT_TEXT(ncid, grounded_id, 'Coordinates', 23, 'pos_lon_deg pos_lat_deg') ! arcGIS 
+  if (status .ne. nf_noerr) call handle_err(status, partit)
+  
+  status = nf_enddef(ncid)
+  if (status .ne. nf_noerr) call handle_err(status, partit)
 
-  status = nf90_close(ncid)
-  if (status .ne. nf90_noerr) call handle_err(status, partit) 
+  status=nf_close(ncid)
+  if (status .ne. nf_noerr) call handle_err(status, partit) 
 
  ! initialize the counter for saving results
   save_count_buoys=1
@@ -1948,12 +1933,13 @@ subroutine write_buoy_props_netcdf(partit)
 
  use g_config
  use g_clock
- use netcdf
 ! use iceberg_params, only : buoy_props, file_icb_netcdf, save_count_buoys, prev_sec_in_year, bvl_mean, lvlv_mean, lvle_mean, lvlb_mean
  use g_forcing_param
 
  implicit none
 
+#include "netcdf.inc" 
+
   integer                   :: status,ncid, istep
   integer                   :: dimid_ib, dimid_rec, dimids(2)
   integer                   :: time_varid, iter_varid
@@ -1961,9 +1947,8 @@ subroutine write_buoy_props_netcdf(partit)
   integer                   :: frozen_id, dudt_id, dvdt_id
   integer                   :: uib_id, vib_id  
   integer                   :: height_id, length_id, width_id
-  integer                   :: bvl_id, lvlv_id, lvle_id, lvlb_id, felem_id
+  integer                   :: bvl_id, lvlv_id, lvle_id, lvlb_id, felem_id, grounded_id
   integer                   :: start(2), count(2)
-  integer                   :: start_1d(1), count_1d(1)
   real(kind=8)              :: sec_in_year
 type(t_partit), intent(inout), target :: partit
 !type(t_ice),    intent(inout), target :: ice
@@ -1979,159 +1964,147 @@ subroutine write_buoy_props_netcdf(partit)
      sec_in_year=dt*istep
 
     ! open files
-     status = nf90_open(trim(file_icb_netcdf), nf90_write, ncid)
-     if (status .ne. nf90_noerr) call handle_err(status, partit)
+     status = nf_open(trim(file_icb_netcdf), nf_write, ncid)
+     if (status .ne. nf_noerr) call handle_err(status, partit)
 
      ! inquire variable id
 
-     status = nf90_inq_varid(ncid, 'time', time_varid)
-     if (status .ne. nf90_noerr) call handle_err(status, partit)
+     status = nf_inq_varid(ncid, 'time', time_varid)
+     if (status .ne. nf_noerr) call handle_err(status, partit)
 
-     status=nf90_inq_varid(ncid, 'iter', iter_varid)
-     if (status .ne. nf90_noerr) call handle_err(status, partit)
+     status=nf_inq_varid(ncid, 'iter', iter_varid)
+     if (status .ne. nf_noerr) call handle_err(status, partit)
 
-     status = nf90_inq_varid(ncid, 'pos_lon_rad', lonrad_id)
-     if (status .ne. nf90_noerr) call handle_err(status, partit)
+     status = nf_inq_varid(ncid, 'pos_lon_rad', lonrad_id)
+     if (status .ne. nf_noerr) call handle_err(status, partit)
 
-     status = nf90_inq_varid(ncid, 'pos_lat_rad', latrad_id)
-     if (status .ne. nf90_noerr) call handle_err(status, partit)
+     status = nf_inq_varid(ncid, 'pos_lat_rad', latrad_id)
+     if (status .ne. nf_noerr) call handle_err(status, partit)
 
-     status = nf90_inq_varid(ncid, 'pos_lon_deg', londeg_id)
-     if (status .ne. nf90_noerr) call handle_err(status, partit)
+     status = nf_inq_varid(ncid, 'pos_lon_deg', londeg_id)
+     if (status .ne. nf_noerr) call handle_err(status, partit)
 
-     status = nf90_inq_varid(ncid, 'pos_lat_deg', latdeg_id)
-     if (status .ne. nf90_noerr) call handle_err(status, partit)
+     status = nf_inq_varid(ncid, 'pos_lat_deg', latdeg_id)
+     if (status .ne. nf_noerr) call handle_err(status, partit)
 
-     status = nf90_inq_varid(ncid, 'frozen_in',  frozen_id)
-     if (status .ne. nf90_noerr) call handle_err(status, partit)
+     status = nf_inq_varid(ncid, 'frozen_in',  frozen_id)
+     if (status .ne. nf_noerr) call handle_err(status, partit)
 
-     status = nf90_inq_varid(ncid, 'du_dt',  dudt_id)
-     if (status .ne. nf90_noerr) call handle_err(status, partit)
+     status = nf_inq_varid(ncid, 'du_dt',  dudt_id)
+     if (status .ne. nf_noerr) call handle_err(status, partit)
 
-     status = nf90_inq_varid(ncid, 'dv_dt',  dvdt_id)
-     if (status .ne. nf90_noerr) call handle_err(status, partit)
+     status = nf_inq_varid(ncid, 'dv_dt',  dvdt_id)
+     if (status .ne. nf_noerr) call handle_err(status, partit)
 
-     status = nf90_inq_varid(ncid, 'icb_vel_u',  uib_id)
-     if (status .ne. nf90_noerr) call handle_err(status, partit)
+     status = nf_inq_varid(ncid, 'icb_vel_u',  uib_id)
+     if (status .ne. nf_noerr) call handle_err(status, partit)
 
-     status = nf90_inq_varid(ncid, 'icb_vel_v',  vib_id)
-     if (status .ne. nf90_noerr) call handle_err(status, partit)  
+     status = nf_inq_varid(ncid, 'icb_vel_v',  vib_id)
+     if (status .ne. nf_noerr) call handle_err(status, partit)  
 
      ! inquire 3 additional IDs for iceberg case, comment for buoy case
 
-     status = nf90_inq_varid(ncid, 'height',  height_id)
-     if (status .ne. nf90_noerr) call handle_err(status, partit)
+     status = nf_inq_varid(ncid, 'height',  height_id)
+     if (status .ne. nf_noerr) call handle_err(status, partit)
 
-     status = nf90_inq_varid(ncid, 'length',  length_id)
-     if (status .ne. nf90_noerr) call handle_err(status, partit)
+     status = nf_inq_varid(ncid, 'length',  length_id)
+     if (status .ne. nf_noerr) call handle_err(status, partit)
 
-     status = nf90_inq_varid(ncid, 'width',  width_id)
-     if (status .ne. nf90_noerr) call handle_err(status, partit)  
+     status = nf_inq_varid(ncid, 'width',  width_id)
+     if (status .ne. nf_noerr) call handle_err(status, partit)  
 
      ! inquire 4 additional IDs for iceberg case, comment for buoy case
 
-     status = nf90_inq_varid(ncid, 'bvl',  bvl_id)
-     if (status .ne. nf90_noerr) call handle_err(status, partit)
+     status = nf_inq_varid(ncid, 'bvl',  bvl_id)
+     if (status .ne. nf_noerr) call handle_err(status, partit)
 
-     status = nf90_inq_varid(ncid, 'lvlv',  lvlv_id)
-     if (status .ne. nf90_noerr) call handle_err(status, partit)
+     status = nf_inq_varid(ncid, 'lvlv',  lvlv_id)
+     if (status .ne. nf_noerr) call handle_err(status, partit)
 
-     status = nf90_inq_varid(ncid, 'lvle',  lvle_id)
-     if (status .ne. nf90_noerr) call handle_err(status, partit)
+     status = nf_inq_varid(ncid, 'lvle',  lvle_id)
+     if (status .ne. nf_noerr) call handle_err(status, partit)
 
-     status = nf90_inq_varid(ncid, 'lvlb',  lvlb_id)
-     if (status .ne. nf90_noerr) call handle_err(status, partit) 
+     status = nf_inq_varid(ncid, 'lvlb',  lvlb_id)
+     if (status .ne. nf_noerr) call handle_err(status, partit) 
 
      ! * LA: include fesom elemt in output
-     status = nf90_inq_varid(ncid, 'felem', felem_id)
-     if (status .ne. nf90_noerr) call handle_err(status, partit) 
-
-     !buoy_props(ib, 1) = lon_rad_out
-     !buoy_props(ib, 2) = lat_rad_out
-     !buoy_props(ib, 3) = lon_deg_out
-     !buoy_props(ib, 4) = lat_deg_out
-     !buoy_props(ib, 5) = frozen_in
-     !buoy_props(ib, 6) = dudt_out
-     !buoy_props(ib, 7) = dvdt_out
-     !buoy_props(ib, 8) = u_ib_out
-     !buoy_props(ib, 9) = v_ib_out
-     !buoy_props(ib,10) = height_ib
-     !buoy_props(ib,11) = length_ib 
-     !buoy_props(ib,12) = width_ib
-
-     !bvl_mean(ib)*step_per_day
-     !lvlv_mean(ib)*step_per_day
-     !lvle_mean(ib)*step_per_day
-     !lvlb_mean(ib)*step_per_day
-
-      ! time and iteration
-      status=nf90_put_var(ncid, time_varid, prev_sec_in_year+sec_in_year) 
-      if (status .ne. nf90_noerr) call handle_err(status, partit)
-      status=nf90_put_var(ncid, iter_varid, istep)
-      if (status .ne. nf90_noerr) call handle_err(status, partit)
-
-      !variables
-      start=(/1,save_count_buoys/)
-      count=(/ib_num, 1/)
-      status=nf90_put_var(ncid, lonrad_id, buoy_props(:, 1), start=start, count=count) 
-      if (status .ne. nf90_noerr) call handle_err(status, partit)
-
-     status=nf90_put_var(ncid, latrad_id, buoy_props(:, 2), start=start, count=count) 
-     if (status .ne. nf90_noerr) call handle_err(status, partit)
-
-     status=nf90_put_var(ncid, londeg_id, buoy_props(:, 3), start=start, count=count) 
-     if (status .ne. nf90_noerr) call handle_err(status, partit)
-
-     status=nf90_put_var(ncid, latdeg_id, buoy_props(:, 4), start=start, count=count) 
-     if (status .ne. nf90_noerr) call handle_err(status, partit)
-
-     status=nf90_put_var(ncid, frozen_id, buoy_props(:, 5), start=start, count=count) 
-     if (status .ne. nf90_noerr) call handle_err(status, partit)
-
-     status=nf90_put_var(ncid, dudt_id, buoy_props(:, 6), start=start, count=count) 
-     if (status .ne. nf90_noerr) call handle_err(status, partit)
-
-     status=nf90_put_var(ncid, dvdt_id, buoy_props(:, 7), start=start, count=count) 
-     if (status .ne. nf90_noerr) call handle_err(status, partit)
-
-     status=nf90_put_var(ncid, uib_id, buoy_props(:, 8), start=start, count=count) 
-     if (status .ne. nf90_noerr) call handle_err(status, partit)
-
-     status=nf90_put_var(ncid, vib_id, buoy_props(:, 9), start=start, count=count) 
-     if (status .ne. nf90_noerr) call handle_err(status, partit)
+     status = nf_inq_varid(ncid, 'felem', felem_id)
+     if (status .ne. nf_noerr) call handle_err(status, partit) 
+
+     status = nf_inq_varid(ncid, 'grounded', grounded_id)
+     if (status .ne. nf_noerr) call handle_err(status, partit) 
+
+     ! time and iteration
+     status=nf_put_vara_double(ncid, time_varid, save_count_buoys, 1, prev_sec_in_year+sec_in_year) 
+     if (status .ne. nf_noerr) call handle_err(status, partit)
+     status=nf_put_vara_int(ncid, iter_varid, save_count_buoys, 1, istep)
+     if (status .ne. nf_noerr) call handle_err(status, partit)
+
+     !variables
+     start=(/1,save_count_buoys/)
+     count=(/ib_num, 1/)
+     status=nf_put_vara_double(ncid, lonrad_id, start, count, buoy_props(:, 1)) 
+     if (status .ne. nf_noerr) call handle_err(status, partit)
+
+     status=nf_put_vara_double(ncid, latrad_id, start, count, buoy_props(:, 2)) 
+     if (status .ne. nf_noerr) call handle_err(status, partit)
+
+     status=nf_put_vara_double(ncid, londeg_id, start, count, buoy_props(:, 3)) 
+     if (status .ne. nf_noerr) call handle_err(status, partit)
+
+     status=nf_put_vara_double(ncid, latdeg_id, start, count, buoy_props(:, 4)) 
+     if (status .ne. nf_noerr) call handle_err(status, partit)
+
+     status=nf_put_vara_double(ncid, frozen_id, start, count, buoy_props(:, 5)) 
+     if (status .ne. nf_noerr) call handle_err(status, partit)
+
+     status=nf_put_vara_double(ncid, dudt_id, start, count, buoy_props(:, 6)) 
+     if (status .ne. nf_noerr) call handle_err(status, partit)
+
+     status=nf_put_vara_double(ncid, dvdt_id, start, count, buoy_props(:, 7)) 
+     if (status .ne. nf_noerr) call handle_err(status, partit)
+
+     status=nf_put_vara_double(ncid, uib_id, start, count, buoy_props(:, 8)) 
+     if (status .ne. nf_noerr) call handle_err(status, partit)
+
+     status=nf_put_vara_double(ncid, vib_id, start, count, buoy_props(:, 9)) 
+     if (status .ne. nf_noerr) call handle_err(status, partit)
 
      ! write 3 additional variables for iceberg case, comment for buoy case
 
-     status=nf90_put_var(ncid, height_id, buoy_props(:,10), start=start, count=count) 
-     if (status .ne. nf90_noerr) call handle_err(status, partit)
+     status=nf_put_vara_double(ncid, height_id, start, count, buoy_props(:,10)) 
+     if (status .ne. nf_noerr) call handle_err(status, partit)
 
-     status=nf90_put_var(ncid, length_id, buoy_props(:,11), start=start, count=count) 
-     if (status .ne. nf90_noerr) call handle_err(status, partit)
+     status=nf_put_vara_double(ncid, length_id, start, count, buoy_props(:,11)) 
+     if (status .ne. nf_noerr) call handle_err(status, partit)
 
-     status=nf90_put_var(ncid, width_id, buoy_props(:,12), start=start, count=count) 
-     if (status .ne. nf90_noerr) call handle_err(status, partit)
+     status=nf_put_vara_double(ncid, width_id, start, count, buoy_props(:,12)) 
+     if (status .ne. nf_noerr) call handle_err(status, partit)
 
      ! write 4 additional variables for iceberg case, comment for buoy case
 
-     status=nf90_put_var(ncid, bvl_id, bvl_mean(:)*step_per_day, start=start, count=count) 
-     if (status .ne. nf90_noerr) call handle_err(status, partit)
+     status=nf_put_vara_double(ncid, bvl_id, start, count, bvl_mean(:)*step_per_day) 
+     if (status .ne. nf_noerr) call handle_err(status, partit)
 
-     status=nf90_put_var(ncid, lvlv_id, lvlv_mean(:)*step_per_day, start=start, count=count) 
-     if (status .ne. nf90_noerr) call handle_err(status, partit)
+     status=nf_put_vara_double(ncid, lvlv_id, start, count, lvlv_mean(:)*step_per_day) 
+     if (status .ne. nf_noerr) call handle_err(status, partit)
 
-     status=nf90_put_var(ncid, lvle_id, lvle_mean(:)*step_per_day, start=start, count=count) 
-     if (status .ne. nf90_noerr) call handle_err(status, partit)
+     status=nf_put_vara_double(ncid, lvle_id, start, count, lvle_mean(:)*step_per_day) 
+     if (status .ne. nf_noerr) call handle_err(status, partit)
 
-     status=nf90_put_var(ncid, lvlb_id, lvlb_mean(:)*step_per_day, start=start, count=count) 
-     if (status .ne. nf90_noerr) call handle_err(status, partit)
+     status=nf_put_vara_double(ncid, lvlb_id, start, count, lvlb_mean(:)*step_per_day) 
+     if (status .ne. nf_noerr) call handle_err(status, partit)
 
      ! LA: add felem
-     status=nf90_put_var(ncid, felem_id, buoy_props(:,13), start=start, count=count) 
-     if (status .ne. nf90_noerr) call handle_err(status, partit)
+     status=nf_put_vara_double(ncid, felem_id, start, count, buoy_props(:,13)) 
+     if (status .ne. nf_noerr) call handle_err(status, partit)
 
+     status=nf_put_vara_double(ncid, grounded_id, start, count, buoy_props(:,14)) 
+     if (status .ne. nf_noerr) call handle_err(status, partit)
+     
      !close file
-     status=nf90_close(ncid)
-     if (status .ne. nf90_noerr) call handle_err(status, partit)
+     status=nf_close(ncid)
+     if (status .ne. nf_noerr) call handle_err(status, partit)
 
      save_count_buoys=save_count_buoys+1
 !==========================================================
diff --git a/src/icb_thermo.F90 b/src/icb_thermo.F90
index 70a42d6e0..b2cf375f4 100644
--- a/src/icb_thermo.F90
+++ b/src/icb_thermo.F90
@@ -91,36 +91,37 @@ subroutine iceberg_meltrates(partit, mesh, M_b, M_v, M_e, M_bv, &
   dz_acc = 0.0
 
 
+  ! Zero out per-level heat flux arrays before filling up to ib_n_lvls.
+  ! prepare_icb2fesom accesses these up to idx_d(i) which can exceed ib_n_lvls
+  ! for element nodes deeper than the shallowest node; stale values from the
+  ! previous iceberg step would otherwise inject spurious heat into the ocean.
+  hfl_flux_ib(ib,:)  = 0.0
+  hfbv_flux_ib(ib,:) = 0.0
+
   n2=elem2D_nodes(1,elem)
-  
-  ! iterate over all layers that contain the iceberg
   do n=1,ib_n_lvls
-    
-    ! 3-eq. formulation for lateral 'basal' melting [m/s]
-    lev_up  = mesh%zbar_3d_n(n, n2)                 ! upper level
-    if( n==nlevels_nod2D(n2) ) then                 ! if level is lowest level ...
-        lev_low = mesh%zbar_n_bot(n2)               ! ... lower level is bottom topography
+  !3-eq. formulation for lateral 'basal' melting [m/s]
+    lev_up  = mesh%zbar_3d_n(n, n2)
+    if( n==nlevels_nod2D(n2) ) then
+        lev_low = mesh%zbar_n_bot(n2)
     else
-        lev_low = mesh%zbar_3d_n(n+1, n2)           ! ... otherwise, lower level is one level below upper level
+        lev_low = mesh%zbar_3d_n(n+1, n2)
     end if
-   
-    ! assign dz for vertical integration
-    if( abs(lev_low)>=abs(depth_ib) ) then          ! if lower level below iceberg base ...
-        dz = abs(lev_up - depth_ib)                 ! ... dz is equal to difference between upper level and iceberg base
-    elseif(lev_low == lev_up) then                  ! if lower level equal to upper level ... (why should this happen?)
-        exit                                        ! ... exit
+    
+    if( abs(lev_low)>=abs(depth_ib) ) then !.AND. (abs(lev_up)<=abs(depth_ib)) ) then
+        dz = abs(lev_up - depth_ib)
+    elseif(lev_low == lev_up) then
+        exit
     else
-        dz = abs(lev_low - lev_up)                  ! ... otherwise, dz is equal to difference between lower level and upper level
+        dz = abs(lev_low - lev_up)
     end if
-    dz_acc = dz_acc + dz                            ! sum up dz along iceberg depth
     
-    v_ibmino  = sqrt( (u_ib - arr_uo_ib(n))**2 + (v_ib - arr_vo_ib(n))**2 )         ! depth-average rel. velocity
-
-    ! calculate freshwater and heat flux densities ...
-    ! why dz_acc+dz/2.0? why not dz_acc-dz/2.0?
-    call iceberg_heat_water_fluxes_3eq(partit, ib, M_bv, H_bv, arr_T_ave_ib(n), arr_S_ave_ib(n), v_ibmino, dz_acc-dz/2.0, tf)
-    M_bv_dz = M_bv_dz + M_bv*dz                     ! ... integrating freshwater flux along iceberg depth
-    hfbv_flux_ib(ib,n) = H_bv * (2*length_ib*dz  + 2*length_ib*dz ) * scaling(ib)   ! ... assign heat flux bv on for layer n 
+    v_ibmino  = sqrt( (u_ib - arr_uo_ib(n))**2 + (v_ib - arr_vo_ib(n))**2 ) ! depth-average rel. velocity
+    call iceberg_heat_water_fluxes_3eq(partit, ib, M_bv, H_bv, arr_T_ave_ib(n), arr_S_ave_ib(n),v_ibmino, -(dz_acc+dz/2.0), tf)
+    dz_acc = dz_acc + dz
+    M_bv_dz = M_bv_dz + M_bv*dz
+    
+    hfbv_flux_ib(ib,n) = H_bv * (2*length_ib*dz  + 2*length_ib*dz ) * scaling(ib)
   
     !'thermal driving', defined as the elevation of ambient water 
     !temperature above freezing point' (Neshyba and Josberger, 1979).
@@ -142,8 +143,15 @@ subroutine iceberg_meltrates(partit, mesh, M_b, M_v, M_e, M_bv, &
     hfl_flux_ib(ib,n) = H_v * (2*length_ib*dz  + 2*length_ib*dz ) * scaling(ib)
     !fwl_flux_ib = M_v
   end do
-  M_bv  = M_bv_dz / abs(depth_ib)
-  M_v   = M_v_dz / abs(depth_ib)
+  ! Guard against depth_ib==0 (should not happen if volume check is passed,
+  ! but fp edge cases and cavity geometry can produce near-zero depth).
+  if (abs(depth_ib) > 0.0) then
+    M_bv  = M_bv_dz / abs(depth_ib)
+    M_v   = M_v_dz / abs(depth_ib)
+  else
+    M_bv  = 0.0
+    M_v   = 0.0
+  end if
 
   !wave erosion
   absamino = sqrt( (ua_ib - uo_ib)**2 + (va_ib - vo_ib)**2 )
@@ -151,6 +159,11 @@ subroutine iceberg_meltrates(partit, mesh, M_b, M_v, M_e, M_bv, &
   damping = 0.5 * (1.0 + cos(conci_ib**3 * Pi))
   M_e = 1./6. * sea_state * (sst_ib + 2.0) * damping
   M_e = M_e/86400.
+  ! The Bigg/Silva wave erosion formula is only defined for SST > -2 degC.
+  ! When SST < -2 degC the formula produces M_e < 0 (nonphysical regrowth),
+  ! which reverses the sign of the FW and heat flux delivered to the ocean,
+  ! creating a runaway cooling feedback via the ALE bc_surface term. Clamp here.
+  if (M_e < 0.0) M_e = 0.0
   H_e = M_e * rho_icb * L
   
   ! check wave erosion potential
@@ -201,6 +214,7 @@ subroutine iceberg_newdimensions(partit, ib, depth_ib,height_ib,length_ib,width_
   character, intent(IN)	:: file_meltrates*80
   
   real			:: dh_b, dh_v, dh_e, dh_bv, bvl, lvl_b, lvl_v, lvl_e, tvl, volume_before, volume_after
+  real			:: hf_scale
   integer		:: icbID
   logical		:: force_last_output
   real, dimension(4)	:: arr
@@ -239,6 +253,12 @@ subroutine iceberg_newdimensions(partit, ib, depth_ib,height_ib,length_ib,width_
     volume_before=height_ib*length_ib*width_ib
 
     if((tvl .ge. volume_before) .OR. (volume_before .le. smallestvol_icb)) then
+        ! Scale heat fluxes when melt exceeds remaining volume
+        hf_scale = volume_before / max(tvl, 1.0e-30)
+        hfb_flux_ib(ib)    = hfb_flux_ib(ib)    * hf_scale
+        hfe_flux_ib(ib)    = hfe_flux_ib(ib)    * hf_scale
+        hfbv_flux_ib(ib,:) = hfbv_flux_ib(ib,:) * hf_scale
+        hfl_flux_ib(ib,:)  = hfl_flux_ib(ib,:)  * hf_scale
     	volume_after=0.0    	
 	    depth_ib = 0.0
     	height_ib= 0.0
@@ -269,6 +289,8 @@ subroutine iceberg_newdimensions(partit, ib, depth_ib,height_ib,length_ib,width_
 
         !iceberg smaller than critical value after melting?
         if (volume_after .le. smallestvol_icb) then
+            ! Add latent heat for remnant volume to erosion heat flux (bug 5)
+            hfe_flux_ib(ib) = hfe_flux_ib(ib) + (volume_before - tvl) * rho_icb * L / dt / REAL(steps_per_ib_step) * scaling(ib)
             volume_after=0.0    	
 	        depth_ib = 0.0
     	    height_ib= 0.0
@@ -447,6 +469,9 @@ subroutine iceberg_heat_water_fluxes_3eq(partit, ib, M_b, H_b, T_ib,S_ib,v_rel,
      gat  = gats1/(gats2+12.5*pr1)
      gas  = gats1/(gats2+12.5*sc1)
 
+     !RG3417 gat  = 1.00e-4   ![m/s]  RT: to be replaced by velocity-dependent equations later
+     !RG3417 gas  = 5.05e-7   ![m/s]  RT: to be replaced by velocity-dependent equations later
+
      ! Calculate
      ! density in the boundary layer: rhow
      ! and interface pressure pg [dbar],
@@ -463,10 +488,23 @@ subroutine iceberg_heat_water_fluxes_3eq(partit, ib, M_b, H_b, T_ib,S_ib,v_rel,
      ep1 = cpw*gat
      ep2 = cpi*gas
      ep3 = lhf*gas
-     ep31 = -rhor*cpi*tdif/zice   !RG4190 / RG44027
+     ! ep31 = -rhor*cpi*tdif/zice is only used in the freezing branch.
+     ! Guard against zice==0 (would be Inf, giving NaN in the quadratic).
+     if (abs(zice) > 0.0) then
+       ep31 = -rhor*cpi*tdif/zice   !RG4190 / RG44027
+     else
+       ep31 = 0.0
+     end if
      ep4 = b+c*zice
      ep5 = gas/rhor
 
+
+!rt RG4190     ! negative heat flux term in the ice (due to -kappa/D)
+!rt RG4190     ex1 = a*(ep1-ep2)
+!rt RG4190     ex2 = ep1*(ep4-tin)+ep2*(tob+a*sal-ep4)-ep3
+!rt RG4190     ex3 = sal*(ep2*(ep4-tob)+ep3)
+
+
 !RT RG4190/RG44027:
 !    In case of melting ice account for changing temperature gradient, i.e. switch from heat conduction to heat capacity approach
 !TR What to do in iceberg case? LEAVE AS IT IS
@@ -506,16 +544,43 @@ subroutine iceberg_heat_water_fluxes_3eq(partit, ib, M_b, H_b, T_ib,S_ib,v_rel,
         sf = sf2
      endif
 
-     t_freeze = tf                                      ! output of freezing temperature
-     
-     H_b = rhow*cpw*gat*(tin-tf)                        ! heat flux [W] positive for upward
-     M_b = gas*(sf-sal)/sf                              ! freshwater flux [m/s] freshwater per second
-     M_b = - (rhow / rhoi) * M_b 		                ! freshwater flux [m (ice) per second], NOW positive for melting
+     t_freeze = tf ! output of freezing temperature
+     ! Calculate the melting/freezing rate [m/s]
+     ! seta = ep5*(1.0-sal/sf)     !rt thinks this is not needed; TR: Why different to M_b? LIQUID vs. ICE
+
+     !rt  t_surf_flux(i,j)=gat*(tf-tin)
+     !rt  s_surf_flux(i,j)=gas*(sf-(s(i,j,N,lrhs)+35.0))
+
+     !hfb_flux_ib(ib)  = rhow*cpw*gat*(tin-tf)*scaling(ib)      ! [W/m2]  ! positive for upward
+     !hfb_flux_ib(ib)  = rhow*cpw*gat*(tin-tf)*length_ib(ib)*width_ib(ib)*scaling(ib)      ! [W]  ! positive for upward
+     H_b  = rhow*cpw*gat*(tin-tf) !*length_ib(ib)*width_ib(ib)*scaling(ib)      ! [W]  ! positive for upward
+
+     !fw_flux_ib(ib) =          gas*(sf-sal)/sf   ! [m/s]   !
+     M_b 	    =          gas*(sf-sal)/sf   ! [m/s]   ! m freshwater per second
+     !fw_flux_ib(ib) = M_b
+     !fw = -M_b
+     M_b = - (rhow / rhoi) * M_b 		 ! [m (ice) per second], positive for melting? NOW positive for melting
 
-     ! avoid basal freezing for grounded icebergs
+     !LA avoid basal freezing for grounded icebergs
      if(grounded(ib) .and. (M_b.lt.0.)) then
          M_b = 0.0
+         H_b = 0.0
      endif
+
+     !      qo=-rhor*seta*oofw
+     !      if(seta.le.0.) then
+     !         qc=rhor*seta*hemw
+     !         qo=rhor*seta*oomw
+     !      endif
+
+     ! write(*,'(a10,i10,9f10.3)') 'ice shelf',n,zice,rhow,temp,sal,tin,tf,sf,heat_flux(n),water_flux(n)*86400.*365.
+
+     !for saving to output:
+     !net_heat_flux(n)=-heat_flux(n)      ! positive down
+     !fresh_wa_flux(n)=-water_flux(n)     ! m freshwater per second
+
+  !enddo
+
 end subroutine iceberg_heat_water_fluxes_3eq
 
 subroutine potit_ib(ib,salz,pt,pres,rfpres,tin)
diff --git a/src/oce_ale_tracer.F90 b/src/oce_ale_tracer.F90
index e39aa5fe4..9c4fb7a37 100644
--- a/src/oce_ale_tracer.F90
+++ b/src/oce_ale_tracer.F90
@@ -1203,9 +1203,7 @@ subroutine diff_ver_part_impl_ale(tr_num, dynamics, tracers, ice, partit, mesh)
         if (use_icebergs .and. (.not. turn_off_hf) .and. tracers%data(tr_num)%ID==1) then
             do nz=nzmin, nzmax-1
                 zinv=1.0_WP*dt  !/(zbar(nz)-zbar(nz+1)) ale!
-                !!PS tr(nz)=tr(nz)+(sw_3d(nz, n)-sw_3d(nz+1, n) * ( area(nz+1,n)/areasvol(nz,n)) ) * zinv
-                tr(nz)=tr(nz)+(ibhf_n(nz, n)) * zinv / vcpw * area(nz+1,n)/areasvol(nz,n) !) * zinv / vcpw
-                !tr(nz)=tr(nz)+(ibhf_n(nz, n)-ibhf_n(nz+1, n)) * zinv / vcpw ! * area(nz+1,n)/areasvol(nz,n)) * zinv / vcpw
+                tr(nz)=tr(nz) + ibhf_n(nz, n) * zinv / vcpw * (area(nz,n)/areasvol(nz,n))
             end do
         end if