Merge branch 'main' into workbench_mainf26_icepack

src/cpl_driver.F90

@@ -16,6 +16,7 @@ module cpl_driver
   use g_config, only : dt, use_icebergs, lwiso
   use o_param,  only : rad
   USE MOD_PARTIT
+  use mpi
   implicit none
   save   
   !
@@ -56,7 +57,7 @@ module cpl_driver
   integer                    :: localRank      ! local MPI rank
   integer                    :: localSize      ! local MPI size
   integer                    :: localComm      ! local MPI size
-  logical                    :: commRank       ! true for ranks doing OASIS communication
+  integer                    :: commRank
   integer                    :: comp_id        ! id returned by oasis_init_comp
 
   logical, save              :: oasis_was_initialized
@@ -773,11 +774,11 @@ subroutine cpl_oasis3mct_define_unstr(partit, mesh)
 !------------------------------------------------------------------
 
    call oasis_enddef(ierror)
-   if (commRank) print *, 'fesom oasis_enddef: COMPLETED'
+   if (ierror .eq. oasis_ok) print *, 'fesom oasis_enddef: COMPLETED'
 #ifndef __oifs
-   if (commRank) print *, 'FESOM: calling exchange_roots'
+   if (ierror .eq. oasis_ok) print *, 'FESOM: calling exchange_roots'
    call exchange_roots(source_root, target_root, 1, partit%MPI_COMM_FESOM, MPI_COMM_WORLD)
-   if (commRank) print *, 'FESOM source/target roots: ', source_root, target_root
+   if (ierror .eq. oasis_ok) print *, 'FESOM source/target roots: ', source_root, target_root
 #endif
 
    ! WAS VOM FOLGENDEN BRAUCHE ICH NOCH ??? 

src/fesom_module.F90

@@ -36,6 +36,7 @@ module fesom_main_storage_module
   use age_tracer_init_interface
   use iceberg_params
   use iceberg_step
+  use iceberg_ocean_coupling
   ! Define icepack module
 
 #if defined (__icepack)
@@ -272,7 +273,7 @@ subroutine fesom_init(fesom_total_nsteps)
         ! --------------
         ! LA icebergs: 2023-05-17 
         if (use_icebergs) then
-            call allocate_icb(f%partit)
+            call allocate_icb(f%partit, f%mesh)
         endif
         ! --------------
 
@@ -570,16 +571,6 @@ subroutine fesom_runloop(current_nsteps)
             if (flag_debug .and. f%mype==0)  print *, achar(27)//'[34m'//' --> call ice_timestep(n)'//achar(27)//'[0m'
             if (f%ice%ice_update) call ice_timestep(n, f%ice, f%partit, f%mesh)  
 
-
-            ! LA commented for debugging
-            ! --------------
-            ! LA icebergs: 2023-05-17 
-            if (use_icebergs .and. mod(n, steps_per_ib_step)==0.0) then
-                call icb2fesom(f%mesh, f%partit, f%ice)
-            end if
-            ! --------------
-
-
             !___compute fluxes to the ocean: heat, freshwater, momentum_________
             if (flag_debug .and. f%mype==0)  print *, achar(27)//'[34m'//' --> call oce_fluxes_mom...'//achar(27)//'[0m'
             call oce_fluxes_mom(f%ice, f%dynamics, f%partit, f%mesh) ! momentum only

src/gen_modules_config.F90

@@ -122,6 +122,10 @@ module g_config
   logical                       :: turn_off_hf=.false.
   logical                       :: turn_off_fw=.false.
   logical                       :: use_icesheet_coupling=.false.  
+  logical                       :: lbalance_fw=.true.
+  integer                       :: cell_saturation=2 ! 0=no cell saturation, 1=one additional iceberg allowed, 2=no daddtional iceberg allowed
+  logical                       :: lmin_latent_hf=.true.
+  logical                       :: lverbose_icb=.false.  
   integer                       :: ib_num=0
   integer                       :: steps_per_ib_step=8
 
@@ -132,7 +136,8 @@ module g_config
   integer                       :: ib_async_mode=0
   integer                       :: thread_support_level_required=3 ! 2 = MPI_THREAD_SERIALIZED, 3 = MPI_THREAD_MULTIPLE
 
-  namelist /icebergs/ use_icebergs, turn_off_hf, turn_off_fw, use_icesheet_coupling, ib_num, steps_per_ib_step, ib_async_mode, thread_support_level_required
+  namelist /icebergs/   use_icebergs, turn_off_hf, turn_off_fw, use_icesheet_coupling, lbalance_fw, cell_saturation, lmin_latent_hf, &
+                        ib_num, steps_per_ib_step, ib_async_mode, thread_support_level_required, lverbose_icb
 
 !wiso-code!!!
   logical                       :: lwiso  =.false.  ! enable isotope?

