Merge pull request #590 from FESOM/feat/cav_plus_icb

Feat/cav plus icb

src/MOD_ICE.F90

@@ -762,15 +762,13 @@ subroutine ice_init(ice, partit, mesh)
     !array of 2D boundary conditions is used in ice_maEVP
 
     ! LA 2023-05-24 initiate bc_index_nod2D also for whichEVP==0
-    !if (ice%whichEVP > 0) then
         allocate(mesh%bc_index_nod2D(myDim_nod2D+eDim_nod2D))
         mesh%bc_index_nod2D=1._WP
         do n=1, myDim_edge2D
             ed=mesh%edges(:, n)
             if (myList_edge2D(n) <= mesh%edge2D_in) cycle
             mesh%bc_index_nod2D(ed)=0._WP
         end do
-    !end if
 
 end subroutine ice_init  
 !

src/icb_dyn.F90

@@ -664,12 +664,10 @@ subroutine iceberg_average_andkeel(mesh, partit, dynamics, uo_dz,vo_dz, uo_keel,
     !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)
+    ! 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
 
-#ifdef use_cavity
-    ! if cavity node ..
-    if( cavity_flag_nod2d(elem2D_nodes(m,iceberg_elem))==1 .AND. abs(depth_ib)<abs(lev_up) ) then
+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)
@@ -689,10 +687,8 @@ subroutine iceberg_average_andkeel(mesh, partit, dynamics, uo_dz,vo_dz, uo_keel,
 
     !****************************************************************
     ! LA 23.11.21 case if depth_ib<lev_up
-    else if( abs(lev_low)>=abs(depth_ib) ) then !.AND. (abs(lev_up)<=abs(depth_ib)) ) then
-#else
+else !# comp cav flag
     if( abs(lev_low)>=abs(depth_ib) ) then !.AND. (abs(lev_up)<=abs(depth_ib)) ) then
-#endif
       if( abs(lev_up)<abs(depth_ib) ) then
         dz = abs ( lev_up - depth_ib )
       else
@@ -795,7 +791,7 @@ subroutine iceberg_average_andkeel(mesh, partit, dynamics, uo_dz,vo_dz, uo_keel,
         S_keel(m)=Sclim_ib(k,n2)
       end if
     end if
-
+end if !cavity
    end do innerloop
 
    ! divide by depth over which was integrated

src/icb_elem.F90

@@ -376,15 +376,14 @@ end subroutine FEM_3eval
 subroutine iceberg_elem4all(mesh, partit, elem, lon_deg, lat_deg) 
  USE MOD_MESH
  use MOD_PARTIT		!for myDim_nod2D, myList_elem2D
-!#ifdef use_cavity
 ! use iceberg_params, only: reject_elem
-!#endif
 
  implicit none
 
  integer, intent(INOUT) :: elem
  real, intent(IN) :: lon_deg, lat_deg
  logical :: i_have_element
+ logical                        :: reject_tmp !LA for debugging
 
 type(t_mesh), intent(in) , target :: mesh
 type(t_partit), intent(inout), target :: partit
@@ -398,17 +397,19 @@ 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.
-#ifdef use_cavity
-   if( reject_elem(mesh, elem) ) then
-    elem=0 !reject element
-    i_have_element=.false.
-    !write(*,*) 'elem4all: iceberg found in shelf region: elem = 0'
-   else 
-    elem=myList_elem2D(elem) !global now
-   end if 
-#else
+   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) )
+      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'
+      else 
+       elem=myList_elem2D(elem) !global now
+      end if 
+   else
    elem=myList_elem2D(elem) !global now
-#endif 
+endif 
   end if
   call com_integer(partit, i_have_element,elem) !max 1 PE sends element here; 
 end subroutine iceberg_elem4all
@@ -419,9 +420,7 @@ end subroutine iceberg_elem4all
 
 subroutine find_new_iceberg_elem(mesh, partit, old_iceberg_elem, pt, left_mype)
   use o_param
-!#ifdef use_cavity
 !  use iceberg_params, only: reject_elem
-!#endif
 
   implicit none
 
@@ -431,6 +430,7 @@ subroutine find_new_iceberg_elem(mesh, partit, old_iceberg_elem, pt, left_mype)
 
   INTEGER :: m, n2, idx_elem_containing_n2, elem_containing_n2, ibelem_tmp
   REAL, DIMENSION(3) :: werte2D
+ logical                        :: reject_tmp !LA for debugging
 
 type(t_mesh), intent(in) , target :: mesh
 type(t_partit), intent(inout), target :: partit
@@ -461,13 +461,15 @@ subroutine find_new_iceberg_elem(mesh, partit, old_iceberg_elem, pt, left_mype)
   if (ALL(werte2D <= 1.+ 1.0e-07) .AND. ALL(werte2D >= 0.0- 1.0e-07) ) then
    old_iceberg_elem=elem_containing_n2
 
