Merge pull request #649 from FESOM/fix/wave_erosion

Fix/wave erosion

src/icb_allocate.F90

@@ -26,6 +26,10 @@ subroutine allocate_icb(partit, mesh)
   allocate(ibhf_n(mesh%nl, n2))
   ibhf_n    = 0.0_WP
 
+  allocate(wave_erosion_potential(elem2D), linit_wave_erosion_pot(elem2D))
+  wave_erosion_potential = 0.0
+  linit_wave_erosion_pot = .true.
+
   allocate(calving_day(ib_num))
   calving_day = 1   !28.0: September 29 for restart in 1 SEP 97 ! 271.0: September 29 for year 1997
   allocate(height_ib(ib_num))

src/icb_coupling.F90

@@ -19,6 +19,8 @@ subroutine reset_ib_fluxes()
     use g_config
     use iceberg_params
 
+    wave_erosion_potential  = 0.0
+    linit_wave_erosion_pot  = .true.
     ibfwbv  = 0.0
     ibfwb   = 0.0
     ibfwl   = 0.0
@@ -124,14 +126,14 @@ subroutine prepare_icb2fesom(mesh, partit, ib,i_have_element,localelement,depth_
                         dz = abs(abs(lev_up) - abs(depth_ib))
                     end if              
 
-                    if( depth_ib==0.0 ) then
+                    if( abs(depth_ib) > 0.0 ) then
+                    !    ibhf_n(j,iceberg_node) = ibhf_n(j,iceberg_node) & 
+                    !                                - (hfbv_flux_ib(ib,j)+hfl_flux_ib(ib,j)) &
+                    !                                / tot_area_nods_in_ib_elem(j)
+                    !else
                         ibhf_n(j,iceberg_node) = ibhf_n(j,iceberg_node) & 
-                                                    - (hfbv_flux_ib(ib,j)+hfl_flux_ib(ib,j)) &
-                                                    / tot_area_nods_in_ib_elem(j)
-                    else
-                        ibhf_n(j,iceberg_node) = ibhf_n(j,iceberg_node) & 
-                                                    - ((hfbv_flux_ib(ib,j)+hfl_flux_ib(ib,j)) * (dz / abs(depth_ib)) & 
-                                                    + hfe_flux_ib(ib) * (dz / abs(height_ib_single))) &
+                                                    - ((hfbv_flux_ib(ib,j)+hfl_flux_ib(ib,j)) * (dz / abs(depth_ib))) & 
+                                                    !+ hfe_flux_ib(ib) * (dz / abs(height_ib_single))) &
                                                     / tot_area_nods_in_ib_elem(j)
                     end if
                 end do
