Origin/tmp/apply patches for icb v3

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
         !!!!!!!!!!!!!!!!!!
 
 

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.

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