-#ifdef use_cavity
-   if( reject_elem(mesh, old_iceberg_elem) ) then
-      left_mype=1.0
-      !write(*,*) 'iceberg found in shelf region: left_mype = 1'
-      old_iceberg_elem=ibelem_tmp
-   end if
-#endif
+   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) )
+      if(reject_tmp) then
+         left_mype=1.0
+         !write(*,*) 'iceberg found in shelf region: left_mype = 1'
+         old_iceberg_elem=ibelem_tmp
+      end if
+   endif
 
    RETURN 
   end if

src/icb_modules.F90

@@ -1,4 +1,7 @@
 module iceberg_params
+use MOD_PARTIT
+USE MOD_MESH
+use g_config, only: use_cavity, use_cavityonelem
 implicit none
 save
   !integer,parameter :: ib_num               ! realistic dataset comprising 6912 icebergs
@@ -110,36 +113,56 @@ module iceberg_params
   real,dimension(:,:), allocatable:: buoy_props
   integer:: save_count_buoys
   real:: prev_sec_in_year
+
 !****************************************************************************************************************************
 !****************************************************************************************************************************
-#ifdef use_cavity
  contains
  ! true if all nodes of the element are either "real" model boundary nodes or shelf nodes
- logical function reject_elem(mesh, elem)
- USE MOD_MESH
- implicit none
- integer, intent(in) :: elem
+ logical function reject_elem(mesh, partit, elem)
+
+integer, intent(in) :: elem
 type(t_mesh), intent(in) , target :: mesh
+type(t_partit), intent(inout), target :: partit
 #include "associate_mesh_def.h"
+#include "associate_part_def.h"
 #include "associate_mesh_ass.h"
+#include "associate_part_ass.h"
 
+write(*,*) "LA DEBUG 3"
+if (use_cavity) then
 ! kh 09.08.21 change index_nod2d -> bc_index_nod2d?
- reject_elem = all( (cavity_flag_nod2d(elem2D_nodes(:,elem))==1) .OR. (index_nod2d(elem2D_nodes(:,elem))==1) )
+ if (.not. use_cavityonelem) then
+   write(*,*) "LA DEBUG: cavity_depth = ", mesh%cavity_depth
+   write(*,*) "LA DEBUG: cavity_flag = ", mesh%cavity_depth(mesh%elem2D_nodes(:,elem))
+   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
+ reject_elem = all( (mesh%bc_index_nod2D(mesh%elem2D_nodes(:,elem))==0.0) )
+endif
  end function reject_elem
 
  ! gives number of "coastal" nodes in cavity setup, i.e. number of nodes that are
  ! either "real" model boundary nodes or shelf nodes
- integer function coastal_nodes(mesh, elem)
- USE MOD_MESH
- implicit none
- integer, intent(in) :: elem
+ integer function coastal_nodes(mesh, partit, elem)
+
+integer, intent(in) :: elem
 type(t_mesh), intent(in) , target :: mesh
+type(t_partit), intent(inout), target :: partit
+#include "associate_mesh_def.h"
 #include "associate_part_def.h"
+#include "associate_mesh_ass.h"
 #include "associate_part_ass.h"
 
+if (use_cavity) then
 ! kh 09.08.21 change index_nod2d -> bc_index_nod2d?
- coastal_nodes = count( (cavity_flag_nod2d(elem2D_nodes(:,elem))==1) .OR. (index_nod2d(elem2D_nodes(:,elem))==1) )
+ if (.not. use_cavityonelem) then
+   coastal_nodes = count( (mesh%cavity_depth(mesh%elem2D_nodes(:,elem))/=0.0) .OR. (mesh%bc_index_nod2D(mesh%elem2D_nodes(:,elem))==0.0) )
+ !else
+ end if
+else 
+ coastal_nodes = count( (mesh%bc_index_nod2D(mesh%elem2D_nodes(:,elem))==0.0) )
+endif
  end function coastal_nodes
-#endif
 
 end module iceberg_params

src/icb_step.F90

@@ -324,7 +324,7 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib,length_ib,wi
  											!=
  real, dimension(2)             :: coords_tmp
 ! integer, pointer  :: mype
-
+ logical                        :: reject_tmp !LA for debugging
 type(t_ice),  intent(inout), target :: ice
 type(t_mesh), intent(in) , target :: mesh
 type(t_partit), intent(inout), target :: partit
@@ -375,21 +375,25 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib,length_ib,wi
   coords_tmp = [lon_deg, lat_deg]
   call point_in_triangle(mesh, partit, iceberg_elem, coords_tmp)
   !call point_in_triangle(mesh, iceberg_elem, (/lon_deg, lat_deg/))
-  i_have_element= (iceberg_elem .ne. 0) !up to 3 PEs possible
+  i_have_element= (iceberg_elem .ne. 0) !up to 3 PEs .true.
 
   if(i_have_element) then
    i_have_element= mesh%elem2D_nodes(1,iceberg_elem) <= partit%myDim_nod2D !1 PE still .true.