@@ -145,7 +147,8 @@ subroutine prepare_icb2fesom(mesh, partit, ib,i_have_element,localelement,depth_
 
                 if( height_ib_single .ne. 0.0 ) then
                     ibhf_n(1,iceberg_node) = ibhf_n(1,iceberg_node) - hfe_flux_ib(ib) & 
-                            * ((abs(height_ib_single)-abs(depth_ib))/abs(height_ib_single)) & 
+                            !* 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

src/icb_dyn.F90

@@ -159,7 +159,7 @@ subroutine iceberg_dyn(mesh, partit, ice, dynamics, ib, new_u_ib, new_v_ib, u_ib
  !========================THERMODYNAMICS============================
  if(l_melt) then
   call FEM_eval(mesh, partit, sst_ib,sss_ib,lon,lat,Tsurf_ib,Ssurf_ib,iceberg_elem)
- 
+
   call iceberg_meltrates(partit, mesh, M_b, M_v, M_e, M_bv, &
 				u_ib,v_ib, arr_uo_ib,arr_vo_ib, ua_ib,va_ib, &
 				sst_ib, length_ib, conci_ib, &

src/icb_modules.F90

@@ -97,7 +97,9 @@ module iceberg_params
   real,dimension(:), allocatable:: fwe_flux_ib, fwl_flux_ib, fwb_flux_ib, fwbv_flux_ib
   real,dimension(:), allocatable:: hfe_flux_ib, hfb_flux_ib, lhfb_flux_ib
   real,dimension(:,:), allocatable:: hfl_flux_ib, hfbv_flux_ib
-
+  real,dimension(:), allocatable:: wave_erosion_potential
+  logical,dimension(:), allocatable:: linit_wave_erosion_pot
+
   !===== FRESHWATER AND HEAT ARRAYS ON FESOM GRID =====
   real,dimension(:), allocatable:: ibhf    !icb heat flux into ocean 
   real(kind=WP),dimension(:,:), allocatable:: ibhf_n    !icb heat flux into ocean

src/icb_thermo.F90

@@ -49,7 +49,7 @@ subroutine iceberg_meltrates(partit, mesh, M_b, M_v, M_e, M_bv, &
   use g_forcing_arrays
   use g_rotate_grid
 
-  use iceberg_params, only: fwe_flux_ib, fwl_flux_ib, fwb_flux_ib, fwbv_flux_ib, hfb_flux_ib, hfbv_flux_ib, hfe_flux_ib, hfl_flux_ib, scaling
+  use iceberg_params, only: fwe_flux_ib, fwl_flux_ib, fwb_flux_ib, fwbv_flux_ib, hfb_flux_ib, hfbv_flux_ib, hfe_flux_ib, hfl_flux_ib, scaling, wave_erosion_potential, linit_wave_erosion_pot
   implicit none
 
   ! LA: include latent heat 2023-04-04
@@ -71,6 +71,7 @@ subroutine iceberg_meltrates(partit, mesh, M_b, M_v, M_e, M_bv, &
   real              :: lev_up, lev_low, dz
   real			:: absamino, damping, sea_state, v_ibmino
   real			:: tf, T_d 				!freezing temp. and 'thermal driving'
+  real          :: ib_wave_erosion_pot
 type(t_mesh), intent(in) , target :: mesh
 type(t_partit), intent(inout), target :: partit
 !==================== MODULES & DECLARATIONS ==========================!= 
@@ -119,7 +120,11 @@ subroutine iceberg_meltrates(partit, mesh, M_b, M_v, M_e, M_bv, &
     !temperature above freezing point' (Neshyba and Josberger, 1979).
     T_d = arr_T_ave_ib(n) - tf
     if(T_d < 0.) T_d = 0.
-
+
+    if (linit_wave_erosion_pot(elem)) then 
+        wave_erosion_potential(elem) = vcpw * T_d * elem_area(elem) * helem(1,elem)
+        linit_wave_erosion_pot(elem) = .false.
+    end if
     !lateral melt (buoyant convection)
     !M_v = 0.00762 * sst_ib + 0.00129 * sst_ib**2
     !M_v = M_v/86400.
@@ -141,11 +146,25 @@ subroutine iceberg_meltrates(partit, mesh, M_b, M_v, M_e, M_bv, &
   M_e = 1./6. * sea_state * (sst_ib + 2.0) * damping
   M_e = M_e/86400.
   H_e = M_e * rho_icb * L
-  hfe_flux_ib(ib) = H_e * (length_ib*abs(height_ib)  + length_ib*abs(height_ib) ) * scaling(ib)
-  !fwe_flux_ib = M_e  
+
+  ! check wave erosion potential
+  ib_wave_erosion_pot = abs(H_e * dt * REAL(steps_per_ib_step) * 2 * length_ib * abs(height_ib) * scaling(ib))
+  if (ib_wave_erosion_pot .le. wave_erosion_potential(elem)) then
+    !if (lverbose_icb) write(*,*) " * wave erosion for ib ", ib, "; wave_erosion_potential=",wave_erosion_potential(elem), "; ib_wave_erosion_pot=",ib_wave_erosion_pot,"; elem=",elem
+    wave_erosion_potential(elem) = wave_erosion_potential(elem) - ib_wave_erosion_pot
+  else
+    !if (lverbose_icb) write(*,*) " * no wave erosion for ib ", ib, "; wave_erosion_potential=",wave_erosion_potential(elem), "; ib_wave_erosion_pot=",ib_wave_erosion_pot,"; elem=",elem
+    H_e = H_e * abs(wave_erosion_potential(elem) / ib_wave_erosion_pot)
+    M_e = M_e * abs(wave_erosion_potential(elem) / ib_wave_erosion_pot)
+    wave_erosion_potential(elem) = 0.0
+!  else
+!    if (lverbose_icb) write(*,*) " * no wave erosion for ib ", ib, "; wave_erosion_potential=",wave_erosion_potential(elem), "; ib_wave_erosion_pot=",ib_wave_erosion_pot
+!    H_e = 0.0
+!    M_e = 0.0
+  end if
+  hfe_flux_ib(ib) = H_e * (length_ib*abs(height_ib) + length_ib*abs(height_ib) ) * scaling(ib)
 end subroutine iceberg_meltrates
 
-
  !***************************************************************************************************************************
  !***************************************************************************************************************************
 

src/write_step_info.F90

@@ -49,7 +49,7 @@ subroutine check_blowup(istep, ice, dynamics, tracers, partit, mesh)
 !
 !===============================================================================
 subroutine write_step_info(istep, outfreq, ice, dynamics, tracers, partit, mesh)
-  use g_config, only: dt, use_ice
+  use g_config, only: dt, use_ice, use_icebergs, ib_num
   use MOD_MESH
   USE MOD_PARTIT
   USE MOD_PARSUP
@@ -61,6 +61,7 @@ subroutine write_step_info(istep, outfreq, ice, dynamics, tracers, partit, mesh)
                  pgf_x, pgf_y, Av, Kv
   use g_comm_auto
   use g_support
+  use iceberg_params
   implicit none
 
   integer                    :: n, istep,outfreq
@@ -283,7 +284,7 @@ subroutine check_blowup(istep, ice, dynamics, tracers, partit, mesh)
     USE MOD_PARTIT
     USE MOD_PARSUP
     USE MOD_MESH
-    use g_config, only: logfile_outfreq, which_ALE, toy_ocean, use_ice
+    use g_config, only: logfile_outfreq, which_ALE, toy_ocean, use_ice, use_icebergs, ib_num
     use o_PARAM
     use o_ARRAYS, only: water_flux, stress_surf, &
                     heat_flux, Kv, Av
@@ -292,6 +293,7 @@ subroutine check_blowup(istep, ice, dynamics, tracers, partit, mesh)
     use g_forcing_arrays
     use diagnostics
     use write_step_info_interface
+    use iceberg_params
     implicit none
 
     type(t_ice)   , intent(in)   , target :: ice
@@ -300,7 +302,7 @@ subroutine check_blowup(istep, ice, dynamics, tracers, partit, mesh)
     type(t_tracer), intent(in)   , target :: tracers
     type(t_mesh)  , intent(in)   , target :: mesh
     !___________________________________________________________________________
-    integer                       :: n, nz, istep, found_blowup_loc=0, found_blowup=0
+    integer                       :: n, nz, istep, found_blowup_loc=0, found_blowup=0, ib
     integer                       :: el, elidx
     !___________________________________________________________________________
     ! pointer on necessary derived types
@@ -388,6 +390,23 @@ subroutine check_blowup(istep, ice, dynamics, tracers, partit, mesh)
           write(*,*)
           write(*,*) 'hnode(:, n)  = ',hnode(ulevels_nod2D(n):nlevels_nod2D(n),n)
           write(*,*)
+          if (use_icebergs) then
+            write(*,*) 'ibhf_n(:, n) = ',ibhf_n(ulevels_nod2D(n):nlevels_nod2D(n),n)
+            write(*,*) 'ibfwb(n) = ',ibfwb(n)
+            write(*,*) 'ibfwl(n) = ',ibfwl(n)
+            write(*,*) 'ibfwe(n) = ',ibfwe(n)
+            write(*,*) 'ibfwbv(n) = ',ibfwbv(n)
+            do ib=1, ib_num
+                if (mesh%elem2d_nodes(1, iceberg_elem(ib)) == n) then
+                    write(*,*) 'ib = ',ib, ', length = ',length_ib(ib), ', height = ', height_ib(ib), ', scaling = ', scaling(ib) 
+                    write(*,*) 'hfb_flux_ib(ib) = ',hfb_flux_ib(ib)
+                    write(*,*) 'hfl_flux_ib(ib,n) = ',hfl_flux_ib(ib,n)
+                    write(*,*) 'hfe_flux_ib(ib) = ',hfe_flux_ib(ib)
+                    write(*,*) 'hfbv_flux_ib(ib,n) = ',hfbv_flux_ib(ib,n)
+                end if
+            end do
+            write(*,*)
+          end if
 !$OMP END CRITICAL
        endif
 
@@ -425,6 +444,23 @@ subroutine check_blowup(istep, ice, dynamics, tracers, partit, mesh)
           write(*,*)
           write(*,*) 'CFL_z(:,n)  = ',CFL_z(:,n)
           write(*,*)
+          if (use_icebergs) then
+            write(*,*) 'ibhf_n(:, n) = ',ibhf_n(ulevels_nod2D(n):nlevels_nod2D(n),n)
+            write(*,*) 'ibfwb(n) = ',ibfwb(n)
+            write(*,*) 'ibfwl(n) = ',ibfwl(n)
+            write(*,*) 'ibfwe(n) = ',ibfwe(n)
+            write(*,*) 'ibfwbv(n) = ',ibfwbv(n)
+            do ib=1, ib_num
+                if (mesh%elem2d_nodes(1, iceberg_elem(ib)) == n) then
+                    write(*,*) 'ib = ',ib, ', length = ',length_ib(ib), ', height = ', height_ib(ib), ', scaling = ', scaling(ib) 
+                    write(*,*) 'hfb_flux_ib(ib) = ',hfb_flux_ib(ib)
+                    write(*,*) 'hfl_flux_ib(ib,n) = ',hfl_flux_ib(ib,n)
+                    write(*,*) 'hfe_flux_ib(ib) = ',hfe_flux_ib(ib)
+                    write(*,*) 'hfbv_flux_ib(ib,n) = ',hfbv_flux_ib(ib,n)
+                end if
+            end do
+            write(*,*)
+          end if
 !$OMP END CRITICAL
        end if ! --> if ( .not. trim(which_ALE)=='linfs' .and. ...
 
@@ -452,6 +488,23 @@ subroutine check_blowup(istep, ice, dynamics, tracers, partit, mesh)
           end if
           write(*,*) 'glon,glat   = ',geo_coord_nod2D(:,n)/rad
           write(*,*)
+          if (use_icebergs) then
+            write(*,*) 'ibhf_n(:, n) = ',ibhf_n(ulevels_nod2D(n):nlevels_nod2D(n),n)
+            write(*,*) 'ibfwb(n) = ',ibfwb(n)
+            write(*,*) 'ibfwl(n) = ',ibfwl(n)
+            write(*,*) 'ibfwe(n) = ',ibfwe(n)
+            write(*,*) 'ibfwbv(n) = ',ibfwbv(n)
+            do ib=1, ib_num
+                if (mesh%elem2d_nodes(1, iceberg_elem(ib)) == n) then
+                    write(*,*) 'ib = ',ib, ', length = ',length_ib(ib), ', height = ', height_ib(ib), ', scaling = ', scaling(ib) 
+                    write(*,*) 'hfb_flux_ib(ib) = ',hfb_flux_ib(ib)
+                    write(*,*) 'hfl_flux_ib(ib,n) = ',hfl_flux_ib(ib,n)
+                    write(*,*) 'hfe_flux_ib(ib) = ',hfe_flux_ib(ib)
+                    write(*,*) 'hfbv_flux_ib(ib,n) = ',hfbv_flux_ib(ib,n)
+                end if
+            end do
+            write(*,*)
+          end if
 !$OMP END CRITICAL
        end if ! --> if ( .not. trim(which_ALE)=='linfs' .and. ...
 
@@ -505,6 +558,23 @@ subroutine check_blowup(istep, ice, dynamics, tracers, partit, mesh)
              write(*,*)
              write(*,*) 'CFL_z(:,n)  = ',CFL_z(:,n)
              write(*,*)
+          if (use_icebergs) then
+            write(*,*) 'ibhf_n(:, n) = ',ibhf_n(ulevels_nod2D(n):nlevels_nod2D(n),n)
+            write(*,*) 'ibfwb(n) = ',ibfwb(n)
+            write(*,*) 'ibfwl(n) = ',ibfwl(n)
+            write(*,*) 'ibfwe(n) = ',ibfwe(n)
+            write(*,*) 'ibfwbv(n) = ',ibfwbv(n)
+            do ib=1, ib_num
+                if (mesh%elem2d_nodes(1, iceberg_elem(ib)) == n) then
+                    write(*,*) 'ib = ',ib, ', length = ',length_ib(ib), ', height = ', height_ib(ib), ', scaling = ', scaling(ib) 
+                    write(*,*) 'hfb_flux_ib(ib) = ',hfb_flux_ib(ib)
+                    write(*,*) 'hfl_flux_ib(ib,n) = ',hfl_flux_ib(ib,n)
+                    write(*,*) 'hfe_flux_ib(ib) = ',hfe_flux_ib(ib)
+                    write(*,*) 'hfbv_flux_ib(ib,n) = ',hfbv_flux_ib(ib,n)
+                end if
+            end do
+            write(*,*)
+          end if
              write(*,*)
 !$OMP END CRITICAL
           endif ! --> if ( (tracers%data(1)%values(nz, n) /= tracers%data(1)%values(nz, n)) .or. & ...
@@ -560,6 +630,23 @@ subroutine check_blowup(istep, ice, dynamics, tracers, partit, mesh)
              write(*,*)
              write(*,*) 'glon,glat   = ',geo_coord_nod2D(:,n)/rad
              write(*,*)
+          if (use_icebergs) then
+            write(*,*) 'ibhf_n(:, n) = ',ibhf_n(ulevels_nod2D(n):nlevels_nod2D(n),n)
+            write(*,*) 'ibfwb(n) = ',ibfwb(n)
+            write(*,*) 'ibfwl(n) = ',ibfwl(n)
+            write(*,*) 'ibfwe(n) = ',ibfwe(n)
+            write(*,*) 'ibfwbv(n) = ',ibfwbv(n)
+            do ib=1, ib_num
+                if (mesh%elem2d_nodes(1, iceberg_elem(ib)) == n) then
+                    write(*,*) 'ib = ',ib, ', length = ',length_ib(ib), ', height = ', height_ib(ib), ', scaling = ', scaling(ib) 
+                    write(*,*) 'hfb_flux_ib(ib) = ',hfb_flux_ib(ib)
+                    write(*,*) 'hfl_flux_ib(ib,n) = ',hfl_flux_ib(ib,n)
+                    write(*,*) 'hfe_flux_ib(ib) = ',hfe_flux_ib(ib)
+                    write(*,*) 'hfbv_flux_ib(ib,n) = ',hfbv_flux_ib(ib,n)
+                end if
+            end do
+            write(*,*)
+          end if
 !$OMP END CRITICAL
           endif ! --> if ( (tracers%data(2)%values(nz, n) /= tracers%data(2)%values(nz, n)) .or.  & ...
        end do ! --> do nz=1,nlevels_nod2D(n)-1