src/icb_allocate.F90

@@ -1,20 +1,30 @@
-subroutine allocate_icb(partit)
+subroutine allocate_icb(partit, mesh)
   use iceberg_params
   use g_config
+  use g_comm
+  use g_comm_auto
+  use o_param
   use MOD_PARTIT
+  use MOD_MESH
 
   integer       :: n2
 type(t_partit), intent(inout), target :: partit
+type(t_mesh), intent(in), target :: mesh
 #include "associate_part_def.h"
+#include "associate_mesh_def.h"
 #include "associate_part_ass.h"
+#include "associate_mesh_ass.h"
   n2=myDim_nod2D+eDim_nod2D
+  icb_outfreq = step_per_day / steps_per_ib_step
 
   allocate(ibhf(n2), ibfwb(n2), ibfwl(n2), ibfwe(n2), ibfwbv(n2))
-  ibhf=0
-  ibfwb=0
-  ibfwl=0
-  ibfwe=0
-  ibfwbv=0
+  ibhf      = 0.0
+  ibfwb     = 0.0
+  ibfwl     = 0.0
+  ibfwe     = 0.0
+  ibfwbv    = 0.0
+  allocate(ibhf_n(mesh%nl, n2))
+  ibhf_n    = 0.0_WP
 
   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
@@ -83,16 +93,29 @@ subroutine allocate_icb(partit)
   allocate(fwl_flux_ib(ib_num))
   allocate(fwb_flux_ib(ib_num))
   allocate(fwbv_flux_ib(ib_num))
-  allocate(heat_flux_ib(ib_num))
-  allocate(lheat_flux_ib(ib_num))
+  allocate(hfe_flux_ib(ib_num))
+  allocate(hfl_flux_ib(ib_num,mesh%nl))
+  allocate(hfb_flux_ib(ib_num))
+  allocate(hfbv_flux_ib(ib_num,mesh%nl))
+  allocate(lhfb_flux_ib(ib_num))
   fwe_flux_ib = 0.0
   fwl_flux_ib = 0.0
   fwb_flux_ib = 0.0
   fwbv_flux_ib = 0.0
-  heat_flux_ib = 0.0
-  lheat_flux_ib = 0.0
+  hfe_flux_ib = 0.0
+  hfl_flux_ib = 0.0
+  hfb_flux_ib = 0.0
+  hfbv_flux_ib = 0.0
+  lhfb_flux_ib = 0.0
   allocate(arr_block(15*ib_num))
   allocate(elem_block(ib_num))
+  allocate(pe_block(ib_num))
+
+  allocate(elem_area_glob(elem2D))
+  elem_area_glob=0.0
+  call gather_elem(elem_area(1:myDim_elem2D), elem_area_glob, partit)
+  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))
   buoy_props = 0.0

src/icb_coupling.F90

@@ -1,19 +1,34 @@
+module iceberg_ocean_coupling
+ USE MOD_MESH
+ use MOD_PARTIT
+ use MOD_ICE
+ USE MOD_DYN
+ use iceberg_params
+
+ public :: reset_ib_fluxes
+ public :: prepare_icb2fesom
+ public :: icb2fesom
+
+ contains
+
+
 subroutine reset_ib_fluxes()
     use o_param
 
 ! kh 18.03.21 not really used here
     use g_config
     use iceberg_params
 
-    ibfwbv = 0
-    ibfwb = 0
-    ibfwl = 0
-    ibfwe = 0
-    ibhf = 0
+    ibfwbv  = 0.0
+    ibfwb   = 0.0
+    ibfwl   = 0.0
+    ibfwe   = 0.0
+    ibhf    = 0.0
+    ibhf_n  = 0.0_WP
 end subroutine
 
 
-subroutine prepare_icb2fesom(mesh, partit, ib,i_have_element,localelement,depth_ib)
+subroutine prepare_icb2fesom(mesh, partit, ib,i_have_element,localelement,depth_ib,height_ib_single)
     !transmits the relevant fields from the iceberg to the ocean model
     !Lars Ackermann, 17.03.2020
 