-#ifdef use_cavity
-   if(reject_elem(mesh, partit, iceberg_elem)) then
-    iceberg_elem=0 !reject element
-    i_have_element=.false.
-   else 
-    iceberg_elem=partit%myList_elem2D(iceberg_elem) !global now
-   end if
-#else
-
-   iceberg_elem=partit%myList_elem2D(iceberg_elem) !global now
-#endif 
+   if (use_cavity) then
+      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_elem(mesh, partit, iceberg_elem)) then
+      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))==1)
+       iceberg_elem=0 !reject element
+       i_have_element=.false.
+      else 
+       iceberg_elem=partit%myList_elem2D(iceberg_elem) !global now
+      end if
+   else
+      iceberg_elem=partit%myList_elem2D(iceberg_elem) !global now
+   endif 
   end if
   call com_integer(partit, i_have_element,iceberg_elem)
 
@@ -983,9 +987,7 @@ end subroutine depth_bathy
 !****************************************************************************************************************************
 
 subroutine parallel2coast(mesh, partit,u, v, lon,lat, elem)
-!#ifdef use_cavity
-! use iceberg_params, only: coastal_nodes
-!#endif
+ use iceberg_params, only: coastal_nodes
  implicit none
 
  real, intent(inout) 	:: u, v 	!velocity
@@ -1004,12 +1006,12 @@ subroutine parallel2coast(mesh, partit,u, v, lon,lat, elem)
 #include "associate_part_ass.h"
 #include "associate_mesh_ass.h"
 
-#ifdef use_cavity
-  SELECT CASE ( coastal_nodes(mesh, elem) ) !num of "coastal" points
-#else
+!if (use_cavity) then
+!  SELECT CASE ( coastal_nodes(mesh, partit, elem) ) !num of "coastal" points
+!else
   SELECT CASE ( sum( mesh%bc_index_nod2D(mesh%elem2D_nodes(:,elem)) ) ) !num of coastal points
   !SELECT CASE ( sum( bc_index_nod2D(elem2D_nodes(:,elem)) ) ) !num of coastal points
-#endif
+!endif
    CASE (0) !...coastal points: do nothing
     return
 
@@ -1020,14 +1022,18 @@ subroutine parallel2coast(mesh, partit,u, v, lon,lat, elem)
     do m = 1, 3
       node = mesh%elem2D_nodes(m,elem)
       !write(*,*) 'index ', m, ':', index_nod2D(node)
-#ifdef use_cavity
-      if( mesh%bc_index_nod2D(node)==1 .OR. cavity_flag_nod2d(node)==1 ) then
-#else
-      if( mesh%bc_index_nod2D(node)==1 ) then
-#endif
-       n(i) = node
-       exit
-      end if 
+      if (use_cavity) then
+            !if( mesh%bc_index_nod2D(node)==1 .OR. mesh%cavity_flag_n(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 
+      endif
     end do 
 
    !write(*,*) 'one coastal node ', n(1)
@@ -1039,7 +1045,7 @@ subroutine parallel2coast(mesh, partit,u, v, lon,lat, elem)
    ! do m = 1, nghbr_nod2D(n(1))%nmb
    !   node = nghbr_nod2D(n(1))%addresses(m) 
 !#ifdef use_cavity
-   !   if ( (node /= n(1)) .and. ( (bc_index_nod2D(node)==1) .OR. (cavity_flag_nod2d(node)==1) ) ) then   
+   !   if ( (node /= n(1)) .and. ( (bc_index_nod2D(node)==1) .OR. (mesh%cavity_flag_n(node)==1) ) ) then   
 !#else
    !   if ( (node /= n(1)) .and. (bc_index_nod2D(node)==1)) then
 !#endif
@@ -1075,14 +1081,18 @@ subroutine parallel2coast(mesh, partit,u, v, lon,lat, elem)
     velocity = [ u, v ]
     do m = 1, 3
       node = mesh%elem2D_nodes(m,elem) 
-#ifdef use_cavity
-      if( (mesh%bc_index_nod2D(node)==1) .OR. (cavity_flag_nod2d(node)==1)) then
-#else
-      if( mesh%bc_index_nod2D(node)==1 ) then
-#endif
+if (use_cavity) then
+      !if( (mesh%bc_index_nod2D(node)==1) .OR. (mesh%cavity_flag_n(node)==1)) then
+      if( (mesh%bc_index_nod2D(node)==0.0) .OR. (mesh%cavity_depth(node)/=0.0) ) then
+       n(i) = node
+       i = i+1
+      end if
+else
+      if( mesh%bc_index_nod2D(node)==0.0 ) then
        n(i) = node
        i = i+1
       end if
+endif
     end do   
     call projection(mesh,partit,  velocity, n(1), n(2))
 

src/icb_thermo.F90

@@ -470,7 +470,7 @@ end subroutine potit_ib
 
 ! if the underlying FESOM is run without cavities, the following routines might be 
 ! missing, so put them here:
-#ifndef use_cavity
+!if (.not. use_cavity) then
 !
 !-------------------------------------------------------------------------------------
 !
@@ -572,7 +572,7 @@ end subroutine potit_ib
 !
 !----------------------------------------------------------------------------------------
 !
-#endif
+!endif
 
 
 ! LA from oce_dens_press for iceberg coupling