@@ -23,46 +38,116 @@ subroutine prepare_icb2fesom(mesh, partit, ib,i_have_element,localelement,depth_
     use MOD_PARTIT
     use iceberg_params
 
+    implicit none
+
     logical                 :: i_have_element
-    real                    :: depth_ib
+    real, intent(in)        :: depth_ib, height_ib_single
+    real                    :: lev_low, lev_up
     integer                 :: localelement
     integer                 :: iceberg_node  
     integer, dimension(3)   :: ib_nods_in_ib_elem
     integer                 :: num_ib_nods_in_ib_elem
-    real                    :: tot_area_nods_in_ib_elem
-    integer                 :: i, ib
+    real                    :: dz
+    real, dimension(:), allocatable      :: tot_area_nods_in_ib_elem
+    integer                 :: i, j, k, ib
+    integer, dimension(3)   :: idx_d
 type(t_mesh), intent(in) , target :: mesh
 type(t_partit), intent(inout), target :: partit
 #include "associate_part_def.h"
 #include "associate_mesh_def.h"
 #include "associate_part_ass.h"
 #include "associate_mesh_ass.h"
 
-    if(i_have_element) then 
+    if(i_have_element) then
+        dz = 0.0
+        allocate(tot_area_nods_in_ib_elem(mesh%nl))
+
         num_ib_nods_in_ib_elem=0
-        tot_area_nods_in_ib_elem=0
+        tot_area_nods_in_ib_elem=0.0
+        idx_d = 0
 
         do i=1,3
             iceberg_node=elem2d_nodes(i,localelement)
 
+            ! 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
+                lev_up  = mesh%zbar_3d_n(k, iceberg_node)
+
+                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)
+                end if
+
+                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
+                    exit
+                else
+                    cycle
+                end if
+            end do
+
             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(1,iceberg_node)
+                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
+                ib_nods_in_ib_elem(i)           = 0
             end if
         end do
-        
+
         do i=1, 3
             iceberg_node=ib_nods_in_ib_elem(i)
+            if (use_cavity .and. ulevels_nod2d(iceberg_node) > 1) cycle
 
             if (iceberg_node>0) then
-                ibfwbv(iceberg_node) = ibfwbv(iceberg_node) - fwbv_flux_ib(ib) / tot_area_nods_in_ib_elem
-                ibfwb(iceberg_node) = ibfwb(iceberg_node) - fwb_flux_ib(ib) / tot_area_nods_in_ib_elem
-                ibfwl(iceberg_node) = ibfwl(iceberg_node) - fwl_flux_ib(ib) / tot_area_nods_in_ib_elem
-                ibfwe(iceberg_node) = ibfwe(iceberg_node) - fwe_flux_ib(ib) / tot_area_nods_in_ib_elem
-                ibhf(iceberg_node) = ibhf(iceberg_node) - heat_flux_ib(ib) / tot_area_nods_in_ib_elem
+                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)
+                !ibhf(iceberg_node) = ibhf(iceberg_node) - hfb_flux_ib(ib) / tot_area_nods_in_ib_elem(1)
+
+                do j=1,idx_d(i)
+                    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)
+                    end if
+                    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))
+                    end if              
+
+                    if( 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)) * (dz / abs(depth_ib)) & 
+                                                    + hfe_flux_ib(ib) * (dz / abs(height_ib_single))) &
+                                                    / tot_area_nods_in_ib_elem(j)
+                    end if
+                end do
+
+                if( idx_d(i) > 1 ) then
+                    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))
+                    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
+                    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( 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)) & 
+                            / tot_area_nods_in_ib_elem(1)
+                end if
             end if
         end do
     end if
@@ -90,18 +175,25 @@ subroutine icb2fesom(mesh, partit, ice)
 #include "associate_mesh_def.h"
 #include "associate_part_ass.h"
 #include "associate_mesh_ass.h"
-
-    fresh_wa_flux => ice%flx_fw(:)
-    net_heat_flux => ice%flx_h(:)
 
+  if (use_cavity) then
+!$OMP PARALLEL DO
     do n=1, myDim_nod2d+eDim_nod2D
-        if (.not.turn_off_hf) then
-                net_heat_flux(n)   = net_heat_flux(n) + ibhf(n) * steps_per_ib_step 
+        if (ulevels_nod2d(n) > 1) cycle
+        if (.not.turn_off_fw) then
+                water_flux(n)  = water_flux(n) - (ibfwb(n)+ibfwl(n)+ibfwe(n)+ibfwbv(n)) !* steps_per_ib_step
         end if
+    end do
+!$OMP END PARALLEL DO
+  else
+!$OMP PARALLEL DO
+    do n=1, myDim_nod2d+eDim_nod2D
         if (.not.turn_off_fw) then
-                fresh_wa_flux(n)   = fresh_wa_flux(n) + (ibfwb(n)+ibfwl(n)+ibfwe(n)+ibfwbv(n)) * steps_per_ib_step
+                water_flux(n)  = water_flux(n) - (ibfwb(n)+ibfwl(n)+ibfwe(n)+ibfwbv(n)) !* steps_per_ib_step
         end if
     end do
+!$OMP END PARALLEL DO
+  end if
 !---wiso-code-begin
     if(lwiso) then
       do n=1, myDim_nod2D+eDim_nod2D
@@ -112,3 +204,4 @@ subroutine icb2fesom(mesh, partit, ice)
     end if
 !---wiso-code-end
 end subroutine icb2fesom
+end module