From 361b9b9d3c0f1f4e773ea41f2c98796e5a50af34 Mon Sep 17 00:00:00 2001 From: ackerlar Date: Fri, 20 Sep 2024 14:45:15 +0200 Subject: [PATCH 01/14] 3D heat flux, fix cell saturation --- src/fesom_module.F90 | 12 +- src/gen_modules_config.F90 | 5 +- src/icb_allocate.F90 | 43 +++-- src/icb_coupling.F90 | 105 ++++++++--- src/icb_coupling.F90.rej | 22 +++ src/icb_dyn.F90 | 30 ++-- src/icb_elem.F90 | 39 ++--- src/icb_elem.F90.rej | 59 +++++++ src/icb_modules.F90 | 24 +-- src/icb_modules.F90.rej | 72 ++++++++ src/icb_step.F90 | 327 +++++++++++++++++++---------------- src/icb_step.F90.rej | 121 +++++++++++++ src/icb_thermo.F90 | 179 ++++++------------- src/icb_thermo.F90.rej | 324 ++++++++++++++++++++++++++++++++++ src/ice_oce_coupling.F90 | 13 ++ src/ice_oce_coupling.F90.rej | 23 +++ src/io_meandata.F90.rej | 11 ++ src/oce_ale_tracer.F90 | 11 ++ src/oce_ale_tracer.F90.rej | 9 + 19 files changed, 1055 insertions(+), 374 deletions(-) create mode 100644 src/icb_coupling.F90.rej create mode 100644 src/icb_elem.F90.rej create mode 100644 src/icb_modules.F90.rej create mode 100644 src/icb_step.F90.rej create mode 100644 src/icb_thermo.F90.rej create mode 100644 src/ice_oce_coupling.F90.rej create mode 100644 src/io_meandata.F90.rej create mode 100644 src/oce_ale_tracer.F90.rej diff --git a/src/fesom_module.F90 b/src/fesom_module.F90 index 3302f61f4..1fb782a0f 100755 --- a/src/fesom_module.F90 +++ b/src/fesom_module.F90 @@ -272,7 +272,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 +570,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 diff --git a/src/gen_modules_config.F90 b/src/gen_modules_config.F90 index a12ff8c9a..f4766764e 100755 --- a/src/gen_modules_config.F90 +++ b/src/gen_modules_config.F90 @@ -122,6 +122,9 @@ module g_config logical :: turn_off_hf=.false. logical :: turn_off_fw=.false. logical :: use_icesheet_coupling=.false. + logical :: lbalance_fw=.true. + logical :: lcell_saturation=.true. + logical :: lmin_latent_hf=.true. integer :: ib_num=0 integer :: steps_per_ib_step=8 @@ -132,7 +135,7 @@ 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, lcell_saturation, lmin_latent_hf, ib_num, steps_per_ib_step, ib_async_mode, thread_support_level_required !wiso-code!!! logical :: lwiso =.false. ! enable isotope? diff --git a/src/icb_allocate.F90 b/src/icb_allocate.F90 index 169f73126..8d235a415 100644 --- a/src/icb_allocate.F90 +++ b/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)) + allocate(hfb_flux_ib(ib_num)) + allocate(hfbv_flux_ib(ib_num)) + 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 diff --git a/src/icb_coupling.F90 b/src/icb_coupling.F90 index e4331ecd7..56004644a 100644 --- a/src/icb_coupling.F90 +++ b/src/icb_coupling.F90 @@ -5,11 +5,12 @@ subroutine reset_ib_fluxes() 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 @@ -23,14 +24,19 @@ 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 + 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" @@ -39,30 +45,83 @@ subroutine prepare_icb2fesom(mesh, partit, ib,i_have_element,localelement,depth_ #include "associate_mesh_ass.h" 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 (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) ice%flx_fw(:) - net_heat_flux => ice%flx_h(:) 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 - end if 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 !---wiso-code-begin diff --git a/src/icb_coupling.F90.rej b/src/icb_coupling.F90.rej new file mode 100644 index 000000000..fd3215347 --- /dev/null +++ b/src/icb_coupling.F90.rej @@ -0,0 +1,22 @@ +diff a/src/icb_coupling.F90 b/src/icb_coupling.F90 (rejected hunks) +@@ -90,16 +149,13 @@ type(t_partit), intent(inout), target :: partit + #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(:) + + do n=1, myDim_nod2d+eDim_nod2D +- if (.not.turn_off_hf) then +- net_heat_flux(n) = net_heat_flux(n) + ibhf(n) +- end if ++ !if (.not.turn_off_hf) then ++ ! heat_flux(n) = heat_flux(n) - ibhf(n) !* steps_per_ib_step ++ !end if + if (.not.turn_off_fw) then +- fresh_wa_flux(n) = fresh_wa_flux(n) + (ibfwb(n)+ibfwl(n)+ibfwe(n)+ibfwbv(n)) ++ water_flux(n) = water_flux(n) - (ibfwb(n)+ibfwl(n)+ibfwe(n)+ibfwbv(n)) !* steps_per_ib_step + end if + end do + !---wiso-code-begin diff --git a/src/icb_dyn.F90 b/src/icb_dyn.F90 index 78e3ad0d7..bc5f785e2 100644 --- a/src/icb_dyn.F90 +++ b/src/icb_dyn.F90 @@ -135,11 +135,12 @@ subroutine iceberg_dyn(mesh, partit, ice, dynamics, ib, new_u_ib, new_v_ib, u_ib if(l_melt) then call FEM_eval(mesh, partit, sst_ib,sss_ib,lon,lat,Tsurf_ib,Ssurf_ib,iceberg_elem) - call iceberg_meltrates( M_b, M_v, M_e, M_bv, & + + call iceberg_meltrates(partit, M_b, M_v, M_e, M_bv, & u_ib,v_ib, uo_ib,vo_ib, ua_ib,va_ib, & sst_ib, length_ib, conci_ib, & - uo_skin_ib, vo_skin_ib, T_keel_ib, S_keel_ib, depth_ib, & - T_ave_ib, S_ave_ib, ib) + uo_skin_ib, vo_skin_ib, T_keel_ib, S_keel_ib, depth_ib, height_ib, & + T_ave_ib, S_ave_ib, ib, rho_icb) call iceberg_newdimensions(partit, ib, depth_ib,height_ib,length_ib,width_ib,M_b,M_v,M_e,M_bv, & rho_h2o, rho_icb, file3) @@ -644,19 +645,26 @@ subroutine iceberg_average_andkeel(mesh, partit, dynamics, uo_dz,vo_dz, uo_keel, ! below n2.. !innerloop: do k=1, nl+1 innerloop: do k=1, nlevels_nod2D(n2) - - 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 + lev_up = mesh%zbar_3d_n(k, n2) if( k==nlevels_nod2D(n2) ) then lev_low = mesh%zbar_n_bot(n2) else - lev_low = mesh%Z_3d_n_ib(k, n2) + lev_low = mesh%zbar_3d_n(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 diff --git a/src/icb_elem.F90 b/src/icb_elem.F90 index b0a1080cc..077630a4f 100644 --- a/src/icb_elem.F90 +++ b/src/icb_elem.F90 @@ -376,14 +376,12 @@ 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 -! use iceberg_params, only: reject_elem 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 @@ -397,30 +395,21 @@ 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) ) - 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 + if( (use_cavity) .and. (reject_elem(mesh, partit, elem) )) then + elem=0 !reject element + i_have_element=.false. + else + elem=myList_elem2D(elem) !global now + end if end if call com_integer(partit, i_have_element,elem) !max 1 PE sends element here; end subroutine iceberg_elem4all !*************************************************************************************************************************** - !*************************************************************************************************************************** subroutine find_new_iceberg_elem(mesh, partit, old_iceberg_elem, pt, left_mype) use o_param -! use iceberg_params, only: reject_elem implicit none @@ -430,7 +419,6 @@ 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,16 +449,11 @@ 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 - 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 - + if( (use_cavity) .and. (reject_elem(mesh, partit, old_iceberg_elem) )) then + left_mype=1.0 + !write(*,*) 'iceberg found in shelf region: left_mype = 1' + old_iceberg_elem=ibelem_tmp + end if RETURN end if end do diff --git a/src/icb_elem.F90.rej b/src/icb_elem.F90.rej new file mode 100644 index 000000000..a34f6a9f6 --- /dev/null +++ b/src/icb_elem.F90.rej @@ -0,0 +1,59 @@ +diff a/src/icb_elem.F90 b/src/icb_elem.F90 (rejected hunks) +@@ -376,9 +376,6 @@ 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 + +@@ -398,30 +395,21 @@ type(t_partit), intent(inout), target :: partit + + 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 ++ if( (use_cavity) .and. (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 + end if + call com_integer(partit, i_have_element,elem) !max 1 PE sends element here; + 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 + +@@ -461,14 +449,11 @@ do m=1, 3 + 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 ++ if( (use_cavity) .and. (reject_elem(mesh, partit, 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 +- + RETURN + end if + end do diff --git a/src/icb_modules.F90 b/src/icb_modules.F90 index e9c0d4260..2b7aa088c 100644 --- a/src/icb_modules.F90 +++ b/src/icb_modules.F90 @@ -2,6 +2,9 @@ module iceberg_params use MOD_PARTIT USE MOD_MESH use g_config, only: use_cavity, use_cavityonelem +use, intrinsic :: iso_fortran_env, only: real64 +USE o_PARAM, only: WP + implicit none save !integer,parameter :: ib_num ! realistic dataset comprising 6912 icebergs @@ -77,7 +80,7 @@ module iceberg_params character(100):: scaling_file='icb_scaling.dat' !scaling factor !===== OUTPUT RELATED SETTINGS ===== - integer :: icb_outfreq = 180 ! 180; for FESOM_dt=2min this is 6 hourly output !120; for FESOM_dt=3min this is 6 hourly output + integer :: icb_outfreq ! 180; for FESOM_dt=2min this is 6 hourly output !120; for FESOM_dt=3min this is 6 hourly output logical :: l_geo_out = .true. ! output in unrotated (.true.) or rotated coordinates logical :: ascii_out = .false. ! old ascii output (slow, more detailed); false: faster nc output !===== NUMERICS (DONT HAVE TO BE CHANGED) ===== @@ -90,11 +93,13 @@ module iceberg_params logical,dimension(:), allocatable:: find_iceberg_elem real,dimension(:), allocatable:: f_u_ib_old, f_v_ib_old real,dimension(:), allocatable:: bvl_mean, lvlv_mean, lvle_mean, lvlb_mean !averaged volume losses - !real,dimension(:), allocatable:: fw_flux_ib, heat_flux_ib - real,dimension(:), allocatable:: fwe_flux_ib, fwl_flux_ib, fwb_flux_ib, fwbv_flux_ib, heat_flux_ib, lheat_flux_ib + !real,dimension(:), allocatable:: fw_flux_ib, hfb_flux_ib + real,dimension(:), allocatable:: fwe_flux_ib, fwl_flux_ib, fwb_flux_ib, fwbv_flux_ib + real,dimension(:), allocatable:: hfe_flux_ib, hfl_flux_ib, hfb_flux_ib, hfbv_flux_ib, lhfb_flux_ib !===== 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 real,dimension(:), allocatable:: ibfwb !freshwater flux into ocean from basal melting real,dimension(:), allocatable:: ibfwbv !freshwater flux into ocean from basal melting real,dimension(:), allocatable:: ibfwl !freshwater flux into ocean from lateral melting @@ -107,6 +112,8 @@ module iceberg_params !for communication real,dimension(:), allocatable:: arr_block integer,dimension(:), allocatable:: elem_block + integer,dimension(:), allocatable:: pe_block + real(real64), dimension(:), allocatable:: elem_area_glob real,dimension(:), allocatable:: vl_block !array for output in netcdf @@ -119,8 +126,8 @@ module iceberg_params contains ! true if all nodes of the element are either "real" model boundary nodes or shelf nodes logical function reject_elem(mesh, partit, elem) - -integer, intent(in) :: elem + implicit none + integer, intent(in) :: elem type(t_mesh), intent(in) , target :: mesh type(t_partit), intent(inout), target :: partit #include "associate_mesh_def.h" @@ -128,12 +135,9 @@ logical function reject_elem(mesh, partit, elem) #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? 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 @@ -145,8 +149,8 @@ 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, partit, elem) - -integer, intent(in) :: elem + implicit none + integer, intent(in) :: elem type(t_mesh), intent(in) , target :: mesh type(t_partit), intent(inout), target :: partit #include "associate_mesh_def.h" diff --git a/src/icb_modules.F90.rej b/src/icb_modules.F90.rej new file mode 100644 index 000000000..469d257ab --- /dev/null +++ b/src/icb_modules.F90.rej @@ -0,0 +1,72 @@ +diff a/src/icb_modules.F90 b/src/icb_modules.F90 (rejected hunks) +@@ -1,4 +1,10 @@ + module iceberg_params ++use MOD_PARTIT ++USE MOD_MESH ++use g_config, only: use_cavity, use_cavityonelem ++use, intrinsic :: iso_fortran_env, only: real64 ++USE o_PARAM, only: WP ++ + implicit none + save + !integer,parameter :: ib_num ! realistic dataset comprising 6912 icebergs +@@ -111,35 +121,50 @@ save + 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 ++ logical function reject_elem(mesh, partit, elem) + implicit none + integer, intent(in) :: elem + 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 (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 ++ 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 ++ integer function coastal_nodes(mesh, partit, elem) + implicit none + integer, intent(in) :: elem + 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 (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 diff --git a/src/icb_step.F90 b/src/icb_step.F90 index 735600e90..ad3af496b 100644 --- a/src/icb_step.F90 +++ b/src/icb_step.F90 @@ -11,7 +11,6 @@ module iceberg_step public :: iceberg_calculation public :: iceberg_step1 - public :: get_total_iceberg_area public :: iceberg_step2 public :: initialize_velo public :: trajectory @@ -53,12 +52,14 @@ subroutine iceberg_calculation(ice, mesh, partit, dynamics, istep) real(kind=8) :: t0, t1, t2, t3, t4, t0_restart, t1_restart != logical :: firstcall=.true. != logical :: lastsubstep != - != + real :: arr_from_block(15) != integer :: elem_from_block != real :: vl_from_block(4) != real,dimension(15*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 != type(t_ice), intent(inout), target :: ice @@ -73,7 +74,6 @@ subroutine iceberg_calculation(ice, mesh, partit, dynamics, istep) ! kh 16.03.21 (asynchronous) iceberg computation starts with the content in common arrays at istep and will merge its results at istep_end_synced istep_end_synced = istep + steps_per_ib_step - 1 - if(firstcall) then !overwrite icb_modules if restart, initialize netcdf output if no restart: @@ -97,6 +97,7 @@ subroutine iceberg_calculation(ice, mesh, partit, dynamics, istep) !faster communication via ALLREDUCE arr_block = 0.0 elem_block = 0 + pe_block = 0 !for communication of averaged volume losses vl_block = 0.0 @@ -138,6 +139,7 @@ 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 @@ -162,6 +164,29 @@ subroutine iceberg_calculation(ice, mesh, partit, dynamics, istep) !$omp end critical end do +!$omp critical + call MPI_IAllREDUCE(pe_block, pe_block_red, ib_num, MPI_INTEGER, MPI_SUM, partit%MPI_COMM_FESOM_IB, req, partit%MPIERR_IB) +!$omp end critical + +completed = .false. + do while (.not. completed) +!$omp critical + CALL MPI_TEST(req, completed, status, partit%MPIERR_IB) +!$omp end critical + end do + +!!$omp critical +! call MPI_IAllREDUCE(elem_area_block, elem_area_block_red, ib_num, MPI_REAL, MPI_SUM, partit%MPI_COMM_FESOM_IB, req, partit%MPIERR_IB) +!!$omp end critical +! +!completed = .false. +! do while (.not. completed) +!!$omp critical +! CALL MPI_TEST(req, completed, status, partit%MPIERR_IB) +!!$omp end critical +! end do + + !$omp critical call MPI_IAllREDUCE(vl_block, vl_block_red, 4*ib_num, MPI_DOUBLE_PRECISION, MPI_SUM, partit%MPI_COMM_FESOM_IB, req, partit%MPIERR_IB) !$omp end critical @@ -206,7 +231,8 @@ subroutine iceberg_calculation(ice, mesh, partit, dynamics, istep) conc_sill(ib),P_sill(ib), rho_h2o(ib),rho_air(ib),rho_ice(ib), & u_ib(ib),v_ib(ib), iceberg_elem(ib), find_iceberg_elem(ib), lastsubstep,& f_u_ib_old(ib), f_v_ib_old(ib), l_semiimplicit, & - semiimplicit_coeff, AB_coeff, istep) + semiimplicit_coeff, AB_coeff, istep, elem_block_red, & + pe_block_red) end if end if end do @@ -243,6 +269,7 @@ subroutine iceberg_calculation(ice, mesh, partit, dynamics, istep) write(*,*) 'reading restart took', t1_restart-t0_restart write(*,*) '*************************************************************' end if + end subroutine iceberg_calculation @@ -250,7 +277,7 @@ end subroutine iceberg_calculation !**************************************************************************************************************************** -subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib,length_ib,width_ib, lon_deg,lat_deg, & +subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib_single,length_ib_single,width_ib_single, lon_deg,lat_deg, & Co,Ca,Ci, Cdo_skin,Cda_skin, rho_icb, & conc_sill,P_sill, rho_h2o,rho_air,rho_ice, & u_ib,v_ib, iceberg_elem, find_iceberg_elem, & @@ -264,8 +291,7 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib,length_ib,wi use g_rotate_grid !for subroutine g2r, logfile_outfreq != use g_config, only: steps_per_ib_step != -!#ifdef use_cavity -! use iceberg_params, only: 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, scaling !, 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 @@ -274,7 +300,7 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib,length_ib,wi integer, intent(in) :: ib, istep - real, intent(inout) :: height_ib,length_ib,width_ib + real, intent(inout) :: height_ib_single,length_ib_single,width_ib_single real, intent(inout) :: lon_deg,lat_deg real, intent(in) :: Co,Ca,Ci, Cdo_skin,Cda_skin real, intent(in) :: rho_icb, conc_sill,P_sill, rho_h2o,rho_air,rho_ice @@ -324,7 +350,6 @@ 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 @@ -334,17 +359,20 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib,length_ib,wi #include "associate_mesh_def.h" #include "associate_part_ass.h" #include "associate_mesh_ass.h" - + + +!write(*,*) "LA DEBUG: 1" mype =>partit%mype istep_end_synced = istep + steps_per_ib_step - 1 - depth_ib = -height_ib * rho_icb/rho_h2o - volume_ib= length_ib * width_ib * height_ib + depth_ib = -height_ib_single * rho_icb/rho_h2o + volume_ib= length_ib_single * width_ib_single * height_ib_single mass_ib = volume_ib * rho_icb !less mass lon_rad = lon_deg*rad lat_rad = lat_deg*rad +!write(*,*) "LA DEBUG: 2" if(volume_ib .le. smallestvol_icb) then melted(ib) = .true. @@ -354,6 +382,7 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib,length_ib,wi return end if +!write(*,*) "LA DEBUG: 3" if (firstcall) then if(mype==0) write(*,*) 'Preparing local_idx_of array...' @@ -363,6 +392,7 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib,length_ib,wi firstcall=.false. if(mype==0) write(*,*) 'Preparing local_idx_of done.' end if +!write(*,*) "LA DEBUG: 4" if (find_iceberg_elem) then lon_rad = lon_deg*rad @@ -375,25 +405,16 @@ 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 .true. + i_have_element= (iceberg_elem .ne. 0) !up to 3 PEs possible if(i_have_element) then i_have_element= mesh%elem2D_nodes(1,iceberg_elem) <= partit%myDim_nod2D !1 PE still .true. - 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 + if( (use_cavity) .and. (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 end if call com_integer(partit, i_have_element,iceberg_elem) @@ -404,7 +425,9 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib,length_ib,wi end if ! initialize the iceberg velocity +!write(*,*) "LA DEBUG: 4d" call initialize_velo(mesh, partit, dynamics, i_have_element, ib, u_ib, v_ib, lon_rad, lat_rad, depth_ib, local_idx_of(iceberg_elem)) +!write(*,*) "LA DEBUG: 4e" !iceberg elem of ib is found find_iceberg_elem = .false. @@ -421,8 +444,10 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib,length_ib,wi endif endif +!write(*,*) "LA DEBUG: 4f" end if +!write(*,*) "LA DEBUG: 5" ! ================== START ICEBERG CALCULATION ==================== @@ -432,7 +457,10 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib,length_ib,wi !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 if( local_idx_of(iceberg_elem) > 0 ) then +write(*,*) "LA DEBUG: 5a" + if( elem2D_nodes(1,local_idx_of(iceberg_elem)) <= partit%myDim_nod2D ) then +write(*,*) "LA DEBUG: 5b" i_have_element=.true. @@ -444,14 +472,13 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib,length_ib,wi 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, length_ib, width_ib, local_idx_of(iceberg_elem), & + height_ib_single, length_ib_single, width_ib_single, local_idx_of(iceberg_elem), & mass_ib, Ci, Ca, Co, Cda_skin, Cdo_skin, & rho_ice, rho_air, rho_h2o, P_sill,conc_sill, frozen_in, & file_forces_u, file_forces_v, P_ib, conci_ib, & dt*REAL(steps_per_ib_step), l_output, f_u_ib_old, & f_v_ib_old, l_semiimplicit, semiimplicit_coeff, & AB_coeff, file_meltrates, rho_icb) - call prepare_icb2fesom(mesh,partit,ib,i_have_element,local_idx_of(iceberg_elem),depth_ib) dudt = (new_u_ib-u_ib)/REAL(steps_per_ib_step) / dt dvdt = (new_v_ib-v_ib)/REAL(steps_per_ib_step) / dt @@ -465,6 +492,7 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib,length_ib,wi !write(*,*) 'depth at iceberg ', ib, 's location:', Zdepth !=================CHECK IF ICEBERG IS GROUNDED...=================== + old_element = iceberg_elem !save if iceberg left model domain if((draft_scale(ib)*abs(depth_ib) .gt. Zdepth) .and. l_allowgrounding ) then !if((draft_scale(ib)*abs(depth_ib) .gt. minval(Zdepth3)) .and. l_allowgrounding ) then !icebergs remains stationary (iceberg can melt above in iceberg_dyn!) @@ -473,25 +501,20 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib,length_ib,wi v_ib = 0.0 old_lon = lon_rad old_lat = lat_rad - - !!########################################### - !! LA: prevent too many icebergs in one element - old_element = iceberg_elem !save if iceberg left model domain - !!########################################### ! 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 - if (mod(istep_end_synced,logfile_outfreq)==0) then - write(*,*) 'iceberg ib ', ib, 'is grounded' - grounded(ib) = .true. - end if + grounded(ib) = .true. + !if (mod(istep_end_synced,logfile_outfreq)==0) then +! if(ib==148) write(*,*) "LA DEBUG: 3 - elem ", iceberg_elem + write(*,*) 'iceberg ib ', ib, 'is grounded' + !end if else !===================...ELSE CALCULATE TRAJECTORY==================== - ! LA: prevent too many icebergs in one element - old_element = iceberg_elem !save if iceberg left model domain t0=MPI_Wtime() +!write(*,*) "LA DEBUG: before_trajectory" call trajectory( lon_rad,lat_rad, u_ib,v_ib, new_u_ib,new_v_ib, & lon_deg,lat_deg,old_lon,old_lat, dt*REAL(steps_per_ib_step)) @@ -518,6 +541,7 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib,length_ib,wi iceberg_elem=local_idx_of(iceberg_elem) !local t7=MPI_Wtime() call find_new_iceberg_elem(mesh,partit, iceberg_elem, (/lon_deg, lat_deg/), left_mype) + t8=MPI_Wtime() iceberg_elem=partit%myList_elem2D(iceberg_elem) !global end if @@ -526,35 +550,53 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib,length_ib,wi !----------------------------- ! LA 2022-11-30 - do idx = 1, size(elem_block) - if (elem_block(idx) == iceberg_elem) then - area_ib_tot = area_ib_tot + length_ib * width_ib * scaling(idx) - end if - end do +!write(*,*) "LA DEBUG: before_saturation" + if(lcell_saturation) then +!write(*,*) "LA DEBUG: start_saturation" + area_ib_tot = length_ib_single*width_ib_single*scaling(ib) +!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(idx, area_ib_tot) +!$OMP DO + 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 +!write(*,*) "LA DEBUG: end_loop_saturation" +!$OMP END DO +!$OMP END PARALLEL !----------------------------- - if((area_ib_tot > elem_area(local_idx_of(iceberg_elem))) .and. & - (iceberg_elem .ne. old_element) .and. & - (old_element .ne. 0) .and. & - (.not.grounded(ib))) then - lon_rad = old_lon - lat_rad = old_lat - lon_deg = lon_rad/rad - lat_deg = lat_rad/rad - iceberg_elem = old_element - u_ib = 0. - v_ib = 0. + if ((area_ib_tot > elem_area_glob(iceberg_elem)) .and. (old_element.ne.0) .and. (left_mype == 0)) 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)) + i_have_element = .true. + left_mype = 0.0 + lon_rad = old_lon + lat_rad = old_lat + lon_deg = lon_rad/rad + lat_deg = lat_rad/rad + iceberg_elem = old_element + u_ib = 0. + v_ib = 0. + end if +!write(*,*) "LA DEBUG: after_cell_saturation" end if !########################################### !values for communication - arr= (/ height_ib,length_ib,width_ib, u_ib,v_ib, lon_rad,lat_rad, & + 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/) !save in larger array arr_block((ib-1)*15+1 : ib*15)=arr elem_block(ib)=iceberg_elem - + pe_block(ib)=mype + + volume_ib=height_ib_single*length_ib_single*width_ib_single +!write(*,*) "LA DEBUG: before_prepare" + call prepare_icb2fesom(mesh,partit,ib,i_have_element,local_idx_of(iceberg_elem),depth_ib) +!write(*,*) "LA DEBUG: after_prepare" end if !processor has element? end if !... and first node belongs to processor? @@ -571,54 +613,25 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib,length_ib,wi end subroutine iceberg_step1 -subroutine get_total_iceberg_area(mesh, partit,iceberg_elem, area_ib_tot) - - use o_param !for rad != - USE MOD_MESH - use MOD_PARTIT !for myDim_elem2D, myList_nod2D != - use g_rotate_grid !for subroutine g2r, logfile_outfreq != - !use iceberg_params, only: arr_block, elem_block, length_ib, width_ib, scaling - - implicit none != - - integer, intent(inout) :: iceberg_elem !global - real, intent(inout) :: area_ib_tot - integer :: idx - -type(t_mesh), intent(in) , target :: mesh -type(t_partit), intent(inout), target :: partit -!========================= MODULES & DECLARATIONS =====================================!= -#include "associate_part_def.h" -#include "associate_mesh_def.h" -#include "associate_part_ass.h" -#include "associate_mesh_ass.h" - area_ib_tot = 0.0 - 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 - !########################################### -end subroutine get_total_iceberg_area - -subroutine iceberg_step2(mesh, partit,arr, elem_from_block, ib, height_ib,length_ib,width_ib, lon_deg,lat_deg, & +subroutine iceberg_step2(mesh, partit,arr, elem_from_block, ib, height_ib_single, length_ib_single, width_ib_single, lon_deg,lat_deg, & Co,Ca,Ci, Cdo_skin,Cda_skin, rho_icb, & conc_sill,P_sill, rho_h2o,rho_air,rho_ice, & u_ib,v_ib, iceberg_elem, find_iceberg_elem, & lastsubstep, f_u_ib_old, & f_v_ib_old, l_semiimplicit, semiimplicit_coeff, & - AB_coeff, istep) + AB_coeff, istep, elem_block_red, pe_block_red) !============================= MODULES & DECLARATIONS =========================================!= != use o_param !for rad != use g_rotate_grid !for subroutine g2r, logfile_outfreq != use g_config, only: steps_per_ib_step + use g_comm_auto != -!#ifdef use_cavity -! 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, reject_elem, elem_block +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 !#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 @@ -628,7 +641,7 @@ subroutine iceberg_step2(mesh, partit,arr, elem_from_block, ib, height_ib,length real, intent(in) :: arr(15) integer, intent(in) :: elem_from_block integer, intent(in) :: ib - real, intent(inout) :: height_ib,length_ib,width_ib + real, intent(inout) :: height_ib_single, length_ib_single, width_ib_single real, intent(inout) :: lon_deg,lat_deg real, intent(in) :: Co,Ca,Ci, Cdo_skin,Cda_skin real, intent(in) :: rho_icb, conc_sill,P_sill, rho_h2o,rho_air,rho_ice @@ -640,7 +653,8 @@ subroutine iceberg_step2(mesh, partit,arr, elem_from_block, ib, height_ib,length logical, intent(in) :: l_semiimplicit real, intent(in) :: semiimplicit_coeff real, intent(in) :: AB_coeff - + integer, intent(in), dimension(ib_num):: elem_block_red != + integer, intent(in), dimension(ib_num):: pe_block_red != integer, dimension(:), save, allocatable :: local_idx_of real :: depth_ib, volume_ib, mass_ib @@ -657,9 +671,11 @@ subroutine iceberg_step2(mesh, partit,arr, elem_from_block, ib, height_ib,length integer :: istep_end_synced ! LA: add threshold for number of icebergs in one elemt + integer status(MPI_STATUS_SIZE) integer :: num_ib_in_elem, idx real :: area_ib_tot - real :: local_elem_area + !real(real64), dimension(:), allocatable :: rbuffer, local_elem_area + real(real64) :: elem_area_tmp !iceberg output character :: ib_char*10 @@ -705,9 +721,9 @@ subroutine iceberg_step2(mesh, partit,arr, elem_from_block, ib, height_ib,length iceberg_elem= elem_from_block !update element as before in com_values old_element = elem_from_block !save if iceberg left model domain - height_ib= arr(1) - length_ib= arr(2) - width_ib = arr(3) + height_ib_single= arr(1) + length_ib_single= arr(2) + width_ib_single = arr(3) u_ib = arr(4) v_ib = arr(5) lon_rad = arr(6) @@ -722,9 +738,9 @@ subroutine iceberg_step2(mesh, partit,arr, elem_from_block, ib, height_ib,length dvdt = arr(13) P_ib = arr(14) conci_ib = arr(15) - + !**** check if iceberg melted in step 1 ****! - volume_ib = height_ib * length_ib * width_ib ! * rho_icb + volume_ib = height_ib_single * length_ib_single * width_ib_single ! * rho_icb if(volume_ib .le. smallestvol_icb) then buoy_props(ib, :) = 0. ! for output: NaN or MissVal could be written here return @@ -742,10 +758,11 @@ subroutine iceberg_step2(mesh, partit,arr, elem_from_block, ib, height_ib,length call global2local(mesh, partit, local_idx_of, elem2D) firstcall=.false. end if - + if(left_mype > 0.) then call iceberg_elem4all(mesh, partit, iceberg_elem, lon_deg, lat_deg) !Just PE changed? if(iceberg_elem == 0 ) then + left_mype = 0.0 lon_rad = old_lon lat_rad = old_lat lon_deg = lon_rad/rad @@ -753,18 +770,38 @@ subroutine iceberg_step2(mesh, partit,arr, elem_from_block, ib, height_ib,length iceberg_elem = old_element u_ib = 0. v_ib = 0. - else + else if(lcell_saturation) then if (mype==0) write(*,*) 'iceberg ',ib, ' changed PE or was very fast' - call get_total_iceberg_area(mesh, partit, iceberg_elem, area_ib_tot) - if(area_ib_tot > elem_area(local_idx_of(iceberg_elem))) then + elem_area_tmp = elem_area_glob(iceberg_elem) + area_ib_tot = length_ib_single * width_ib_single * scaling(ib) +!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(idx, area_ib_tot) +!$OMP DO + 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 + if((area_ib_tot > elem_area_tmp) .and. (elem_area_tmp > 0.0)) then + if(mype==pe_block_red(ib)) then + write(*,*) " *******************************************************" + write(*,*) " * iceberg changed PE and saturation" + write(*,*) " * set iceberg ", ib, " back to elem ", old_element, " from elem ", iceberg_elem + write(*,*) " * area_ib = ", length_ib_single * width_ib_single, "; area_ib_tot = ", area_ib_tot, "; elem_area = ", elem_area_tmp + end if + left_mype = 0.0 lon_rad = old_lon - lat_rad = old_lat + lat_rad = old_lat lon_deg = lon_rad/rad lat_deg = lat_rad/rad iceberg_elem = old_element u_ib = 0. - v_ib = 0. + v_ib = 0. end if + else if(lcell_saturation) then + if (mype==0) write(*,*) 'iceberg ',ib, ' changed PE or was very fast' end if end if @@ -826,9 +863,9 @@ subroutine iceberg_step2(mesh, partit,arr, elem_from_block, ib, height_ib,length 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 + buoy_props(ib,10) = height_ib_single + buoy_props(ib,11) = length_ib_single + buoy_props(ib,12) = width_ib_single buoy_props(ib,13) = iceberg_elem ! end if @@ -836,14 +873,6 @@ subroutine iceberg_step2(mesh, partit,arr, elem_from_block, ib, height_ib,length end if t4=MPI_Wtime() - -! if (mod(istep,logfile_outfreq)==0 .and. mype==0 .and. lastsubstep) then -! write(*,*) '*** step2 ***' -! write(*,*) 'comvalues took', t2-t1 -! write(*,*) 'left mype took', t3-t2 -! write(*,*) 'track out took', t4-t3 -! end if - end subroutine iceberg_step2 @@ -994,6 +1023,7 @@ subroutine parallel2coast(mesh, partit,u, v, lon,lat, elem) real, intent(in) :: lon, lat !radiant integer, intent(in) :: elem + integer :: fld_tmp integer, dimension(3) :: n integer :: node, m, i real, dimension(2) :: velocity, velocity1, velocity2 @@ -1006,12 +1036,13 @@ subroutine parallel2coast(mesh, partit,u, v, lon,lat, elem) #include "associate_part_ass.h" #include "associate_mesh_ass.h" -!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 + if( use_cavity ) then + fld_tmp = coastal_nodes(mesh, partit, elem) + else + fld_tmp = sum( mesh%bc_index_nod2D(mesh%elem2D_nodes(:,elem)) ) + end if + + SELECT CASE ( fld_tmp ) !num of coastal points CASE (0) !...coastal points: do nothing return @@ -1022,18 +1053,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) - 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 + 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 - endif + if( mesh%bc_index_nod2D(node)==1 ) then + n(i) = node + exit + end if + end if end do !write(*,*) 'one coastal node ', n(1) @@ -1045,7 +1076,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. (mesh%cavity_flag_n(node)==1) ) ) then + ! if ( (node /= n(1)) .and. ( (bc_index_nod2D(node)==1) .OR. (cavity_flag_nod2d(node)==1) ) ) then !#else ! if ( (node /= n(1)) .and. (bc_index_nod2D(node)==1)) then !#endif @@ -1081,18 +1112,18 @@ subroutine parallel2coast(mesh, partit,u, v, lon,lat, elem) velocity = [ u, v ] do m = 1, 3 node = mesh%elem2D_nodes(m,elem) -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 + 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 + i = i+1 + end if + else + if( mesh%bc_index_nod2D(node)==1 ) then + n(i) = node + i = i+1 + end if end if -endif end do call projection(mesh,partit, velocity, n(1), n(2)) @@ -1654,7 +1685,7 @@ subroutine init_buoy_output(partit) longname='time' ! use NetCDF Climate and Forecast (CF) Metadata Convention status = nf_PUT_ATT_TEXT(ncid, time_varid, 'long_name', len_trim(longname), trim(longname)) if (status .ne. nf_noerr) call handle_err(status) - write(att_text, '(a14,I4.4,a1,I2.2,a1,I2.2,a6)') 'seconds since ', year_start, '-', month_start, '-', day_start, ' 00:00:00' + write(att_text, '(a14,I4.4,a1,I2.2,a1,I2.2,a6)'), 'seconds since ', year_start, '-', month_start, '-', day_start, ' 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) if (include_fleapyear) then diff --git a/src/icb_step.F90.rej b/src/icb_step.F90.rej new file mode 100644 index 000000000..f9c2a01af --- /dev/null +++ b/src/icb_step.F90.rej @@ -0,0 +1,121 @@ +diff a/src/icb_step.F90 b/src/icb_step.F90 (rejected hunks) +@@ -376,22 +407,23 @@ type(t_dyn) , intent(inout), target :: dynamics + 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 ++!write(*,*) "LA DEBUG: 4a" + + if(i_have_element) then ++!write(*,*) "LA DEBUG: 4a1" + 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 ++!write(*,*) "LA DEBUG: 4a2" ++ if( (use_cavity) .and. (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 ++!write(*,*) "LA DEBUG: 4a3" + end if ++!write(*,*) "LA DEBUG: 4b" + call com_integer(partit, i_have_element,iceberg_elem) ++!write(*,*) "LA DEBUG: 4c" + + if(iceberg_elem .EQ. 0) then + write(*,*) 'IB ',ib,' rot. coords:', lon_deg, lat_deg !,lon_rad, lat_rad +@@ -983,15 +1023,14 @@ 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 + real, intent(in) :: lon, lat !radiant + integer, intent(in) :: elem + ++ integer :: fld_tmp + integer, dimension(3) :: n + integer :: node, m, i + real, dimension(2) :: velocity, velocity1, velocity2 +@@ -1004,12 +1043,13 @@ type(t_partit), intent(inout), target :: partit + #include "associate_part_ass.h" + #include "associate_mesh_ass.h" + +-#ifdef use_cavity +- SELECT CASE ( coastal_nodes(mesh, 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 ++ if( use_cavity ) then ++ fld_tmp = coastal_nodes(mesh, partit, elem) ++ else ++ fld_tmp = sum( mesh%bc_index_nod2D(mesh%elem2D_nodes(:,elem)) ) ++ end if ++ ++ SELECT CASE ( fld_tmp ) !num of coastal points + CASE (0) !...coastal points: do nothing + return + +@@ -1020,14 +1060,18 @@ type(t_partit), intent(inout), target :: partit + 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. 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)==1 ) then ++ n(i) = node ++ exit ++ end if ++ end if + end do + + !write(*,*) 'one coastal node ', n(1) +@@ -1075,13 +1119,17 @@ type(t_partit), intent(inout), target :: partit + 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 +- n(i) = node +- i = i+1 ++ 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 ++ i = i+1 ++ end if ++ else ++ if( mesh%bc_index_nod2D(node)==1 ) then ++ n(i) = node ++ i = i+1 ++ end if + end if + end do + call projection(mesh,partit, velocity, n(1), n(2)) diff --git a/src/icb_thermo.F90 b/src/icb_thermo.F90 index edf124e6c..57cb85f63 100644 --- a/src/icb_thermo.F90 +++ b/src/icb_thermo.F90 @@ -1,4 +1,4 @@ -!============================================================================== +!!============================================================================= ! calculates the empirical melt rates of the iceberg as in ! Martin: 'Parameterizing the fresh-water flux from land ice to ocean ! with interactive icebergs in a coupled climate model'(2010) @@ -16,51 +16,61 @@ ! use 3D information for T,S and velocities ! instead of SSTs; M_v depends on 'thermal driving') !============================================================================== -subroutine iceberg_meltrates( M_b, M_v, M_e, M_bv, & +subroutine iceberg_meltrates(partit, M_b, M_v, M_e, M_bv, & u_ib,v_ib, uo_ib,vo_ib, ua_ib,va_ib, & sst_ib, length_ib, conci_ib, & - uo_keel_ib, vo_keel_ib, T_keel_ib, S_keel_ib, depth_ib, & - T_ave_ib, S_ave_ib, ib) + uo_keel_ib, vo_keel_ib, T_keel_ib, S_keel_ib, depth_ib, height_ib, & + T_ave_ib, S_ave_ib, ib, rho_icb) use o_param + use MOD_PARTIT use g_clock use g_forcing_arrays use g_rotate_grid - use iceberg_params, only: fwe_flux_ib, fwl_flux_ib, fwb_flux_ib, fwbv_flux_ib, heat_flux_ib - + 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 implicit none + ! LA: include latent heat 2023-04-04 + real(kind=8),parameter :: L = 334000. ! [J/Kg] + real, intent(IN) :: u_ib,v_ib, uo_ib,vo_ib, ua_ib,va_ib !iceberg velo, (int.) ocean & atm velo real, intent(IN) :: uo_keel_ib, vo_keel_ib !ocean velo at iceberg's draft - real, intent(IN) :: sst_ib, length_ib, conci_ib !SST, length and sea ice conc. - real, intent(IN) :: T_keel_ib, S_keel_ib, depth_ib !T & S at depth 'depth_ib' + real, intent(IN) :: sst_ib, length_ib, conci_ib, rho_icb !SST, length and sea ice conc. + real, intent(IN) :: T_keel_ib, S_keel_ib, depth_ib, height_ib !T & S at depth 'depth_ib' real, intent(IN) :: T_ave_ib, S_ave_ib !T & S averaged, i.e. at 'depth_ib/2' integer, intent(IN) :: ib !iceberg ID real, intent(OUT) :: M_b, M_v, M_e, M_bv !melt rates [m (ice) per s] + real :: H_b, H_v, H_e, H_bv !melt rates [m (ice) per s] real :: absamino, damping, sea_state, v_ibmino real :: tf, T_d !freezing temp. and 'thermal driving' +type(t_partit), intent(inout), target :: partit +#include "associate_part_def.h" +#include "associate_part_ass.h" !3-eq. formulation for bottom melting [m/s] v_ibmino = sqrt( (u_ib - uo_keel_ib)**2 + (v_ib - vo_keel_ib)**2 ) - call iceberg_heat_water_fluxes_3eq(ib, M_b, T_keel_ib,S_keel_ib,v_ibmino, depth_ib, tf) + call iceberg_heat_water_fluxes_3eq(partit, ib, M_b, H_b, T_keel_ib,S_keel_ib,v_ibmino, depth_ib, tf) + hfb_flux_ib = H_b * length_ib*length_ib*scaling(ib) !3-eq. formulation for lateral 'basal' melting [m/s] v_ibmino = sqrt( (u_ib - uo_ib)**2 + (v_ib - vo_ib)**2 ) ! depth-average rel. velocity - call iceberg_heat_water_fluxes_3eq(ib, M_bv, T_ave_ib,S_ave_ib,v_ibmino, depth_ib/2.0, tf) + call iceberg_heat_water_fluxes_3eq(partit, ib, M_bv, H_bv, T_ave_ib,S_ave_ib,v_ibmino, depth_ib/2.0, tf) + hfbv_flux_ib = H_bv * (2*length_ib*abs(depth_ib) + 2*length_ib*abs(depth_ib) ) * scaling(ib) !'thermal driving', defined as the elevation of ambient water !temperature above freezing point' (Neshyba and Josberger, 1979). T_d = T_ave_ib - tf if(T_d < 0.) T_d = 0. - !write(*,*) 'thermal driving:',T_d,'; Tf:',tf,'T_ave:',T_ave_ib !lateral melt (buoyant convection) !M_v is a function of the 'thermal driving', NOT just sst! Cf. Neshyba and Josberger (1979) M_v = 0.00762 * T_d + 0.00129 * T_d**2 M_v = M_v/86400. + H_v = M_v * rho_icb * L + hfl_flux_ib = H_v * (2*length_ib*abs(depth_ib) + 2*length_ib*abs(depth_ib) ) * scaling(ib) !wave erosion absamino = sqrt( (ua_ib - uo_ib)**2 + (va_ib - vo_ib)**2 ) @@ -68,6 +78,8 @@ subroutine iceberg_meltrates( 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. + H_e = M_e * rho_icb * L + hfe_flux_ib = H_e * (length_ib*abs(height_ib) + length_ib*abs(height_ib) ) * scaling(ib) !fwe_flux_ib = M_e end subroutine iceberg_meltrates @@ -90,7 +102,7 @@ subroutine iceberg_newdimensions(partit, ib, depth_ib,height_ib,length_ib,width_ use g_clock use g_forcing_arrays use g_rotate_grid - use iceberg_params, only: l_weeksmellor, ascii_out, icb_outfreq, vl_block, bvl_mean, lvlv_mean, lvle_mean, lvlb_mean, smallestvol_icb, fwb_flux_ib, fwe_flux_ib, fwbv_flux_ib, fwl_flux_ib, scaling, heat_flux_ib, lheat_flux_ib + use iceberg_params, only: l_weeksmellor, ascii_out, icb_outfreq, vl_block, bvl_mean, lvlv_mean, lvle_mean, lvlb_mean, smallestvol_icb, fwb_flux_ib, fwe_flux_ib, fwbv_flux_ib, fwl_flux_ib, scaling, hfb_flux_ib, hfbv_flux_ib, hfe_flux_ib, hfl_flux_ib, lhfb_flux_ib use g_config, only: steps_per_ib_step implicit none @@ -98,6 +110,7 @@ subroutine iceberg_newdimensions(partit, ib, depth_ib,height_ib,length_ib,width_ integer, intent(IN) :: ib real, intent(INOUT) :: depth_ib, height_ib, length_ib, width_ib real, intent(IN) :: M_b, M_v, M_e, M_bv, rho_h2o, rho_icb + real :: hf_tot_tmp 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 @@ -140,7 +153,7 @@ subroutine iceberg_newdimensions(partit, ib, depth_ib,height_ib,length_ib,width_ if((tvl .ge. volume_before) .OR. (volume_before .le. smallestvol_icb)) then volume_after=0.0 - depth_ib = 0.0 + depth_ib = 0.0 height_ib= 0.0 length_ib= 0.0 width_ib = 0.0 @@ -170,7 +183,7 @@ 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 volume_after=0.0 - depth_ib = 0.0 + depth_ib = 0.0 height_ib= 0.0 length_ib= 0.0 width_ib = 0.0 @@ -226,9 +239,18 @@ subroutine iceberg_newdimensions(partit, ib, depth_ib,height_ib,length_ib,width_ ! ----------------------- ! LA: set iceberg heatflux at least to latent heat 2023-04-04 ! Latent heat flux at base and sides also changes lines 475/476 - lheat_flux_ib(ib) = rho_icb*L*tvl*scaling(ib)/dt/REAL(steps_per_ib_step) - if( (heat_flux_ib(ib).gt.0.0) .and. (heat_flux_ib(ib).lt.lheat_flux_ib(ib))) then - heat_flux_ib(ib)=lheat_flux_ib(ib) + ! Lateral heat flux set to latent heat and basal heat flux set to zero + if( lmin_latent_hf ) then + lhfb_flux_ib(ib) = rho_icb*L*tvl*scaling(ib)/dt/REAL(steps_per_ib_step) + + hf_tot_tmp = hfb_flux_ib(ib)+hfbv_flux_ib(ib)+hfl_flux_ib(ib)+hfe_flux_ib(ib) + + if( (hf_tot_tmp >= 0.0) .and. (hf_tot_tmp < lhfb_flux_ib(ib))) then + hfe_flux_ib(ib)=-lhfb_flux_ib(ib) * abs(hfe_flux_ib(ib)/hf_tot_tmp) + hfl_flux_ib(ib)=-lhfb_flux_ib(ib) * abs(hfl_flux_ib(ib)/hf_tot_tmp) + hfb_flux_ib(ib)=-lhfb_flux_ib(ib) * abs(hfb_flux_ib(ib)/hf_tot_tmp) + hfbv_flux_ib(ib)=-lhfb_flux_ib(ib) * abs(hfbv_flux_ib(ib)/hf_tot_tmp) + end if end if ! ----------------------- end subroutine iceberg_newdimensions @@ -263,7 +285,7 @@ end subroutine weeksmellor !*************************************************************************************************************************** !*************************************************************************************************************************** -subroutine iceberg_heat_water_fluxes_3eq(ib, M_b, T_ib,S_ib,v_rel, depth_ib, t_freeze) +subroutine iceberg_heat_water_fluxes_3eq(partit, ib, M_b, H_b, T_ib,S_ib,v_rel, depth_ib, t_freeze) ! The three-equation model of ice-shelf ocean interaction (Hellmer et al., 1997) ! Code derived from BRIOS subroutine iceshelf (which goes back to H.Hellmer's 2D ice shelf model code) ! adjusted for use in FESOM by Ralph Timmermann, 16.02.2011 @@ -276,7 +298,7 @@ subroutine iceberg_heat_water_fluxes_3eq(ib, M_b, T_ib,S_ib,v_rel, depth_ib, t_f implicit none integer, INTENT(IN) :: ib - real(kind=8),INTENT(OUT) :: M_b, t_freeze + real(kind=8),INTENT(OUT) :: M_b, H_b, t_freeze real(kind=8),INTENT(IN) :: T_ib, S_ib ! ocean temperature & salinity (at depth 'depth_ib') real(kind=8),INTENT(IN) :: v_rel, depth_ib ! relative velocity iceberg-ocean (at depth 'depth_ib') @@ -312,6 +334,10 @@ subroutine iceberg_heat_water_fluxes_3eq(ib, M_b, T_ib,S_ib,v_rel, depth_ib, t_f real(kind=8),parameter :: cpi = 152.5+7.122*(atk+tob) !Paterson:"The Physics of Glaciers" real(kind=8),parameter :: L = 334000. ! [J/Kg] +type(t_partit), intent(inout), target :: partit +!==================== MODULES & DECLARATIONS ==========================!= +#include "associate_part_def.h" +#include "associate_part_ass.h" temp = T_ib sal = S_ib @@ -404,17 +430,18 @@ subroutine iceberg_heat_water_fluxes_3eq(ib, M_b, T_ib,S_ib,v_rel, depth_ib, t_f endif 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)) - !heat_flux_ib(ib) = rhow*cpw*gat*(tin-tf)*scaling(ib) ! [W/m2] ! positive for upward - heat_flux_ib(ib) = rhow*cpw*gat*(tin-tf)*length_ib(ib)*width_ib(ib)*scaling(ib) ! [W] ! positive for upward + !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 + 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 @@ -468,112 +495,6 @@ subroutine potit_ib(ib,salz,pt,pres,rfpres,tin) return end subroutine potit_ib -! if the underlying FESOM is run without cavities, the following routines might be -! missing, so put them here: -!if (.not. use_cavity) then -! -!------------------------------------------------------------------------------------- -! -!subroutine potit(salz,pt,pres,rfpres,tin) -! ! Berechnet aus dem Salzgehalt[psu] (SALZ), der pot. Temperatur[oC] -! ! (PT) und dem Referenzdruck[dbar] (REFPRES) die in-situ Temperatur -! ! [oC] (TIN) bezogen auf den in-situ Druck[dbar] (PRES) mit Hilfe -! ! eines Iterationsverfahrens aus. -! -! integer iter -! real salz,pt,pres,rfpres,tin -! real epsi,tpmd,pt1,ptd,pttmpr -! -! data tpmd / 0.001 / -! -! epsi = 0. -! do iter=1,100 -! tin = pt+epsi -! pt1 = pttmpr(salz,tin,pres,rfpres) -! ptd = pt1-pt -! if(abs(ptd).lt.tpmd) return -! epsi = epsi-ptd -! enddo -! write(6,*) ' WARNING!' -! write(6,*) ' in-situ temperature calculation has not converged.' -! stop -! return -!end subroutine potit -! -!------------------------------------------------------------------------------------- -! -!real function pttmpr(salz,temp,pres,rfpres) -! ! Berechnet aus dem Salzgehalt/psu (SALZ), der in-situ Temperatur/degC -! ! (TEMP) und dem in-situ Druck/dbar (PRES) die potentielle Temperatur/ -! ! degC (PTTMPR) bezogen auf den Referenzdruck/dbar (RFPRES). Es wird -! ! ein Runge-Kutta Verfahren vierter Ordnung verwendet. -! ! Checkwert: PTTMPR = 36.89073 DegC -! ! fuer SALZ = 40.0 psu -! ! TEMP = 40.0 DegC -! ! PRES = 10000.000 dbar -! ! RFPRES = 0.000 dbar -! -! data ct2 ,ct3 /0.29289322 , 1.707106781/ -! data cq2a,cq2b /0.58578644 , 0.121320344/ -! data cq3a,cq3b /3.414213562, -4.121320344/ -! -! real salz,temp,pres,rfpres -! real p,t,dp,dt,q,ct2,ct3,cq2a,cq2b,cq3a,cq3b -! real adlprt -! -! p = pres -! t = temp -! dp = rfpres-pres -! dt = dp*adlprt(salz,t,p) -! t = t +0.5*dt -! q = dt -! p = p +0.5*dp -! dt = dp*adlprt(salz,t,p) -! t = t + ct2*(dt-q) -! q = cq2a*dt + cq2b*q -! dt = dp*adlprt(salz,t,p) -! t = t + ct3*(dt-q) -! q = cq3a*dt + cq3b*q -! p = rfpres -! dt = dp*adlprt(salz,t,p) -! -! pttmpr = t + (dt-q-q)/6.0 -! -!end function pttmpr -! -!------------------------------------------------------------------------------------- -! -!real function adlprt(salz,temp,pres) -! ! Berechnet aus dem Salzgehalt/psu (SALZ), der in-situ Temperatur/degC -! ! (TEMP) und dem in-situ Druck/dbar (PRES) den adiabatischen Temperatur- -! ! gradienten/(K Dbar^-1) ADLPRT. -! ! Checkwert: ADLPRT = 3.255976E-4 K dbar^-1 -! ! fuer SALZ = 40.0 psu -! ! TEMP = 40.0 DegC -! ! PRES = 10000.000 dbar -! -! real salz,temp,pres -! real s0,a0,a1,a2,a3,b0,b1,c0,c1,c2,c3,d0,d1,e0,e1,e2,ds -! -! data s0 /35.0/ -! data a0,a1,a2,a3 /3.5803E-5, 8.5258E-6, -6.8360E-8, 6.6228E-10/ -! data b0,b1 /1.8932E-6, -4.2393E-8/ -! data c0,c1,c2,c3 /1.8741E-8, -6.7795E-10, 8.7330E-12, -5.4481E-14/ -! data d0,d1 /-1.1351E-10, 2.7759E-12/ -! data e0,e1,e2 /-4.6206E-13, 1.8676E-14, -2.1687E-16/ -! -! ds = salz-s0 -! adlprt = ( ( (e2*temp + e1)*temp + e0 )*pres & -! + ( (d1*temp + d0)*ds & -! + ( (c3*temp + c2)*temp + c1 )*temp + c0 ) )*pres & -! + (b1*temp + b0)*ds + ( (a3*temp + a2)*temp + a1 )*temp + a0 -! -!END function adlprt -! -!---------------------------------------------------------------------------------------- -! -!endif - ! LA from oce_dens_press for iceberg coupling subroutine fcn_density(t,s,z,rho) diff --git a/src/icb_thermo.F90.rej b/src/icb_thermo.F90.rej new file mode 100644 index 000000000..14d566b0e --- /dev/null +++ b/src/icb_thermo.F90.rej @@ -0,0 +1,324 @@ +diff a/src/icb_thermo.F90 b/src/icb_thermo.F90 (rejected hunks) +@@ -16,45 +16,39 @@ + ! use 3D information for T,S and velocities + ! instead of SSTs; M_v depends on 'thermal driving') + !============================================================================== +-subroutine iceberg_meltrates( M_b, M_v, M_e, M_bv, & ++subroutine iceberg_meltrates(partit, M_b, M_v, M_e, M_bv, & + u_ib,v_ib, uo_ib,vo_ib, ua_ib,va_ib, & + sst_ib, length_ib, conci_ib, & + uo_keel_ib, vo_keel_ib, T_keel_ib, S_keel_ib, depth_ib, & +- T_ave_ib, S_ave_ib, ib) ++ T_ave_ib, S_ave_ib, ib, rho_icb) + +-! use o_mesh + use o_param +-! use i_therm_param +-! use i_param +-! use MOD_ICE +-! use i_arrays +-! use MOD_PARTIT +- +-! kh 18.03.21 not really used here +-! use o_arrays +- ++ use MOD_PARTIT + use g_clock + use g_forcing_arrays + use g_rotate_grid + +- use iceberg_params, only: fwe_flux_ib, fwl_flux_ib, fwb_flux_ib, fwbv_flux_ib, heat_flux_ib +- ++ 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, height_ib, scaling + implicit none + ++ ! LA: include latent heat 2023-04-04 ++ real(kind=8),parameter :: L = 334000. ! [J/Kg] ++ + real, intent(IN) :: u_ib,v_ib, uo_ib,vo_ib, ua_ib,va_ib !iceberg velo, (int.) ocean & atm velo + real, intent(IN) :: uo_keel_ib, vo_keel_ib !ocean velo at iceberg's draft +- real, intent(IN) :: sst_ib, length_ib, conci_ib !SST, length and sea ice conc. +- real, intent(IN) :: T_keel_ib, S_keel_ib, depth_ib !T & S at depth 'depth_ib' ++ real, intent(IN) :: sst_ib, length_ib, conci_ib, rho_icb !SST, length and sea ice conc. ++ real, intent(IN) :: T_keel_ib, S_keel_ib, depth_ib !T & S at depth 'depth_ib' + real, intent(IN) :: T_ave_ib, S_ave_ib !T & S averaged, i.e. at 'depth_ib/2' + integer, intent(IN) :: ib !iceberg ID + real, intent(OUT) :: M_b, M_v, M_e, M_bv !melt rates [m (ice) per s] ++ real :: H_b, H_v, H_e, H_bv !melt rates [m (ice) per s] + + + real :: absamino, damping, sea_state, v_ibmino + real :: tf, T_d !freezing temp. and 'thermal driving' +-!type(t_partit), intent(inout), target :: partit +-!#include "associate_part_def.h" +-!#include "associate_part_ass.h" ++type(t_partit), intent(inout), target :: partit ++#include "associate_part_def.h" ++#include "associate_part_ass.h" + + !bottom melt (basal turbulent melting rate) + !M_b = 0.58 * sqrt( (u_ib - uo_ib)**2 + (v_ib - vo_ib)**2 )**0.8 & +@@ -63,17 +57,18 @@ subroutine iceberg_meltrates( M_b, M_v, M_e, M_bv, & + + !3-eq. formulation for bottom melting [m/s] + v_ibmino = sqrt( (u_ib - uo_keel_ib)**2 + (v_ib - vo_keel_ib)**2 ) +- call iceberg_heat_water_fluxes_3eq(ib, M_b, T_keel_ib,S_keel_ib,v_ibmino, depth_ib, tf) ++ call iceberg_heat_water_fluxes_3eq(partit, ib, M_b, H_b, T_keel_ib,S_keel_ib,v_ibmino, depth_ib, tf) ++ hfb_flux_ib = H_b * length_ib*length_ib*scaling(ib) + + !3-eq. formulation for lateral 'basal' melting [m/s] + v_ibmino = sqrt( (u_ib - uo_ib)**2 + (v_ib - vo_ib)**2 ) ! depth-average rel. velocity +- call iceberg_heat_water_fluxes_3eq(ib, M_bv, T_ave_ib,S_ave_ib,v_ibmino, depth_ib/2.0, tf) ++ call iceberg_heat_water_fluxes_3eq(partit, ib, M_bv, H_bv, T_ave_ib,S_ave_ib,v_ibmino, depth_ib/2.0, tf) ++ hfbv_flux_ib = H_bv * (2*length_ib*abs(depth_ib) + 2*length_ib*abs(depth_ib) ) * scaling(ib) + + !'thermal driving', defined as the elevation of ambient water + !temperature above freezing point' (Neshyba and Josberger, 1979). + T_d = T_ave_ib - tf + if(T_d < 0.) T_d = 0. +- !write(*,*) 'thermal driving:',T_d,'; Tf:',tf,'T_ave:',T_ave_ib + + !lateral melt (buoyant convection) + !M_v = 0.00762 * sst_ib + 0.00129 * sst_ib**2 +@@ -81,6 +76,8 @@ subroutine iceberg_meltrates( M_b, M_v, M_e, M_bv, & + !M_v is a function of the 'thermal driving', NOT just sst! Cf. Neshyba and Josberger (1979) + M_v = 0.00762 * T_d + 0.00129 * T_d**2 + M_v = M_v/86400. ++ H_v = M_v * rho_icb * L ++ hfl_flux_ib = H_v * (2*length_ib*abs(depth_ib) + 2*length_ib*abs(depth_ib) ) * scaling(ib) + !fwl_flux_ib = M_v + + !wave erosion +@@ -106,21 +105,12 @@ end subroutine iceberg_meltrates + subroutine iceberg_newdimensions(partit, ib, depth_ib,height_ib,length_ib,width_ib,M_b,M_v,M_e,M_bv, & + rho_h2o, rho_icb, file_meltrates) + +-! use o_mesh + use o_param !for step_per_day +-! use i_therm_param +-! use i_param +-! use MOD_ICE +-! use i_arrays + use MOD_PARTIT !for mype +- +-! kh 18.03.21 not really used here +-! use o_arrays +- + use g_clock + use g_forcing_arrays + use g_rotate_grid +- use iceberg_params, only: l_weeksmellor, ascii_out, icb_outfreq, vl_block, bvl_mean, lvlv_mean, lvle_mean, lvlb_mean, smallestvol_icb, fwb_flux_ib, fwe_flux_ib, fwbv_flux_ib, fwl_flux_ib, scaling, heat_flux_ib, lheat_flux_ib ++ use iceberg_params, only: l_weeksmellor, ascii_out, icb_outfreq, vl_block, bvl_mean, lvlv_mean, lvle_mean, lvlb_mean, smallestvol_icb, fwb_flux_ib, fwe_flux_ib, fwbv_flux_ib, fwl_flux_ib, scaling, hfb_flux_ib, hfbv_flux_ib, hfe_flux_ib, hfl_flux_ib, lhfb_flux_ib + use g_config, only: steps_per_ib_step + + implicit none +@@ -185,8 +176,6 @@ type(t_partit), intent(inout), target :: partit + volume_after=volume_before-tvl + + !calculating the new iceberg dimensions +- !depth_ib = (abs(depth_ib)-dh_b)*(-1.) +- !height_ib= abs(depth_ib) * rho_h2o/rho_icb + height_ib= height_ib - dh_b + depth_ib = -height_ib * rho_icb/rho_h2o + +@@ -196,15 +185,13 @@ type(t_partit), intent(inout), target :: partit + + !distribute dh_e equally between length and width + !as in code of michael schodlok, but not dh_v? +- !length_ib= length_ib - dh_v -dh_e/2. +- !width_ib = width_ib - dh_v -dh_e/2. + + volume_after=height_ib*length_ib*width_ib + + !iceberg smaller than critical value after melting? + if (volume_after .le. smallestvol_icb) then + volume_after=0.0 +- depth_ib = 0.0 ++ depth_ib = 0.0 + height_ib= 0.0 + length_ib= 0.0 + width_ib = 0.0 +@@ -297,25 +293,20 @@ end subroutine weeksmellor + !*************************************************************************************************************************** + !*************************************************************************************************************************** + +-subroutine iceberg_heat_water_fluxes_3eq(ib, M_b, T_ib,S_ib,v_rel, depth_ib, t_freeze) ++subroutine iceberg_heat_water_fluxes_3eq(partit, ib, M_b, H_b, T_ib,S_ib,v_rel, depth_ib, t_freeze) + ! The three-equation model of ice-shelf ocean interaction (Hellmer et al., 1997) + ! Code derived from BRIOS subroutine iceshelf (which goes back to H.Hellmer's 2D ice shelf model code) + ! adjusted for use in FESOM by Ralph Timmermann, 16.02.2011 + ! adopted and modified for iceberg basal melting by Thomas Rackow, 11.06.2014 + !---------------------------------------------------------------- + +- !use o_mesh +- !use o_param +- !use o_arrays +- !use i_arrays +- !use MOD_PARTIT + use iceberg_params + use g_config + + implicit none + + integer, INTENT(IN) :: ib +- real(kind=8),INTENT(OUT) :: M_b, t_freeze ++ real(kind=8),INTENT(OUT) :: M_b, H_b, t_freeze + real(kind=8),INTENT(IN) :: T_ib, S_ib ! ocean temperature & salinity (at depth 'depth_ib') + real(kind=8),INTENT(IN) :: v_rel, depth_ib ! relative velocity iceberg-ocean (at depth 'depth_ib') + +@@ -351,22 +342,11 @@ subroutine iceberg_heat_water_fluxes_3eq(ib, M_b, T_ib,S_ib,v_rel, depth_ib, t_f + real(kind=8),parameter :: cpi = 152.5+7.122*(atk+tob) !Paterson:"The Physics of Glaciers" + + real(kind=8),parameter :: L = 334000. ! [J/Kg] ++type(t_partit), intent(inout), target :: partit ++!==================== MODULES & DECLARATIONS ==========================!= ++#include "associate_part_def.h" ++#include "associate_part_ass.h" + +- ! hemw = helium content of the glacial meltwater +- ! oomw = isotopic fractionation due to melting +- ! oofw = isotopic fractionation due to freezing +- ! hemw= 4.02*14. +- ! oomw= -30. +- ! oofw= -2.5 +- +- !n3=myDim_nod3d+eDim_nod3d +- +- !do n=1,myDim_nod2D+eDim_nod2D +- !if(cavity_flag_nod2d(n)==0) cycle +- !nk=nod3d_below_nod2d(1,n) +- !temp = tracer(nk,1) +- !sal = tracer(nk,2) +- !zice = coord_nod3d(3,nk) !(<0) + temp = T_ib + sal = S_ib + zice = depth_ib !(<0) +@@ -378,15 +358,7 @@ subroutine iceberg_heat_water_fluxes_3eq(ib, M_b, T_ib,S_ib,v_rel, depth_ib, t_f + ! Calculate or prescribe the turbulent heat and salt transfer coeff. GAT and GAS + ! velocity-dependent approach of Jenkins (1991) + +- !rt vt1 = 0.25*sqrt((u(i,j,N,lrhs)+u(i+1,j,N,lrhs))**2 +- !rt & +(v(i,j,N,lrhs)+v(i,j+1,N,lrhs))**2) +- ! if(vt1.eq.0.) vt1=0.001 +- !rt re = Hz_r(i,j,N)*ds/un !Reynolds number +- +- !vt1 = sqrt(uf(nk)*uf(nk)+uf(nk+n3)*uf(nk+n3)) ! relative velocity ice-ocean + vt1 = v_rel ! relative velocity iceberg-ocean (at depth 'depth_ib') +- +-!rt RG44030 vt1 = max(vt1,0.001) + vt1 = max(vt1,0.005) ! RG44030 + + re = 10./un !vt1*re (=velocity times length scale over kinematic viscosity) is the Reynolds number +@@ -530,112 +503,6 @@ subroutine potit_ib(ib,salz,pt,pres,rfpres,tin) + return + 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 +-! +-!------------------------------------------------------------------------------------- +-! +-!subroutine potit(salz,pt,pres,rfpres,tin) +-! ! Berechnet aus dem Salzgehalt[psu] (SALZ), der pot. Temperatur[oC] +-! ! (PT) und dem Referenzdruck[dbar] (REFPRES) die in-situ Temperatur +-! ! [oC] (TIN) bezogen auf den in-situ Druck[dbar] (PRES) mit Hilfe +-! ! eines Iterationsverfahrens aus. +-! +-! integer iter +-! real salz,pt,pres,rfpres,tin +-! real epsi,tpmd,pt1,ptd,pttmpr +-! +-! data tpmd / 0.001 / +-! +-! epsi = 0. +-! do iter=1,100 +-! tin = pt+epsi +-! pt1 = pttmpr(salz,tin,pres,rfpres) +-! ptd = pt1-pt +-! if(abs(ptd).lt.tpmd) return +-! epsi = epsi-ptd +-! enddo +-! write(6,*) ' WARNING!' +-! write(6,*) ' in-situ temperature calculation has not converged.' +-! stop +-! return +-!end subroutine potit +-! +-!------------------------------------------------------------------------------------- +-! +-!real function pttmpr(salz,temp,pres,rfpres) +-! ! Berechnet aus dem Salzgehalt/psu (SALZ), der in-situ Temperatur/degC +-! ! (TEMP) und dem in-situ Druck/dbar (PRES) die potentielle Temperatur/ +-! ! degC (PTTMPR) bezogen auf den Referenzdruck/dbar (RFPRES). Es wird +-! ! ein Runge-Kutta Verfahren vierter Ordnung verwendet. +-! ! Checkwert: PTTMPR = 36.89073 DegC +-! ! fuer SALZ = 40.0 psu +-! ! TEMP = 40.0 DegC +-! ! PRES = 10000.000 dbar +-! ! RFPRES = 0.000 dbar +-! +-! data ct2 ,ct3 /0.29289322 , 1.707106781/ +-! data cq2a,cq2b /0.58578644 , 0.121320344/ +-! data cq3a,cq3b /3.414213562, -4.121320344/ +-! +-! real salz,temp,pres,rfpres +-! real p,t,dp,dt,q,ct2,ct3,cq2a,cq2b,cq3a,cq3b +-! real adlprt +-! +-! p = pres +-! t = temp +-! dp = rfpres-pres +-! dt = dp*adlprt(salz,t,p) +-! t = t +0.5*dt +-! q = dt +-! p = p +0.5*dp +-! dt = dp*adlprt(salz,t,p) +-! t = t + ct2*(dt-q) +-! q = cq2a*dt + cq2b*q +-! dt = dp*adlprt(salz,t,p) +-! t = t + ct3*(dt-q) +-! q = cq3a*dt + cq3b*q +-! p = rfpres +-! dt = dp*adlprt(salz,t,p) +-! +-! pttmpr = t + (dt-q-q)/6.0 +-! +-!end function pttmpr +-! +-!------------------------------------------------------------------------------------- +-! +-!real function adlprt(salz,temp,pres) +-! ! Berechnet aus dem Salzgehalt/psu (SALZ), der in-situ Temperatur/degC +-! ! (TEMP) und dem in-situ Druck/dbar (PRES) den adiabatischen Temperatur- +-! ! gradienten/(K Dbar^-1) ADLPRT. +-! ! Checkwert: ADLPRT = 3.255976E-4 K dbar^-1 +-! ! fuer SALZ = 40.0 psu +-! ! TEMP = 40.0 DegC +-! ! PRES = 10000.000 dbar +-! +-! real salz,temp,pres +-! real s0,a0,a1,a2,a3,b0,b1,c0,c1,c2,c3,d0,d1,e0,e1,e2,ds +-! +-! data s0 /35.0/ +-! data a0,a1,a2,a3 /3.5803E-5, 8.5258E-6, -6.8360E-8, 6.6228E-10/ +-! data b0,b1 /1.8932E-6, -4.2393E-8/ +-! data c0,c1,c2,c3 /1.8741E-8, -6.7795E-10, 8.7330E-12, -5.4481E-14/ +-! data d0,d1 /-1.1351E-10, 2.7759E-12/ +-! data e0,e1,e2 /-4.6206E-13, 1.8676E-14, -2.1687E-16/ +-! +-! ds = salz-s0 +-! adlprt = ( ( (e2*temp + e1)*temp + e0 )*pres & +-! + ( (d1*temp + d0)*ds & +-! + ( (c3*temp + c2)*temp + c1 )*temp + c0 ) )*pres & +-! + (b1*temp + b0)*ds + ( (a3*temp + a2)*temp + a1 )*temp + a0 +-! +-!END function adlprt +-! +-!---------------------------------------------------------------------------------------- +-! +-#endif +- + + ! LA from oce_dens_press for iceberg coupling + subroutine fcn_density(t,s,z,rho) diff --git a/src/ice_oce_coupling.F90 b/src/ice_oce_coupling.F90 index 0d4575668..c33824352 100755 --- a/src/ice_oce_coupling.F90 +++ b/src/ice_oce_coupling.F90 @@ -264,6 +264,7 @@ subroutine oce_fluxes(ice, dynamics, tracers, partit, mesh) use icedrv_main, only: icepack_to_fesom, & init_flux_atm_ocn #endif + use iceberg_params use cavity_interfaces !---fwf-code use g_clock @@ -567,6 +568,18 @@ subroutine oce_fluxes(ice, dynamics, tracers, partit, mesh) where (ulevels_nod2d > 1) flux = -water_flux end if end if + + !___________________________________________________________________________ + if (use_icebergs .and. lbalance_fw) then + call icb2fesom(mesh, partit, ice) + if (.not.turn_off_fw) then + flux = flux + (ibfwb + ibfwe + ibfwl + ibfwbv) !* steps_per_ib_step + end if + + call integrate_nod(ibfwb + ibfwe + ibfwl + ibfwbv, net, partit, mesh) + if (mype==0) write(*,*) " * total iceberg fw flux: ", net + end if + !___________________________________________________________________________ ! compute total global net freshwater flux into the ocean call integrate_nod(flux, net, partit, mesh) diff --git a/src/ice_oce_coupling.F90.rej b/src/ice_oce_coupling.F90.rej new file mode 100644 index 000000000..2d095512b --- /dev/null +++ b/src/ice_oce_coupling.F90.rej @@ -0,0 +1,23 @@ +diff a/src/ice_oce_coupling.F90 b/src/ice_oce_coupling.F90 (rejected hunks) +@@ -563,9 +564,21 @@ subroutine oce_fluxes(ice, dynamics, tracers, partit, mesh) + end if + end if + ++ !___________________________________________________________________________ ++ if (use_icebergs .and. lbalance_fw) then ++ call icb2fesom(mesh, partit, ice) ++ if (.not.turn_off_fw) then ++ flux = flux + (ibfwb + ibfwe + ibfwl + ibfwbv) !* steps_per_ib_step ++ end if ++ ++ call integrate_nod(ibfwb + ibfwe + ibfwl + ibfwbv, net, partit, mesh) ++ if (mype==0) write(*,*) " * total iceberg fw flux: ", net ++ end if ++ + !___________________________________________________________________________ + ! compute total global net freshwater flux into the ocean + call integrate_nod(flux, net, partit, mesh) ++ + + !___________________________________________________________________________ + ! here the + sign must be used because we switched up the sign of the diff --git a/src/io_meandata.F90.rej b/src/io_meandata.F90.rej new file mode 100644 index 000000000..3dc502956 --- /dev/null +++ b/src/io_meandata.F90.rej @@ -0,0 +1,11 @@ +diff a/src/io_meandata.F90 b/src/io_meandata.F90 (rejected hunks) +@@ -596,7 +596,8 @@ CASE ('icb ') + call def_stream(nod2D, myDim_nod2D, 'ibfwbv', 'basal iceberg melting', 'm/s', ibfwbv(:), 1, 'm', i_real4, partit, mesh) + call def_stream(nod2D, myDim_nod2D, 'ibfwl', 'lateral iceberg melting', 'm/s', ibfwl(:), 1, 'm', i_real4, partit, mesh) + call def_stream(nod2D, myDim_nod2D, 'ibfwe', 'iceberg erosion', 'm/s', ibfwe(:), 1, 'm', i_real4, partit, mesh) +- call def_stream(nod2D, myDim_nod2D, 'ibhf', 'heat flux from iceberg melting', 'm/s', ibhf(:), 1, 'm', i_real4, partit, mesh) ++ !call def_stream(nod2D, myDim_nod2D, 'ibhf', 'heat flux from iceberg melting', 'm/s', ibhf(:), 1, 'm', i_real4, partit, mesh) ++ call def_stream((/nl,nod2D/), (/nl,myDim_nod2D/), 'ibhf', 'heat flux from iceberg melting', 'm/s', ibhf_n(:,:), 1, 'm', i_real4, partit, mesh) + end if + !------------------------------------------ + diff --git a/src/oce_ale_tracer.F90 b/src/oce_ale_tracer.F90 index 3e64502ad..6704714ff 100644 --- a/src/oce_ale_tracer.F90 +++ b/src/oce_ale_tracer.F90 @@ -573,6 +573,7 @@ subroutine diff_ver_part_impl_ale(tr_num, dynamics, tracers, ice, partit, mesh) use bc_surface_interface use transit_bc_surface_interface use mod_ice + use iceberg_params implicit none integer , intent(in) , target :: tr_num type(t_dyn) , intent(inout), target :: dynamics @@ -988,6 +989,16 @@ subroutine diff_ver_part_impl_ale(tr_num, dynamics, tracers, ice, partit, mesh) end do end if + !_______________________________________________________________________ + ! case of activated shortwave penetration into the ocean, ad 3d contribution + 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)-ibhf_n(nz+1, n) * area(nz+1,n)/areasvol(nz,n)) * zinv / vcpw + end do + end if + !_______________________________________________________________________ ! The first row contains also the boundary condition from heatflux, ! freshwaterflux and relaxation terms diff --git a/src/oce_ale_tracer.F90.rej b/src/oce_ale_tracer.F90.rej new file mode 100644 index 000000000..f7caf0eb0 --- /dev/null +++ b/src/oce_ale_tracer.F90.rej @@ -0,0 +1,9 @@ +diff a/src/oce_ale_tracer.F90 b/src/oce_ale_tracer.F90 (rejected hunks) +@@ -452,6 +452,7 @@ subroutine diff_ver_part_impl_ale(tr_num, dynamics, tracers, partit, mesh) + use o_mixing_KPP_mod !for ghats _GO_ + use g_cvmix_kpp, only: kpp_nonlcltranspT, kpp_nonlcltranspS, kpp_oblmixc + use bc_surface_interface ++ use iceberg_params + implicit none + integer , intent(in) , target :: tr_num + type(t_dyn) , intent(inout), target :: dynamics From ca9b39a9dbdddf07a8088d1c1fc68a52b35de743 Mon Sep 17 00:00:00 2001 From: ackerlar Date: Fri, 27 Sep 2024 18:38:41 +0200 Subject: [PATCH 02/14] enable icebergs with cavities --- src/gen_modules_config.F90 | 6 +- src/icb_allocate.F90 | 4 +- src/icb_coupling.F90 | 55 ++++-- src/icb_coupling.F90.rej | 22 --- src/icb_dyn.F90 | 304 ++++++++++++++++++++++++++++++-- src/icb_elem.F90 | 43 +++-- src/icb_elem.F90.rej | 59 ------- src/icb_modules.F90 | 5 +- src/icb_modules.F90.rej | 72 -------- src/icb_step.F90 | 99 ++++++----- src/icb_step.F90.rej | 121 ------------- src/icb_thermo.F90 | 124 +++++++++----- src/icb_thermo.F90.rej | 324 ----------------------------------- src/ice_oce_coupling.F90 | 25 +-- src/ice_oce_coupling.F90.rej | 23 --- src/io_meandata.F90.rej | 11 -- src/oce_ale_tracer.F90.rej | 9 - 17 files changed, 522 insertions(+), 784 deletions(-) delete mode 100644 src/icb_coupling.F90.rej delete mode 100644 src/icb_elem.F90.rej delete mode 100644 src/icb_modules.F90.rej delete mode 100644 src/icb_step.F90.rej delete mode 100644 src/icb_thermo.F90.rej delete mode 100644 src/ice_oce_coupling.F90.rej delete mode 100644 src/io_meandata.F90.rej delete mode 100644 src/oce_ale_tracer.F90.rej diff --git a/src/gen_modules_config.F90 b/src/gen_modules_config.F90 index f4766764e..8916e845e 100755 --- a/src/gen_modules_config.F90 +++ b/src/gen_modules_config.F90 @@ -123,8 +123,9 @@ module g_config logical :: turn_off_fw=.false. logical :: use_icesheet_coupling=.false. logical :: lbalance_fw=.true. - logical :: lcell_saturation=.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 @@ -135,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, lbalance_fw, lcell_saturation, lmin_latent_hf, 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? diff --git a/src/icb_allocate.F90 b/src/icb_allocate.F90 index 8d235a415..202214397 100644 --- a/src/icb_allocate.F90 +++ b/src/icb_allocate.F90 @@ -94,9 +94,9 @@ subroutine allocate_icb(partit, mesh) allocate(fwb_flux_ib(ib_num)) allocate(fwbv_flux_ib(ib_num)) allocate(hfe_flux_ib(ib_num)) - allocate(hfl_flux_ib(ib_num)) + allocate(hfl_flux_ib(ib_num,mesh%nl)) allocate(hfb_flux_ib(ib_num)) - allocate(hfbv_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 diff --git a/src/icb_coupling.F90 b/src/icb_coupling.F90 index 56004644a..51671c8d1 100644 --- a/src/icb_coupling.F90 +++ b/src/icb_coupling.F90 @@ -14,7 +14,7 @@ subroutine reset_ib_fluxes() 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 @@ -27,7 +27,7 @@ subroutine prepare_icb2fesom(mesh, partit, ib,i_have_element,localelement,depth_ implicit none logical :: i_have_element - real, intent(in) :: depth_ib + real, intent(in) :: depth_ib, height_ib_single real :: lev_low, lev_up integer :: localelement integer :: iceberg_node @@ -44,7 +44,7 @@ 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 + if(i_have_element) then dz = 0.0 allocate(tot_area_nods_in_ib_elem(mesh%nl)) @@ -89,6 +89,7 @@ subroutine prepare_icb2fesom(mesh, partit, ib,i_have_element,localelement,depth_ 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(1) @@ -105,23 +106,34 @@ subroutine prepare_icb2fesom(mesh, partit, ib,i_have_element,localelement,depth_ 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)) .and. (abs(lev_up) 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 @@ -150,11 +162,24 @@ subroutine icb2fesom(mesh, partit, ice) #include "associate_part_ass.h" #include "associate_mesh_ass.h" + if (use_cavity) then +!$OMP PARALLEL DO + do n=1, myDim_nod2d+eDim_nod2D + 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 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 diff --git a/src/icb_coupling.F90.rej b/src/icb_coupling.F90.rej deleted file mode 100644 index fd3215347..000000000 --- a/src/icb_coupling.F90.rej +++ /dev/null @@ -1,22 +0,0 @@ -diff a/src/icb_coupling.F90 b/src/icb_coupling.F90 (rejected hunks) -@@ -90,16 +149,13 @@ type(t_partit), intent(inout), target :: partit - #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(:) - - do n=1, myDim_nod2d+eDim_nod2D -- if (.not.turn_off_hf) then -- net_heat_flux(n) = net_heat_flux(n) + ibhf(n) -- end if -+ !if (.not.turn_off_hf) then -+ ! heat_flux(n) = heat_flux(n) - ibhf(n) !* steps_per_ib_step -+ !end if - if (.not.turn_off_fw) then -- fresh_wa_flux(n) = fresh_wa_flux(n) + (ibfwb(n)+ibfwl(n)+ibfwe(n)+ibfwbv(n)) -+ water_flux(n) = water_flux(n) - (ibfwb(n)+ibfwl(n)+ibfwe(n)+ibfwbv(n)) !* steps_per_ib_step - end if - end do - !---wiso-code-begin diff --git a/src/icb_dyn.F90 b/src/icb_dyn.F90 index bc5f785e2..5184e2a2b 100644 --- a/src/icb_dyn.F90 +++ b/src/icb_dyn.F90 @@ -25,7 +25,7 @@ module iceberg_dynamics ! Thomas Rackow, 29.06.2010 !============================================================================== subroutine iceberg_dyn(mesh, partit, ice, dynamics, ib, new_u_ib, new_v_ib, u_ib, v_ib, lon,lat, depth_ib, & - height_ib, length_ib, width_ib, iceberg_elem, & + height_ib, length_ib, width_ib, iceberg_elem, & mass_ib, Ci, Ca, Co, Cda_skin, Cdo_skin, & rho_ice, rho_air, rho_h2o, P_sill, conc_sill, frozen_in, & file1, file2, P_ib, conci_ib, dt_ib, lastsubstep, & @@ -35,8 +35,8 @@ subroutine iceberg_dyn(mesh, partit, ice, dynamics, ib, new_u_ib, new_v_ib, u_ib use g_forcing_arrays !for u_wind, v_wind or u_wind_ib, v_wind_ib respectively use o_arrays, only: Tsurf_ib, Ssurf_ib use o_param !for dt - !use iceberg_params,only: l_melt, coriolis_scale !are icebergs allowed to melt? + integer :: ib_n_lvls, m integer, intent(IN) :: ib !current iceberg's index real, intent(OUT) :: new_u_ib, new_v_ib real, intent(IN) :: u_ib, v_ib @@ -61,21 +61,21 @@ subroutine iceberg_dyn(mesh, partit, ice, dynamics, ib, new_u_ib, new_v_ib, u_ib !LA 2023-03-07 real, dimension(:), pointer :: hi_ib3, conci_ib3, coriolis - real, dimension(3) :: uo_dz, vo_dz, uo_keel, vo_keel, T_dz,S_dz, T_keel,S_keel !hi_ib3, conci_ib3, + real, dimension(3) :: uo_keel, vo_keel, T_keel,S_keel, uo_dz, vo_dz, T_dz,S_dz!hi_ib3, conci_ib3, + real, dimension(:,:), allocatable :: arr_uo_dz, arr_vo_dz, arr_T_dz,arr_S_dz !hi_ib3, conci_ib3, real :: uo_ib, vo_ib, ua_ib, va_ib, ui_ib, vi_ib, hi_ib, uo_skin_ib, vo_skin_ib + real, dimension(:), allocatable :: arr_uo_ib, arr_vo_ib, arr_T_ave_ib, arr_S_ave_ib real :: Ao, Aa, Ai, Ad, fcoriolis real :: au_ib, av_ib real, dimension(2,2) :: SI_matrix real, dimension(2) :: SI_velo real :: u_ib_tmp, v_ib_tmp, normold, normnew, abs_omib, abs_omib_skin, ocean_drag - integer :: iter_ib + integer :: iter_ib, n, n2 real :: M_b, M_v, M_e, M_bv, sst_ib, sss_ib ! meltrates (basal, lateral, erosion, lateral 'basal'), temp. & salinity real :: T_ave_ib, S_ave_ib, T_keel_ib, S_keel_ib character, intent(IN) :: file3*80 real, intent(IN) :: rho_icb -! integer, dimension(3) :: tmp_arr - type(t_ice) , intent(inout), target :: ice type(t_mesh), intent(in) , target :: mesh type(t_partit), intent(inout), target :: partit @@ -85,19 +85,44 @@ subroutine iceberg_dyn(mesh, partit, ice, dynamics, ib, new_u_ib, new_v_ib, u_ib #include "associate_part_ass.h" #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))) +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))) +arr_uo_ib = 0.0 +arr_vo_ib = 0.0 +arr_T_ave_ib = 0.0 +arr_S_ave_ib = 0.0 + !OCEAN VELOCITIES: ! - (uo_ib, vo_ib) : integrated mean velocity at location of iceberg ! - (uo_skin_ib, vo_skin_ib) : velocity below the draft of the iceberg ! call iceberg_avvelo_ufkeel(uo_dz,vo_dz, uo_keel,vo_keel, depth_ib,iceberg_elem) call 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) + + !OCEANIC VELOCITY uo_ib, vo_ib call FEM_3eval(mesh, partit, uo_ib,vo_ib,lon,lat,uo_dz,vo_dz,iceberg_elem) + + call iceberg_levelwise_andkeel(mesh, partit, dynamics, arr_uo_dz,arr_vo_dz, uo_keel,vo_keel, arr_T_dz,arr_S_dz, T_keel,S_keel, depth_ib,iceberg_elem, ib, ib_n_lvls) + do n=1, ib_n_lvls + call FEM_3eval(mesh, partit, arr_uo_ib(n),arr_vo_ib(n),lon,lat,arr_uo_dz(:,n),arr_vo_dz(:,n),iceberg_elem) + call FEM_3eval(mesh, partit, arr_T_ave_ib(n),arr_S_ave_ib(n),lon,lat,arr_T_dz(:,n),arr_S_dz(:,n),iceberg_elem) + end do call FEM_3eval(mesh, partit, uo_skin_ib,vo_skin_ib,lon,lat,uo_keel,vo_keel,iceberg_elem) - !TEMPERATURE AND SALINITY: ! - T_ave_ib, S_ave_ib : Mean T & S (integrated) at location of iceberg ! - T_keel_ib, S_keel_ib : T & S below the draft of the iceberg (depth_ib) - call FEM_3eval(mesh, partit, T_ave_ib,S_ave_ib,lon,lat,T_dz,S_dz,iceberg_elem) call FEM_3eval(mesh, partit, T_keel_ib,S_keel_ib,lon,lat,T_keel,S_keel,iceberg_elem) @@ -133,14 +158,13 @@ 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, M_b, M_v, M_e, M_bv, & - u_ib,v_ib, uo_ib,vo_ib, ua_ib,va_ib, & + + 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, & uo_skin_ib, vo_skin_ib, T_keel_ib, S_keel_ib, depth_ib, height_ib, & - T_ave_ib, S_ave_ib, ib, rho_icb) + arr_T_ave_ib, arr_S_ave_ib, ib, rho_icb, uo_ib, vo_ib, ib_n_lvls, iceberg_elem, nlevels_nod2D(n2)) call iceberg_newdimensions(partit, ib, depth_ib,height_ib,length_ib,width_ib,M_b,M_v,M_e,M_bv, & rho_h2o, rho_icb, file3) @@ -580,6 +604,260 @@ subroutine compute_areas(Ao, Aa, Ai, Ad, depth_ib, & 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) + USE MOD_MESH + use o_param + use MOD_PARTIT + use MOD_DYN + + use o_arrays, only: Tclim_ib, Sclim_ib !, UV_ib, Z_3d_n_ib + + use g_clock + use g_forcing_arrays + use g_rotate_grid + + implicit none + + REAL, DIMENSION(:,:), allocatable, INTENT(OUT) :: uo_dz + REAL, DIMENSION(:,:), allocatable, INTENT(OUT) :: vo_dz + REAL, DIMENSION(3), INTENT(OUT) :: uo_keel + REAL, DIMENSION(3), INTENT(OUT) :: vo_keel + REAL, DIMENSION(:,:), allocatable, INTENT(OUT) :: T_dz + REAL, DIMENSION(:,:), allocatable, INTENT(OUT) :: S_dz + REAL, DIMENSION(3), INTENT(OUT) :: T_keel + REAL, DIMENSION(3), INTENT(OUT) :: S_keel + REAl, INTENT(IN) :: depth_ib + INTEGER :: ib_n_lvls_old + INTEGER, INTENT(IN) :: iceberg_elem, ib + INTEGER, INTENT(OUT) :: ib_n_lvls + INTEGER, dimension(3) :: arr_ib_n_lvls + REAL, dimension(:,:,:), pointer :: UV_ib + + real :: lev_up, lev_low + integer :: m, k, n2, n_up, n_low, cavity_count, max_node_level_count + ! depth over which is integrated (layer and sum) + real :: dz, ufkeel1, ufkeel2, Temkeel, Salkeel, ldepth_up, ldepth_low, dz_depth + +type(t_mesh), intent(in) , target :: mesh +type(t_dyn), intent(in) , target :: dynamics +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" + + UV_IB => dynamics%uv_ib(:,:,:) + cavity_count=0 + + do m=1,3 + if(m==1) then + max_node_level_count = nlevels_nod2D(elem2D_nodes(m,iceberg_elem)) + else + max_node_level_count = max(max_node_level_count, nlevels_nod2D(elem2D_nodes(m,iceberg_elem))) + end if + end do + + allocate(uo_dz(3,max_node_level_count)) + allocate(vo_dz(3,max_node_level_count)) + allocate(T_dz(3,max_node_level_count)) + allocate(S_dz(3,max_node_level_count)) + + ib_n_lvls_old = 0 + ib_n_lvls = 0 + arr_ib_n_lvls = 0 + + uo_dz = 0.0 + vo_dz = 0.0 + uo_keel = 0.0 + vo_keel = 0.0 + T_dz = 0.0 + S_dz = 0.0 + T_keel = 0.0 + S_keel = 0.0 + + !LOOP: over all nodes of the iceberg element + nodeloop: do m=1, 3 + !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) + 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) + 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 ) + 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 + + 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=abs(depth_ib) ) then !.AND. (abs(lev_up)<=abs(depth_ib)) ) then + if( abs(lev_up)= 0.0- 1.0e-07) ) then old_iceberg_elem=elem_containing_n2 - - if( (use_cavity) .and. (reject_elem(mesh, partit, old_iceberg_elem) )) then - left_mype=1.0 - !write(*,*) 'iceberg found in shelf region: left_mype = 1' - old_iceberg_elem=ibelem_tmp - end if + 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 end do @@ -576,7 +587,7 @@ SUBROUTINE locbafu_2D(mesh, partit, values, elem, coords) DO i=1, 2 TRANS(:,i) = local_cart_coords(:,i+1)-local_cart_coords(:,1) END DO - call matrix_inverse_2x2(TRANS, TRANS_inv, DET) + call matrix_inverse_2x2(TRANS, TRANS_inv, DET, elem, coords) vec=x_cart-local_cart_coords(:,1) stdel_coords = MATMUL(TRANS_inv, vec) @@ -835,13 +846,16 @@ end subroutine com_integer !END SUBROUTINE tides_distr !LA from oce_mesh_setup ofr iceberg coupling -subroutine matrix_inverse_2x2 (A, AINV, DET) +subroutine matrix_inverse_2x2 (A, AINV, DET, elem, coords) ! ! Coded by Sergey Danilov ! Reviewed by Qiang Wang !------------------------------------------------------------- implicit none + + integer :: elem + REAL, DIMENSION(2), INTENT(IN) :: coords real(kind=8), dimension(2,2), intent(IN) :: A real(kind=8), dimension(2,2), intent(OUT) :: AINV @@ -853,6 +867,7 @@ subroutine matrix_inverse_2x2 (A, AINV, DET) do j=1,2 write(*,*) (A(i,j),i=1,2) end do + write(*,*) " * elem: ", elem, ", coords: ", coords stop 'SINGULAR 2X2 MATRIX' else AINV(1,1) = A(2,2)/DET diff --git a/src/icb_elem.F90.rej b/src/icb_elem.F90.rej deleted file mode 100644 index a34f6a9f6..000000000 --- a/src/icb_elem.F90.rej +++ /dev/null @@ -1,59 +0,0 @@ -diff a/src/icb_elem.F90 b/src/icb_elem.F90 (rejected hunks) -@@ -376,9 +376,6 @@ 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 - -@@ -398,30 +395,21 @@ type(t_partit), intent(inout), target :: partit - - 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 -+ if( (use_cavity) .and. (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 - end if - call com_integer(partit, i_have_element,elem) !max 1 PE sends element here; - 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 - -@@ -461,14 +449,11 @@ do m=1, 3 - 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 -+ if( (use_cavity) .and. (reject_elem(mesh, partit, 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 -- - RETURN - end if - end do diff --git a/src/icb_modules.F90 b/src/icb_modules.F90 index 2b7aa088c..1162182db 100644 --- a/src/icb_modules.F90 +++ b/src/icb_modules.F90 @@ -95,7 +95,8 @@ module iceberg_params real,dimension(:), allocatable:: bvl_mean, lvlv_mean, lvle_mean, lvlb_mean !averaged volume losses !real,dimension(:), allocatable:: fw_flux_ib, hfb_flux_ib real,dimension(:), allocatable:: fwe_flux_ib, fwl_flux_ib, fwb_flux_ib, fwbv_flux_ib - real,dimension(:), allocatable:: hfe_flux_ib, hfl_flux_ib, hfb_flux_ib, hfbv_flux_ib, lhfb_flux_ib + real,dimension(:), allocatable:: hfe_flux_ib, hfb_flux_ib, lhfb_flux_ib + real,dimension(:,:), allocatable:: hfl_flux_ib, hfbv_flux_ib !===== FRESHWATER AND HEAT ARRAYS ON FESOM GRID ===== real,dimension(:), allocatable:: ibhf !icb heat flux into ocean @@ -129,7 +130,7 @@ logical function reject_elem(mesh, partit, elem) implicit none integer, intent(in) :: elem type(t_mesh), intent(in) , target :: mesh -type(t_partit), intent(inout), target :: partit +type(t_partit), intent(in), target :: partit #include "associate_mesh_def.h" #include "associate_part_def.h" #include "associate_mesh_ass.h" diff --git a/src/icb_modules.F90.rej b/src/icb_modules.F90.rej deleted file mode 100644 index 469d257ab..000000000 --- a/src/icb_modules.F90.rej +++ /dev/null @@ -1,72 +0,0 @@ -diff a/src/icb_modules.F90 b/src/icb_modules.F90 (rejected hunks) -@@ -1,4 +1,10 @@ - module iceberg_params -+use MOD_PARTIT -+USE MOD_MESH -+use g_config, only: use_cavity, use_cavityonelem -+use, intrinsic :: iso_fortran_env, only: real64 -+USE o_PARAM, only: WP -+ - implicit none - save - !integer,parameter :: ib_num ! realistic dataset comprising 6912 icebergs -@@ -111,35 +121,50 @@ save - 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 -+ logical function reject_elem(mesh, partit, elem) - implicit none - integer, intent(in) :: elem - 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 (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 -+ 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 -+ integer function coastal_nodes(mesh, partit, elem) - implicit none - integer, intent(in) :: elem - 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 (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 diff --git a/src/icb_step.F90 b/src/icb_step.F90 index ad3af496b..3517817ba 100644 --- a/src/icb_step.F90 +++ b/src/icb_step.F90 @@ -298,7 +298,7 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib_single,lengt != implicit none != - + logical :: reject_tmp integer, intent(in) :: ib, istep real, intent(inout) :: height_ib_single,length_ib_single,width_ib_single real, intent(inout) :: lon_deg,lat_deg @@ -361,7 +361,6 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib_single,lengt #include "associate_mesh_ass.h" -!write(*,*) "LA DEBUG: 1" mype =>partit%mype istep_end_synced = istep + steps_per_ib_step - 1 @@ -372,7 +371,6 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib_single,lengt lon_rad = lon_deg*rad lat_rad = lat_deg*rad -!write(*,*) "LA DEBUG: 2" if(volume_ib .le. smallestvol_icb) then melted(ib) = .true. @@ -382,7 +380,6 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib_single,lengt return end if -!write(*,*) "LA DEBUG: 3" if (firstcall) then if(mype==0) write(*,*) 'Preparing local_idx_of array...' @@ -392,7 +389,6 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib_single,lengt firstcall=.false. if(mype==0) write(*,*) 'Preparing local_idx_of done.' end if -!write(*,*) "LA DEBUG: 4" if (find_iceberg_elem) then lon_rad = lon_deg*rad @@ -409,12 +405,23 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib_single,lengt if(i_have_element) then i_have_element= mesh%elem2D_nodes(1,iceberg_elem) <= partit%myDim_nod2D !1 PE still .true. - if( (use_cavity) .and. (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 + + + + 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_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) @@ -425,10 +432,11 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib_single,lengt end if ! initialize the iceberg velocity -!write(*,*) "LA DEBUG: 4d" - call initialize_velo(mesh, partit, dynamics, i_have_element, ib, u_ib, v_ib, lon_rad, lat_rad, depth_ib, local_idx_of(iceberg_elem)) -!write(*,*) "LA DEBUG: 4e" - + if(local_idx_of(iceberg_elem) <= partit%myDim_elem2D ) then + call initialize_velo(mesh, partit, dynamics, i_have_element, ib, u_ib, v_ib, lon_rad, lat_rad, depth_ib, local_idx_of(iceberg_elem)) + else + write(*,*) " * skip initialize_velo" + end if !iceberg elem of ib is found find_iceberg_elem = .false. @@ -444,10 +452,8 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib_single,lengt endif endif -!write(*,*) "LA DEBUG: 4f" end if -!write(*,*) "LA DEBUG: 5" ! ================== START ICEBERG CALCULATION ==================== @@ -457,10 +463,8 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib_single,lengt !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 if( local_idx_of(iceberg_elem) > 0 ) then -write(*,*) "LA DEBUG: 5a" if( elem2D_nodes(1,local_idx_of(iceberg_elem)) <= partit%myDim_nod2D ) then -write(*,*) "LA DEBUG: 5b" i_have_element=.true. @@ -505,8 +509,7 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib_single,lengt ! 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. !if (mod(istep_end_synced,logfile_outfreq)==0) then -! if(ib==148) write(*,*) "LA DEBUG: 3 - elem ", iceberg_elem - write(*,*) 'iceberg ib ', ib, 'is grounded' + if (lverbose_icb) write(*,*) 'iceberg ib ', ib, 'is grounded' !end if else @@ -514,7 +517,6 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib_single,lengt t0=MPI_Wtime() -!write(*,*) "LA DEBUG: before_trajectory" call trajectory( lon_rad,lat_rad, u_ib,v_ib, new_u_ib,new_v_ib, & lon_deg,lat_deg,old_lon,old_lat, dt*REAL(steps_per_ib_step)) @@ -550,10 +552,16 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib_single,lengt !----------------------------- ! LA 2022-11-30 -!write(*,*) "LA DEBUG: before_saturation" - if(lcell_saturation) then -!write(*,*) "LA DEBUG: start_saturation" - area_ib_tot = length_ib_single*width_ib_single*scaling(ib) + if(cell_saturation > 0) then + if( lverbose_icb) then + write(*,*) " * checking for cell saturation - ib: ", ib, ", old_elem: ", old_element, ", new_elem: ", iceberg_elem + end if + select case(cell_saturation) !num of coastal points + case(1) + area_ib_tot = 0.0 + case(2) + area_ib_tot = length_ib_single*width_ib_single*scaling(ib) + end select !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(idx, area_ib_tot) !$OMP DO do idx = 1, size(elem_block) @@ -561,15 +569,17 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib_single,lengt area_ib_tot = area_ib_tot + length_ib(idx) * width_ib(idx) * scaling(idx) end if end do -!write(*,*) "LA DEBUG: end_loop_saturation" !$OMP END DO !$OMP END PARALLEL !----------------------------- - if ((area_ib_tot > elem_area_glob(iceberg_elem)) .and. (old_element.ne.0) .and. (left_mype == 0)) 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)) + if ((area_ib_tot > elem_area_glob(iceberg_elem)) .and. (old_element.ne.0) .and. (iceberg_elem.ne.old_element)) 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)) + + end if i_have_element = .true. left_mype = 0.0 lon_rad = old_lon @@ -580,7 +590,6 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib_single,lengt u_ib = 0. v_ib = 0. end if -!write(*,*) "LA DEBUG: after_cell_saturation" end if !########################################### @@ -594,9 +603,7 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib_single,lengt pe_block(ib)=mype volume_ib=height_ib_single*length_ib_single*width_ib_single -!write(*,*) "LA DEBUG: before_prepare" - call prepare_icb2fesom(mesh,partit,ib,i_have_element,local_idx_of(iceberg_elem),depth_ib) -!write(*,*) "LA DEBUG: after_prepare" + call prepare_icb2fesom(mesh,partit,ib,i_have_element,local_idx_of(iceberg_elem),depth_ib,height_ib_single) end if !processor has element? end if !... and first node belongs to processor? @@ -770,10 +777,15 @@ subroutine iceberg_step2(mesh, partit,arr, elem_from_block, ib, height_ib_single iceberg_elem = old_element u_ib = 0. v_ib = 0. - else if(lcell_saturation) then + else if(cell_saturation > 0) then if (mype==0) write(*,*) 'iceberg ',ib, ' changed PE or was very fast' elem_area_tmp = elem_area_glob(iceberg_elem) - area_ib_tot = length_ib_single * width_ib_single * scaling(ib) + select case(cell_saturation) !num of coastal points + case(1) + area_ib_tot = 0.0 + case(2) + area_ib_tot = length_ib_single*width_ib_single*scaling(ib) + end select !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(idx, area_ib_tot) !$OMP DO do idx = 1, size(elem_block_red) @@ -784,8 +796,8 @@ subroutine iceberg_step2(mesh, partit,arr, elem_from_block, ib, height_ib_single end do !$OMP END DO !$OMP END PARALLEL - if((area_ib_tot > elem_area_tmp) .and. (elem_area_tmp > 0.0)) then - if(mype==pe_block_red(ib)) then + if((area_ib_tot > elem_area_tmp) .and. (elem_area_tmp > 0.0) .and. (iceberg_elem.ne.old_element)) then + if(mype==pe_block_red(ib) .and. lverbose_icb) then write(*,*) " *******************************************************" write(*,*) " * iceberg changed PE and saturation" write(*,*) " * set iceberg ", ib, " back to elem ", old_element, " from elem ", iceberg_elem @@ -800,7 +812,7 @@ subroutine iceberg_step2(mesh, partit,arr, elem_from_block, ib, height_ib_single u_ib = 0. v_ib = 0. end if - else if(lcell_saturation) then + else if (mype==0) write(*,*) 'iceberg ',ib, ' changed PE or was very fast' end if end if @@ -1035,11 +1047,12 @@ subroutine parallel2coast(mesh, partit,u, v, lon,lat, elem) #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 = 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 = sum( mesh%bc_index_nod2D(mesh%elem2D_nodes(:,elem)) ) + fld_tmp = count( (mesh%bc_index_nod2D(mesh%elem2D_nodes(:,elem))==0.0) ) end if SELECT CASE ( fld_tmp ) !num of coastal points diff --git a/src/icb_step.F90.rej b/src/icb_step.F90.rej deleted file mode 100644 index f9c2a01af..000000000 --- a/src/icb_step.F90.rej +++ /dev/null @@ -1,121 +0,0 @@ -diff a/src/icb_step.F90 b/src/icb_step.F90 (rejected hunks) -@@ -376,22 +407,23 @@ type(t_dyn) , intent(inout), target :: dynamics - 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 -+!write(*,*) "LA DEBUG: 4a" - - if(i_have_element) then -+!write(*,*) "LA DEBUG: 4a1" - 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 -+!write(*,*) "LA DEBUG: 4a2" -+ if( (use_cavity) .and. (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 -+!write(*,*) "LA DEBUG: 4a3" - end if -+!write(*,*) "LA DEBUG: 4b" - call com_integer(partit, i_have_element,iceberg_elem) -+!write(*,*) "LA DEBUG: 4c" - - if(iceberg_elem .EQ. 0) then - write(*,*) 'IB ',ib,' rot. coords:', lon_deg, lat_deg !,lon_rad, lat_rad -@@ -983,15 +1023,14 @@ 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 - real, intent(in) :: lon, lat !radiant - integer, intent(in) :: elem - -+ integer :: fld_tmp - integer, dimension(3) :: n - integer :: node, m, i - real, dimension(2) :: velocity, velocity1, velocity2 -@@ -1004,12 +1043,13 @@ type(t_partit), intent(inout), target :: partit - #include "associate_part_ass.h" - #include "associate_mesh_ass.h" - --#ifdef use_cavity -- SELECT CASE ( coastal_nodes(mesh, 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 -+ if( use_cavity ) then -+ fld_tmp = coastal_nodes(mesh, partit, elem) -+ else -+ fld_tmp = sum( mesh%bc_index_nod2D(mesh%elem2D_nodes(:,elem)) ) -+ end if -+ -+ SELECT CASE ( fld_tmp ) !num of coastal points - CASE (0) !...coastal points: do nothing - return - -@@ -1020,14 +1060,18 @@ type(t_partit), intent(inout), target :: partit - 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. 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)==1 ) then -+ n(i) = node -+ exit -+ end if -+ end if - end do - - !write(*,*) 'one coastal node ', n(1) -@@ -1075,13 +1119,17 @@ type(t_partit), intent(inout), target :: partit - 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 -- n(i) = node -- i = i+1 -+ 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 -+ i = i+1 -+ end if -+ else -+ if( mesh%bc_index_nod2D(node)==1 ) then -+ n(i) = node -+ i = i+1 -+ end if - end if - end do - call projection(mesh,partit, velocity, n(1), n(2)) diff --git a/src/icb_thermo.F90 b/src/icb_thermo.F90 index 57cb85f63..7b9b38afe 100644 --- a/src/icb_thermo.F90 +++ b/src/icb_thermo.F90 @@ -16,14 +16,15 @@ ! use 3D information for T,S and velocities ! instead of SSTs; M_v depends on 'thermal driving') !============================================================================== -subroutine iceberg_meltrates(partit, M_b, M_v, M_e, M_bv, & - u_ib,v_ib, uo_ib,vo_ib, ua_ib,va_ib, & +subroutine 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, & uo_keel_ib, vo_keel_ib, T_keel_ib, S_keel_ib, depth_ib, height_ib, & - T_ave_ib, S_ave_ib, ib, rho_icb) + arr_T_ave_ib, arr_S_ave_ib, ib, rho_icb,uo_ib,vo_ib, ib_n_lvls, elem, n_lvls) use o_param use MOD_PARTIT + use MOD_MESH use g_clock use g_forcing_arrays use g_rotate_grid @@ -34,43 +35,84 @@ subroutine iceberg_meltrates(partit, M_b, M_v, M_e, M_bv, & ! LA: include latent heat 2023-04-04 real(kind=8),parameter :: L = 334000. ! [J/Kg] - real, intent(IN) :: u_ib,v_ib, uo_ib,vo_ib, ua_ib,va_ib !iceberg velo, (int.) ocean & atm velo + real, intent(IN) :: u_ib,v_ib, uo_ib,vo_ib,ua_ib,va_ib !iceberg velo, (int.) ocean & atm velo real, intent(IN) :: uo_keel_ib, vo_keel_ib !ocean velo at iceberg's draft real, intent(IN) :: sst_ib, length_ib, conci_ib, rho_icb !SST, length and sea ice conc. real, intent(IN) :: T_keel_ib, S_keel_ib, depth_ib, height_ib !T & S at depth 'depth_ib' - real, intent(IN) :: T_ave_ib, S_ave_ib !T & S averaged, i.e. at 'depth_ib/2' - integer, intent(IN) :: ib !iceberg ID + integer, intent(IN) :: ib, n_lvls, ib_n_lvls !iceberg ID + real, intent(IN), dimension(n_lvls) :: arr_T_ave_ib, arr_S_ave_ib !T & S averaged, i.e. at 'depth_ib/2' + real, intent(IN), dimension(n_lvls) :: arr_uo_ib, arr_vo_ib !iceberg velo, (int.) ocean & atm velo real, intent(OUT) :: M_b, M_v, M_e, M_bv !melt rates [m (ice) per s] real :: H_b, H_v, H_e, H_bv !melt rates [m (ice) per s] + real :: M_bv_dz, M_v_dz, dz_acc - + integer :: n, n2 + integer, intent(IN) :: elem + real :: lev_up, lev_low, dz real :: absamino, damping, sea_state, v_ibmino real :: tf, T_d !freezing temp. and 'thermal driving' +type(t_mesh), intent(in) , target :: mesh type(t_partit), intent(inout), target :: partit +!==================== MODULES & DECLARATIONS ==========================!= #include "associate_part_def.h" +#include "associate_mesh_def.h" #include "associate_part_ass.h" +#include "associate_mesh_ass.h" + !3-eq. formulation for bottom melting [m/s] v_ibmino = sqrt( (u_ib - uo_keel_ib)**2 + (v_ib - vo_keel_ib)**2 ) call iceberg_heat_water_fluxes_3eq(partit, ib, M_b, H_b, T_keel_ib,S_keel_ib,v_ibmino, depth_ib, tf) - hfb_flux_ib = H_b * length_ib*length_ib*scaling(ib) + hfb_flux_ib(ib) = H_b * length_ib*length_ib*scaling(ib) + + M_bv_dz = 0.0 + M_v_dz = 0.0 + dz_acc = 0.0 + + n2=elem2D_nodes(1,elem) + do n=1,ib_n_lvls !3-eq. formulation for lateral 'basal' melting [m/s] - v_ibmino = sqrt( (u_ib - uo_ib)**2 + (v_ib - vo_ib)**2 ) ! depth-average rel. velocity - call iceberg_heat_water_fluxes_3eq(partit, ib, M_bv, H_bv, T_ave_ib,S_ave_ib,v_ibmino, depth_ib/2.0, tf) - hfbv_flux_ib = H_bv * (2*length_ib*abs(depth_ib) + 2*length_ib*abs(depth_ib) ) * scaling(ib) + 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) + end if + + 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) + end if + dz_acc = dz_acc + dz + + 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) + 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). - T_d = T_ave_ib - tf - if(T_d < 0.) T_d = 0. - - !lateral melt (buoyant convection) - !M_v is a function of the 'thermal driving', NOT just sst! Cf. Neshyba and Josberger (1979) - M_v = 0.00762 * T_d + 0.00129 * T_d**2 - M_v = M_v/86400. - H_v = M_v * rho_icb * L - hfl_flux_ib = H_v * (2*length_ib*abs(depth_ib) + 2*length_ib*abs(depth_ib) ) * scaling(ib) + !'thermal driving', defined as the elevation of ambient water + !temperature above freezing point' (Neshyba and Josberger, 1979). + T_d = arr_T_ave_ib(n) - tf + if(T_d < 0.) T_d = 0. + + !lateral melt (buoyant convection) + !M_v = 0.00762 * sst_ib + 0.00129 * sst_ib**2 + !M_v = M_v/86400. + !M_v is a function of the 'thermal driving', NOT just sst! Cf. Neshyba and Josberger (1979) + M_v = 0.00762 * T_d + 0.00129 * T_d**2 + M_v = M_v/86400. + H_v = M_v * rho_icb * L + M_v_dz = M_v_dz + M_v*dz + 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) !wave erosion absamino = sqrt( (ua_ib - uo_ib)**2 + (va_ib - vo_ib)**2 ) @@ -79,7 +121,7 @@ subroutine iceberg_meltrates(partit, 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 = H_e * (length_ib*abs(height_ib) + length_ib*abs(height_ib) ) * scaling(ib) + hfe_flux_ib(ib) = H_e * (length_ib*abs(height_ib) + length_ib*abs(height_ib) ) * scaling(ib) !fwe_flux_ib = M_e end subroutine iceberg_meltrates @@ -236,23 +278,23 @@ subroutine iceberg_newdimensions(partit, ib, depth_ib,height_ib,length_ib,width_ !save in larger array vl_block((ib-1)*4+1 : ib*4)=arr - ! ----------------------- - ! LA: set iceberg heatflux at least to latent heat 2023-04-04 - ! Latent heat flux at base and sides also changes lines 475/476 - ! Lateral heat flux set to latent heat and basal heat flux set to zero - if( lmin_latent_hf ) then - lhfb_flux_ib(ib) = rho_icb*L*tvl*scaling(ib)/dt/REAL(steps_per_ib_step) - - hf_tot_tmp = hfb_flux_ib(ib)+hfbv_flux_ib(ib)+hfl_flux_ib(ib)+hfe_flux_ib(ib) - - if( (hf_tot_tmp >= 0.0) .and. (hf_tot_tmp < lhfb_flux_ib(ib))) then - hfe_flux_ib(ib)=-lhfb_flux_ib(ib) * abs(hfe_flux_ib(ib)/hf_tot_tmp) - hfl_flux_ib(ib)=-lhfb_flux_ib(ib) * abs(hfl_flux_ib(ib)/hf_tot_tmp) - hfb_flux_ib(ib)=-lhfb_flux_ib(ib) * abs(hfb_flux_ib(ib)/hf_tot_tmp) - hfbv_flux_ib(ib)=-lhfb_flux_ib(ib) * abs(hfbv_flux_ib(ib)/hf_tot_tmp) - end if - end if - ! ----------------------- + !! ----------------------- + !! LA: set iceberg heatflux at least to latent heat 2023-04-04 + !! Latent heat flux at base and sides also changes lines 475/476 + !! Lateral heat flux set to latent heat and basal heat flux set to zero + !if( lmin_latent_hf ) then + ! lhfb_flux_ib(ib) = rho_icb*L*tvl*scaling(ib)/dt/REAL(steps_per_ib_step) + + ! hf_tot_tmp = hfb_flux_ib(ib)+hfbv_flux_ib(ib)+hfl_flux_ib(ib)+hfe_flux_ib(ib) + + ! if( (hf_tot_tmp >= 0.0) .and. (abs(hf_tot_tmp) < abs(lhfb_flux_ib(ib)))) then + ! hfe_flux_ib(ib)=-lhfb_flux_ib(ib) * abs(hfe_flux_ib(ib)/hf_tot_tmp) + ! hfl_flux_ib(ib)=-lhfb_flux_ib(ib) * abs(hfl_flux_ib(ib)/hf_tot_tmp) + ! hfb_flux_ib(ib)=-lhfb_flux_ib(ib) * abs(hfb_flux_ib(ib)/hf_tot_tmp) + ! hfbv_flux_ib(ib)=-lhfb_flux_ib(ib) * abs(hfbv_flux_ib(ib)/hf_tot_tmp) + ! end if + !end if + !! ----------------------- end subroutine iceberg_newdimensions @@ -447,7 +489,7 @@ subroutine iceberg_heat_water_fluxes_3eq(partit, ib, M_b, H_b, T_ib,S_ib,v_rel, M_b = - (rhow / rhoi) * M_b ! [m (ice) per second], positive for melting? NOW positive for melting !LA avoid basal freezing for grounded icebergs - if(M_b.lt.0.) then + if(grounded(ib) .and. (M_b.lt.0.)) then M_b = 0.0 endif diff --git a/src/icb_thermo.F90.rej b/src/icb_thermo.F90.rej deleted file mode 100644 index 14d566b0e..000000000 --- a/src/icb_thermo.F90.rej +++ /dev/null @@ -1,324 +0,0 @@ -diff a/src/icb_thermo.F90 b/src/icb_thermo.F90 (rejected hunks) -@@ -16,45 +16,39 @@ - ! use 3D information for T,S and velocities - ! instead of SSTs; M_v depends on 'thermal driving') - !============================================================================== --subroutine iceberg_meltrates( M_b, M_v, M_e, M_bv, & -+subroutine iceberg_meltrates(partit, M_b, M_v, M_e, M_bv, & - u_ib,v_ib, uo_ib,vo_ib, ua_ib,va_ib, & - sst_ib, length_ib, conci_ib, & - uo_keel_ib, vo_keel_ib, T_keel_ib, S_keel_ib, depth_ib, & -- T_ave_ib, S_ave_ib, ib) -+ T_ave_ib, S_ave_ib, ib, rho_icb) - --! use o_mesh - use o_param --! use i_therm_param --! use i_param --! use MOD_ICE --! use i_arrays --! use MOD_PARTIT -- --! kh 18.03.21 not really used here --! use o_arrays -- -+ use MOD_PARTIT - use g_clock - use g_forcing_arrays - use g_rotate_grid - -- use iceberg_params, only: fwe_flux_ib, fwl_flux_ib, fwb_flux_ib, fwbv_flux_ib, heat_flux_ib -- -+ 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, height_ib, scaling - implicit none - -+ ! LA: include latent heat 2023-04-04 -+ real(kind=8),parameter :: L = 334000. ! [J/Kg] -+ - real, intent(IN) :: u_ib,v_ib, uo_ib,vo_ib, ua_ib,va_ib !iceberg velo, (int.) ocean & atm velo - real, intent(IN) :: uo_keel_ib, vo_keel_ib !ocean velo at iceberg's draft -- real, intent(IN) :: sst_ib, length_ib, conci_ib !SST, length and sea ice conc. -- real, intent(IN) :: T_keel_ib, S_keel_ib, depth_ib !T & S at depth 'depth_ib' -+ real, intent(IN) :: sst_ib, length_ib, conci_ib, rho_icb !SST, length and sea ice conc. -+ real, intent(IN) :: T_keel_ib, S_keel_ib, depth_ib !T & S at depth 'depth_ib' - real, intent(IN) :: T_ave_ib, S_ave_ib !T & S averaged, i.e. at 'depth_ib/2' - integer, intent(IN) :: ib !iceberg ID - real, intent(OUT) :: M_b, M_v, M_e, M_bv !melt rates [m (ice) per s] -+ real :: H_b, H_v, H_e, H_bv !melt rates [m (ice) per s] - - - real :: absamino, damping, sea_state, v_ibmino - real :: tf, T_d !freezing temp. and 'thermal driving' --!type(t_partit), intent(inout), target :: partit --!#include "associate_part_def.h" --!#include "associate_part_ass.h" -+type(t_partit), intent(inout), target :: partit -+#include "associate_part_def.h" -+#include "associate_part_ass.h" - - !bottom melt (basal turbulent melting rate) - !M_b = 0.58 * sqrt( (u_ib - uo_ib)**2 + (v_ib - vo_ib)**2 )**0.8 & -@@ -63,17 +57,18 @@ subroutine iceberg_meltrates( M_b, M_v, M_e, M_bv, & - - !3-eq. formulation for bottom melting [m/s] - v_ibmino = sqrt( (u_ib - uo_keel_ib)**2 + (v_ib - vo_keel_ib)**2 ) -- call iceberg_heat_water_fluxes_3eq(ib, M_b, T_keel_ib,S_keel_ib,v_ibmino, depth_ib, tf) -+ call iceberg_heat_water_fluxes_3eq(partit, ib, M_b, H_b, T_keel_ib,S_keel_ib,v_ibmino, depth_ib, tf) -+ hfb_flux_ib = H_b * length_ib*length_ib*scaling(ib) - - !3-eq. formulation for lateral 'basal' melting [m/s] - v_ibmino = sqrt( (u_ib - uo_ib)**2 + (v_ib - vo_ib)**2 ) ! depth-average rel. velocity -- call iceberg_heat_water_fluxes_3eq(ib, M_bv, T_ave_ib,S_ave_ib,v_ibmino, depth_ib/2.0, tf) -+ call iceberg_heat_water_fluxes_3eq(partit, ib, M_bv, H_bv, T_ave_ib,S_ave_ib,v_ibmino, depth_ib/2.0, tf) -+ hfbv_flux_ib = H_bv * (2*length_ib*abs(depth_ib) + 2*length_ib*abs(depth_ib) ) * scaling(ib) - - !'thermal driving', defined as the elevation of ambient water - !temperature above freezing point' (Neshyba and Josberger, 1979). - T_d = T_ave_ib - tf - if(T_d < 0.) T_d = 0. -- !write(*,*) 'thermal driving:',T_d,'; Tf:',tf,'T_ave:',T_ave_ib - - !lateral melt (buoyant convection) - !M_v = 0.00762 * sst_ib + 0.00129 * sst_ib**2 -@@ -81,6 +76,8 @@ subroutine iceberg_meltrates( M_b, M_v, M_e, M_bv, & - !M_v is a function of the 'thermal driving', NOT just sst! Cf. Neshyba and Josberger (1979) - M_v = 0.00762 * T_d + 0.00129 * T_d**2 - M_v = M_v/86400. -+ H_v = M_v * rho_icb * L -+ hfl_flux_ib = H_v * (2*length_ib*abs(depth_ib) + 2*length_ib*abs(depth_ib) ) * scaling(ib) - !fwl_flux_ib = M_v - - !wave erosion -@@ -106,21 +105,12 @@ end subroutine iceberg_meltrates - subroutine iceberg_newdimensions(partit, ib, depth_ib,height_ib,length_ib,width_ib,M_b,M_v,M_e,M_bv, & - rho_h2o, rho_icb, file_meltrates) - --! use o_mesh - use o_param !for step_per_day --! use i_therm_param --! use i_param --! use MOD_ICE --! use i_arrays - use MOD_PARTIT !for mype -- --! kh 18.03.21 not really used here --! use o_arrays -- - use g_clock - use g_forcing_arrays - use g_rotate_grid -- use iceberg_params, only: l_weeksmellor, ascii_out, icb_outfreq, vl_block, bvl_mean, lvlv_mean, lvle_mean, lvlb_mean, smallestvol_icb, fwb_flux_ib, fwe_flux_ib, fwbv_flux_ib, fwl_flux_ib, scaling, heat_flux_ib, lheat_flux_ib -+ use iceberg_params, only: l_weeksmellor, ascii_out, icb_outfreq, vl_block, bvl_mean, lvlv_mean, lvle_mean, lvlb_mean, smallestvol_icb, fwb_flux_ib, fwe_flux_ib, fwbv_flux_ib, fwl_flux_ib, scaling, hfb_flux_ib, hfbv_flux_ib, hfe_flux_ib, hfl_flux_ib, lhfb_flux_ib - use g_config, only: steps_per_ib_step - - implicit none -@@ -185,8 +176,6 @@ type(t_partit), intent(inout), target :: partit - volume_after=volume_before-tvl - - !calculating the new iceberg dimensions -- !depth_ib = (abs(depth_ib)-dh_b)*(-1.) -- !height_ib= abs(depth_ib) * rho_h2o/rho_icb - height_ib= height_ib - dh_b - depth_ib = -height_ib * rho_icb/rho_h2o - -@@ -196,15 +185,13 @@ type(t_partit), intent(inout), target :: partit - - !distribute dh_e equally between length and width - !as in code of michael schodlok, but not dh_v? -- !length_ib= length_ib - dh_v -dh_e/2. -- !width_ib = width_ib - dh_v -dh_e/2. - - volume_after=height_ib*length_ib*width_ib - - !iceberg smaller than critical value after melting? - if (volume_after .le. smallestvol_icb) then - volume_after=0.0 -- depth_ib = 0.0 -+ depth_ib = 0.0 - height_ib= 0.0 - length_ib= 0.0 - width_ib = 0.0 -@@ -297,25 +293,20 @@ end subroutine weeksmellor - !*************************************************************************************************************************** - !*************************************************************************************************************************** - --subroutine iceberg_heat_water_fluxes_3eq(ib, M_b, T_ib,S_ib,v_rel, depth_ib, t_freeze) -+subroutine iceberg_heat_water_fluxes_3eq(partit, ib, M_b, H_b, T_ib,S_ib,v_rel, depth_ib, t_freeze) - ! The three-equation model of ice-shelf ocean interaction (Hellmer et al., 1997) - ! Code derived from BRIOS subroutine iceshelf (which goes back to H.Hellmer's 2D ice shelf model code) - ! adjusted for use in FESOM by Ralph Timmermann, 16.02.2011 - ! adopted and modified for iceberg basal melting by Thomas Rackow, 11.06.2014 - !---------------------------------------------------------------- - -- !use o_mesh -- !use o_param -- !use o_arrays -- !use i_arrays -- !use MOD_PARTIT - use iceberg_params - use g_config - - implicit none - - integer, INTENT(IN) :: ib -- real(kind=8),INTENT(OUT) :: M_b, t_freeze -+ real(kind=8),INTENT(OUT) :: M_b, H_b, t_freeze - real(kind=8),INTENT(IN) :: T_ib, S_ib ! ocean temperature & salinity (at depth 'depth_ib') - real(kind=8),INTENT(IN) :: v_rel, depth_ib ! relative velocity iceberg-ocean (at depth 'depth_ib') - -@@ -351,22 +342,11 @@ subroutine iceberg_heat_water_fluxes_3eq(ib, M_b, T_ib,S_ib,v_rel, depth_ib, t_f - real(kind=8),parameter :: cpi = 152.5+7.122*(atk+tob) !Paterson:"The Physics of Glaciers" - - real(kind=8),parameter :: L = 334000. ! [J/Kg] -+type(t_partit), intent(inout), target :: partit -+!==================== MODULES & DECLARATIONS ==========================!= -+#include "associate_part_def.h" -+#include "associate_part_ass.h" - -- ! hemw = helium content of the glacial meltwater -- ! oomw = isotopic fractionation due to melting -- ! oofw = isotopic fractionation due to freezing -- ! hemw= 4.02*14. -- ! oomw= -30. -- ! oofw= -2.5 -- -- !n3=myDim_nod3d+eDim_nod3d -- -- !do n=1,myDim_nod2D+eDim_nod2D -- !if(cavity_flag_nod2d(n)==0) cycle -- !nk=nod3d_below_nod2d(1,n) -- !temp = tracer(nk,1) -- !sal = tracer(nk,2) -- !zice = coord_nod3d(3,nk) !(<0) - temp = T_ib - sal = S_ib - zice = depth_ib !(<0) -@@ -378,15 +358,7 @@ subroutine iceberg_heat_water_fluxes_3eq(ib, M_b, T_ib,S_ib,v_rel, depth_ib, t_f - ! Calculate or prescribe the turbulent heat and salt transfer coeff. GAT and GAS - ! velocity-dependent approach of Jenkins (1991) - -- !rt vt1 = 0.25*sqrt((u(i,j,N,lrhs)+u(i+1,j,N,lrhs))**2 -- !rt & +(v(i,j,N,lrhs)+v(i,j+1,N,lrhs))**2) -- ! if(vt1.eq.0.) vt1=0.001 -- !rt re = Hz_r(i,j,N)*ds/un !Reynolds number -- -- !vt1 = sqrt(uf(nk)*uf(nk)+uf(nk+n3)*uf(nk+n3)) ! relative velocity ice-ocean - vt1 = v_rel ! relative velocity iceberg-ocean (at depth 'depth_ib') -- --!rt RG44030 vt1 = max(vt1,0.001) - vt1 = max(vt1,0.005) ! RG44030 - - re = 10./un !vt1*re (=velocity times length scale over kinematic viscosity) is the Reynolds number -@@ -530,112 +503,6 @@ subroutine potit_ib(ib,salz,pt,pres,rfpres,tin) - return - 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 --! --!------------------------------------------------------------------------------------- --! --!subroutine potit(salz,pt,pres,rfpres,tin) --! ! Berechnet aus dem Salzgehalt[psu] (SALZ), der pot. Temperatur[oC] --! ! (PT) und dem Referenzdruck[dbar] (REFPRES) die in-situ Temperatur --! ! [oC] (TIN) bezogen auf den in-situ Druck[dbar] (PRES) mit Hilfe --! ! eines Iterationsverfahrens aus. --! --! integer iter --! real salz,pt,pres,rfpres,tin --! real epsi,tpmd,pt1,ptd,pttmpr --! --! data tpmd / 0.001 / --! --! epsi = 0. --! do iter=1,100 --! tin = pt+epsi --! pt1 = pttmpr(salz,tin,pres,rfpres) --! ptd = pt1-pt --! if(abs(ptd).lt.tpmd) return --! epsi = epsi-ptd --! enddo --! write(6,*) ' WARNING!' --! write(6,*) ' in-situ temperature calculation has not converged.' --! stop --! return --!end subroutine potit --! --!------------------------------------------------------------------------------------- --! --!real function pttmpr(salz,temp,pres,rfpres) --! ! Berechnet aus dem Salzgehalt/psu (SALZ), der in-situ Temperatur/degC --! ! (TEMP) und dem in-situ Druck/dbar (PRES) die potentielle Temperatur/ --! ! degC (PTTMPR) bezogen auf den Referenzdruck/dbar (RFPRES). Es wird --! ! ein Runge-Kutta Verfahren vierter Ordnung verwendet. --! ! Checkwert: PTTMPR = 36.89073 DegC --! ! fuer SALZ = 40.0 psu --! ! TEMP = 40.0 DegC --! ! PRES = 10000.000 dbar --! ! RFPRES = 0.000 dbar --! --! data ct2 ,ct3 /0.29289322 , 1.707106781/ --! data cq2a,cq2b /0.58578644 , 0.121320344/ --! data cq3a,cq3b /3.414213562, -4.121320344/ --! --! real salz,temp,pres,rfpres --! real p,t,dp,dt,q,ct2,ct3,cq2a,cq2b,cq3a,cq3b --! real adlprt --! --! p = pres --! t = temp --! dp = rfpres-pres --! dt = dp*adlprt(salz,t,p) --! t = t +0.5*dt --! q = dt --! p = p +0.5*dp --! dt = dp*adlprt(salz,t,p) --! t = t + ct2*(dt-q) --! q = cq2a*dt + cq2b*q --! dt = dp*adlprt(salz,t,p) --! t = t + ct3*(dt-q) --! q = cq3a*dt + cq3b*q --! p = rfpres --! dt = dp*adlprt(salz,t,p) --! --! pttmpr = t + (dt-q-q)/6.0 --! --!end function pttmpr --! --!------------------------------------------------------------------------------------- --! --!real function adlprt(salz,temp,pres) --! ! Berechnet aus dem Salzgehalt/psu (SALZ), der in-situ Temperatur/degC --! ! (TEMP) und dem in-situ Druck/dbar (PRES) den adiabatischen Temperatur- --! ! gradienten/(K Dbar^-1) ADLPRT. --! ! Checkwert: ADLPRT = 3.255976E-4 K dbar^-1 --! ! fuer SALZ = 40.0 psu --! ! TEMP = 40.0 DegC --! ! PRES = 10000.000 dbar --! --! real salz,temp,pres --! real s0,a0,a1,a2,a3,b0,b1,c0,c1,c2,c3,d0,d1,e0,e1,e2,ds --! --! data s0 /35.0/ --! data a0,a1,a2,a3 /3.5803E-5, 8.5258E-6, -6.8360E-8, 6.6228E-10/ --! data b0,b1 /1.8932E-6, -4.2393E-8/ --! data c0,c1,c2,c3 /1.8741E-8, -6.7795E-10, 8.7330E-12, -5.4481E-14/ --! data d0,d1 /-1.1351E-10, 2.7759E-12/ --! data e0,e1,e2 /-4.6206E-13, 1.8676E-14, -2.1687E-16/ --! --! ds = salz-s0 --! adlprt = ( ( (e2*temp + e1)*temp + e0 )*pres & --! + ( (d1*temp + d0)*ds & --! + ( (c3*temp + c2)*temp + c1 )*temp + c0 ) )*pres & --! + (b1*temp + b0)*ds + ( (a3*temp + a2)*temp + a1 )*temp + a0 --! --!END function adlprt --! --!---------------------------------------------------------------------------------------- --! --#endif -- - - ! LA from oce_dens_press for iceberg coupling - subroutine fcn_density(t,s,z,rho) diff --git a/src/ice_oce_coupling.F90 b/src/ice_oce_coupling.F90 index c33824352..270118675 100755 --- a/src/ice_oce_coupling.F90 +++ b/src/ice_oce_coupling.F90 @@ -385,6 +385,10 @@ subroutine oce_fluxes(ice, dynamics, tracers, partit, mesh) !$OMP END PARALLEL DO #endif + if (use_icebergs) then + call icb2fesom(mesh, partit, ice) + end if + !___________________________________________________________________________ ! add heat and fresh water flux from cavity if (use_cavity) then @@ -554,6 +558,16 @@ subroutine oce_fluxes(ice, dynamics, tracers, partit, mesh) !$OMP END PARALLEL DO end if + !___________________________________________________________________________ + if (use_icebergs) then + if (lbalance_fw .and. (.not. turn_off_fw)) then + flux = flux + (ibfwb + ibfwe + ibfwl + ibfwbv) !* steps_per_ib_step + end if + + call integrate_nod(ibfwb + ibfwe + ibfwl + ibfwbv, net, partit, mesh) + if (mype==0) write(*,*) " * total iceberg fw flux: ", net + end if + !___________________________________________________________________________ if (use_cavity) then ! with zstar we do not balance the freshwater flux under the cavity since its @@ -568,17 +582,6 @@ subroutine oce_fluxes(ice, dynamics, tracers, partit, mesh) where (ulevels_nod2d > 1) flux = -water_flux end if end if - - !___________________________________________________________________________ - if (use_icebergs .and. lbalance_fw) then - call icb2fesom(mesh, partit, ice) - if (.not.turn_off_fw) then - flux = flux + (ibfwb + ibfwe + ibfwl + ibfwbv) !* steps_per_ib_step - end if - - call integrate_nod(ibfwb + ibfwe + ibfwl + ibfwbv, net, partit, mesh) - if (mype==0) write(*,*) " * total iceberg fw flux: ", net - end if !___________________________________________________________________________ ! compute total global net freshwater flux into the ocean diff --git a/src/ice_oce_coupling.F90.rej b/src/ice_oce_coupling.F90.rej deleted file mode 100644 index 2d095512b..000000000 --- a/src/ice_oce_coupling.F90.rej +++ /dev/null @@ -1,23 +0,0 @@ -diff a/src/ice_oce_coupling.F90 b/src/ice_oce_coupling.F90 (rejected hunks) -@@ -563,9 +564,21 @@ subroutine oce_fluxes(ice, dynamics, tracers, partit, mesh) - end if - end if - -+ !___________________________________________________________________________ -+ if (use_icebergs .and. lbalance_fw) then -+ call icb2fesom(mesh, partit, ice) -+ if (.not.turn_off_fw) then -+ flux = flux + (ibfwb + ibfwe + ibfwl + ibfwbv) !* steps_per_ib_step -+ end if -+ -+ call integrate_nod(ibfwb + ibfwe + ibfwl + ibfwbv, net, partit, mesh) -+ if (mype==0) write(*,*) " * total iceberg fw flux: ", net -+ end if -+ - !___________________________________________________________________________ - ! compute total global net freshwater flux into the ocean - call integrate_nod(flux, net, partit, mesh) -+ - - !___________________________________________________________________________ - ! here the + sign must be used because we switched up the sign of the diff --git a/src/io_meandata.F90.rej b/src/io_meandata.F90.rej deleted file mode 100644 index 3dc502956..000000000 --- a/src/io_meandata.F90.rej +++ /dev/null @@ -1,11 +0,0 @@ -diff a/src/io_meandata.F90 b/src/io_meandata.F90 (rejected hunks) -@@ -596,7 +596,8 @@ CASE ('icb ') - call def_stream(nod2D, myDim_nod2D, 'ibfwbv', 'basal iceberg melting', 'm/s', ibfwbv(:), 1, 'm', i_real4, partit, mesh) - call def_stream(nod2D, myDim_nod2D, 'ibfwl', 'lateral iceberg melting', 'm/s', ibfwl(:), 1, 'm', i_real4, partit, mesh) - call def_stream(nod2D, myDim_nod2D, 'ibfwe', 'iceberg erosion', 'm/s', ibfwe(:), 1, 'm', i_real4, partit, mesh) -- call def_stream(nod2D, myDim_nod2D, 'ibhf', 'heat flux from iceberg melting', 'm/s', ibhf(:), 1, 'm', i_real4, partit, mesh) -+ !call def_stream(nod2D, myDim_nod2D, 'ibhf', 'heat flux from iceberg melting', 'm/s', ibhf(:), 1, 'm', i_real4, partit, mesh) -+ call def_stream((/nl,nod2D/), (/nl,myDim_nod2D/), 'ibhf', 'heat flux from iceberg melting', 'm/s', ibhf_n(:,:), 1, 'm', i_real4, partit, mesh) - end if - !------------------------------------------ - diff --git a/src/oce_ale_tracer.F90.rej b/src/oce_ale_tracer.F90.rej deleted file mode 100644 index f7caf0eb0..000000000 --- a/src/oce_ale_tracer.F90.rej +++ /dev/null @@ -1,9 +0,0 @@ -diff a/src/oce_ale_tracer.F90 b/src/oce_ale_tracer.F90 (rejected hunks) -@@ -452,6 +452,7 @@ subroutine diff_ver_part_impl_ale(tr_num, dynamics, tracers, partit, mesh) - use o_mixing_KPP_mod !for ghats _GO_ - use g_cvmix_kpp, only: kpp_nonlcltranspT, kpp_nonlcltranspS, kpp_oblmixc - use bc_surface_interface -+ use iceberg_params - implicit none - integer , intent(in) , target :: tr_num - type(t_dyn) , intent(inout), target :: dynamics From ec878cb1a6aba112f707539d41c88314ac5ad5e6 Mon Sep 17 00:00:00 2001 From: ackerlar Date: Fri, 27 Sep 2024 19:54:04 +0200 Subject: [PATCH 03/14] add module for icb_thermo and icb_coupling --- src/fesom_module.F90 | 1 + src/icb_coupling.F90 | 15 +++++++++++++++ src/icb_elem.F90 | 2 ++ src/icb_step.F90 | 2 ++ src/icb_thermo.F90 | 21 +++++++++++++++++++++ src/ice_oce_coupling.F90 | 1 + 6 files changed, 42 insertions(+) diff --git a/src/fesom_module.F90 b/src/fesom_module.F90 index 1fb782a0f..dffef9186 100755 --- a/src/fesom_module.F90 +++ b/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) diff --git a/src/icb_coupling.F90 b/src/icb_coupling.F90 index 51671c8d1..35b29ab9f 100644 --- a/src/icb_coupling.F90 +++ b/src/icb_coupling.F90 @@ -1,3 +1,17 @@ +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 @@ -190,3 +204,4 @@ subroutine icb2fesom(mesh, partit, ice) end if !---wiso-code-end end subroutine icb2fesom +end module diff --git a/src/icb_elem.F90 b/src/icb_elem.F90 index c79f2365b..e5de4e046 100644 --- a/src/icb_elem.F90 +++ b/src/icb_elem.F90 @@ -3,6 +3,8 @@ module iceberg_element USE MOD_MESH USE MOD_DYN use iceberg_params + use iceberg_thermodynamics + use iceberg_ocean_coupling ! use iceberg_dynamics ! use iceberg_step diff --git a/src/icb_step.F90 b/src/icb_step.F90 index 3517817ba..17ae8b31c 100644 --- a/src/icb_step.F90 +++ b/src/icb_step.F90 @@ -6,6 +6,8 @@ module iceberg_step use iceberg_params use iceberg_dynamics use iceberg_element + use iceberg_thermodynamics + use iceberg_ocean_coupling implicit none diff --git a/src/icb_thermo.F90 b/src/icb_thermo.F90 index 7b9b38afe..6aec637e5 100644 --- a/src/icb_thermo.F90 +++ b/src/icb_thermo.F90 @@ -1,3 +1,23 @@ +module iceberg_thermodynamics + USE MOD_MESH + use MOD_PARTIT + use MOD_ICE + USE MOD_DYN + use iceberg_params + !use iceberg_element + !use iceberg_step + use iceberg_ocean_coupling + + implicit none + + public :: iceberg_meltrates + public :: iceberg_newdimensions + public :: iceberg_heat_water_fluxes_3eq + public :: potit_ib + public :: fcn_density + + contains + !!============================================================================= ! calculates the empirical melt rates of the iceberg as in ! Martin: 'Parameterizing the fresh-water flux from land ice to ocean @@ -596,3 +616,4 @@ subroutine fcn_density(t,s,z,rho) + 4.8314e-4*s**2) rho = rhopot / (1.0 + 0.1*z/bulk) end subroutine fcn_density +end module diff --git a/src/ice_oce_coupling.F90 b/src/ice_oce_coupling.F90 index 270118675..55a259930 100755 --- a/src/ice_oce_coupling.F90 +++ b/src/ice_oce_coupling.F90 @@ -265,6 +265,7 @@ subroutine oce_fluxes(ice, dynamics, tracers, partit, mesh) init_flux_atm_ocn #endif use iceberg_params + use iceberg_ocean_coupling use cavity_interfaces !---fwf-code use g_clock From c8f287a9acbab0ed543b3f72f6b35f2fd6ebfc66 Mon Sep 17 00:00:00 2001 From: ackerlar Date: Fri, 27 Sep 2024 23:05:18 +0200 Subject: [PATCH 04/14] avoid unneccesary loop --- src/icb_step.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/icb_step.F90 b/src/icb_step.F90 index 17ae8b31c..f15b28c91 100644 --- a/src/icb_step.F90 +++ b/src/icb_step.F90 @@ -554,9 +554,10 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib_single,lengt !----------------------------- ! LA 2022-11-30 - if(cell_saturation > 0) then + if((cell_saturation>0) .and. (iceberg_elem.ne.old_element)) then if( lverbose_icb) then - write(*,*) " * checking for cell saturation - ib: ", ib, ", old_elem: ", old_element, ", new_elem: ", iceberg_elem + write(*,*) " * checking for cell saturation" + write(*,*) " * ib: ", ib, ", old_elem: ", old_element, ", new_elem: ", iceberg_elem end if select case(cell_saturation) !num of coastal points case(1) @@ -575,12 +576,11 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib_single,lengt !$OMP END PARALLEL !----------------------------- - if ((area_ib_tot > elem_area_glob(iceberg_elem)) .and. (old_element.ne.0) .and. (iceberg_elem.ne.old_element)) then ! .and. (left_mype == 0)) then + 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)) - end if i_have_element = .true. left_mype = 0.0 @@ -779,7 +779,7 @@ subroutine iceberg_step2(mesh, partit,arr, elem_from_block, ib, height_ib_single iceberg_elem = old_element u_ib = 0. v_ib = 0. - else if(cell_saturation > 0) then + else if((cell_saturation>0) .and. (iceberg_elem.ne.old_element)) then if (mype==0) write(*,*) 'iceberg ',ib, ' changed PE or was very fast' elem_area_tmp = elem_area_glob(iceberg_elem) select case(cell_saturation) !num of coastal points @@ -798,7 +798,7 @@ subroutine iceberg_step2(mesh, partit,arr, elem_from_block, ib, height_ib_single end do !$OMP END DO !$OMP END PARALLEL - if((area_ib_tot > elem_area_tmp) .and. (elem_area_tmp > 0.0) .and. (iceberg_elem.ne.old_element)) then + if((area_ib_tot > elem_area_tmp) .and. (elem_area_tmp > 0.0)) then if(mype==pe_block_red(ib) .and. lverbose_icb) then write(*,*) " *******************************************************" write(*,*) " * iceberg changed PE and saturation" From bfee944c50345de1537e8bb1bdd173ff1f994f93 Mon Sep 17 00:00:00 2001 From: ackerlar Date: Mon, 7 Oct 2024 14:48:02 +0200 Subject: [PATCH 05/14] fix segfault --- src/icb_step.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/icb_step.F90 b/src/icb_step.F90 index f15b28c91..74b38d293 100644 --- a/src/icb_step.F90 +++ b/src/icb_step.F90 @@ -798,7 +798,7 @@ subroutine iceberg_step2(mesh, partit,arr, elem_from_block, ib, height_ib_single end do !$OMP END DO !$OMP END PARALLEL - if((area_ib_tot > elem_area_tmp) .and. (elem_area_tmp > 0.0)) then + 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(*,*) " *******************************************************" write(*,*) " * iceberg changed PE and saturation" From 97d2441a3190eb580108e56e6210a7fc0fd32841 Mon Sep 17 00:00:00 2001 From: Pengyang Song Date: Tue, 8 Oct 2024 15:07:33 +0200 Subject: [PATCH 06/14] bug fix: use_cavity_fw2press does not work --- src/oce_ale.F90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/oce_ale.F90 b/src/oce_ale.F90 index 5f5f05470..8006f7632 100644 --- a/src/oce_ale.F90 +++ b/src/oce_ale.F90 @@ -1784,10 +1784,12 @@ subroutine compute_ssh_rhs_ale(dynamics, partit, mesh) !$OMP DO do n=1,myDim_nod2D nzmin = ulevels_nod2D(n) - if (ulevels_nod2D(n)>1 .and. use_cavity_fw2press ) then + if (ulevels_nod2D(n)>1) then ! use_cavity_fw2press=.true.: adds freshwater under the cavity thereby ! increasing the local pressure - ssh_rhs(n)=ssh_rhs(n)-alpha*water_flux(n)*areasvol(nzmin,n) + if (use_cavity_fw2press) then + ssh_rhs(n)=ssh_rhs(n)-alpha*water_flux(n)*areasvol(nzmin,n) + end if else ssh_rhs(n)=ssh_rhs(n)-alpha*water_flux(n)*areasvol(nzmin,n)+(1.0_WP-alpha)*ssh_rhs_old(n) end if From 3a5d6bae3d2ff621aa6f3b1f19494f6589a8d79f Mon Sep 17 00:00:00 2001 From: Jan Streffing Date: Tue, 8 Oct 2024 15:48:03 +0200 Subject: [PATCH 07/14] Revert "update mpif.h to use mpi" This reverts commit 0e5bfdc1b402347cc1b707cfb6396142e8aba447. --- src/MOD_PARTIT.F90 | 4 ++-- src/fortran_utils.F90 | 2 +- src/gen_modules_partitioning.F90 | 3 ++- src/ifs_interface/mpp_io.F90 | 3 ++- src/io_restart.F90 | 6 ++++-- src/temp/MOD_PARTIT.F90 | 4 ++-- 6 files changed, 13 insertions(+), 9 deletions(-) diff --git a/src/MOD_PARTIT.F90 b/src/MOD_PARTIT.F90 index 2e8330a4b..6603b8922 100644 --- a/src/MOD_PARTIT.F90 +++ b/src/MOD_PARTIT.F90 @@ -5,12 +5,12 @@ module MOD_PARTIT USE, intrinsic :: ISO_FORTRAN_ENV, only : int32 USE MOD_WRITE_BINARY_ARRAYS USE MOD_READ_BINARY_ARRAYS -USE mpi #if defined(_OPENMP) USE OMP_LIB #endif IMPLICIT NONE SAVE +include 'mpif.h' integer, parameter :: MAX_LAENDERECK=16 integer, parameter :: MAX_NEIGHBOR_PARTITIONS=32 @@ -217,4 +217,4 @@ subroutine READ_T_PARTIT(partit, unit, iostat, iomsg) read(unit, iostat=iostat, iomsg=iomsg) partit%pe_status end subroutine READ_T_PARTIT -end module MOD_PARTIT \ No newline at end of file +end module MOD_PARTIT diff --git a/src/fortran_utils.F90 b/src/fortran_utils.F90 index 86a93d226..35efdb742 100644 --- a/src/fortran_utils.F90 +++ b/src/fortran_utils.F90 @@ -1,6 +1,5 @@ ! synopsis: basic Fortran utilities, no MPI, dependencies only to INTRINSIC modules module fortran_utils - use mpi implicit none contains @@ -49,6 +48,7 @@ function mpirank_to_txt(mpicomm) result(txt) integer mype integer npes integer mpierr + include 'mpif.h' call MPI_Comm_Rank(mpicomm, mype, mpierr) call MPI_Comm_Size(mpicomm, npes, mpierr) diff --git a/src/gen_modules_partitioning.F90 b/src/gen_modules_partitioning.F90 index 2e701e427..04ecd94d3 100644 --- a/src/gen_modules_partitioning.F90 +++ b/src/gen_modules_partitioning.F90 @@ -113,7 +113,7 @@ subroutine par_ex(COMM, mype, abort) ! finalizes MPI #ifndef __oasis if (present(abort)) then if (mype==0) write(*,*) 'Run finished unexpectedly!' - call MPI_ABORT(MPI_COMM_WORLD, 1, error) + call MPI_ABORT(MPI_COMM_WORLD, 1 ) else ! TODO: this is where fesom standalone, ifsinterface etc get to !1. there no abort actually even when model calls abort, and barrier may hang @@ -580,3 +580,4 @@ subroutine status_check(partit) call par_ex(partit%MPI_COMM_FESOM, partit%mype, 1) endif end subroutine status_check + diff --git a/src/ifs_interface/mpp_io.F90 b/src/ifs_interface/mpp_io.F90 index d7838cd32..6c57995e3 100644 --- a/src/ifs_interface/mpp_io.F90 +++ b/src/ifs_interface/mpp_io.F90 @@ -6,7 +6,6 @@ !----------------------------------------------------- MODULE mpp_io - USE mpi #if defined(__MULTIO) USE iom, only : iom_enable_multio, iom_initialize, iom_init_server, iom_finalize #endif @@ -31,6 +30,7 @@ MODULE mpp_io SUBROUTINE mpp_io_init( iicomm, lio, irequired, iprovided, lmpi1 ) + INCLUDE "mpif.h" INTEGER, INTENT(INOUT) :: iicomm LOGICAL, INTENT(INOUT) :: lio INTEGER, INTENT(INOUT) :: irequired, iprovided @@ -126,6 +126,7 @@ SUBROUTINE mpp_io_init_2( iicomm ) INTEGER :: icode, ierr, icolor, iicommx, iicommm, iicommo INTEGER :: ji,inum LOGICAL :: lcompp + INCLUDE "mpif.h" ! Construct multio server communicator diff --git a/src/io_restart.F90 b/src/io_restart.F90 index fd1121278..9077e27af 100644 --- a/src/io_restart.F90 +++ b/src/io_restart.F90 @@ -12,7 +12,6 @@ MODULE io_RESTART use MOD_PARTIT use MOD_PARSUP use fortran_utils - use mpi #if defined(__recom) use recom_glovar use recom_config @@ -772,6 +771,7 @@ subroutine read_all_raw_restarts(mpicomm, mype) integer fileunit integer status integer mpierr + include 'mpif.h' if(mype == RAW_RESTART_METADATA_RANK) then ! read metadata info for the raw restart @@ -860,6 +860,7 @@ subroutine finalize_restart() ! !_______________________________________________________________________________ subroutine read_restart(path, filegroup, mpicomm, mype) + include 'mpif.h' character(len=*), intent(in) :: path type(restart_file_group), intent(inout) :: filegroup integer, intent(in) :: mpicomm @@ -1003,6 +1004,7 @@ function is_due(unit, frequency, istep) result(d) ! integer mype ! integer npes ! integer mpierr +! include 'mpif.h' ! ! call MPI_Comm_Rank(mpicomm, mype, mpierr) ! call MPI_Comm_Size(mpicomm, npes, mpierr) @@ -1010,4 +1012,4 @@ function is_due(unit, frequency, istep) result(d) ! end function !!PS --> move this function also to fortran_utils.F90 -end module \ No newline at end of file +end module diff --git a/src/temp/MOD_PARTIT.F90 b/src/temp/MOD_PARTIT.F90 index 818b46541..bd3b7dec2 100644 --- a/src/temp/MOD_PARTIT.F90 +++ b/src/temp/MOD_PARTIT.F90 @@ -5,9 +5,9 @@ module MOD_PARTIT USE, intrinsic :: ISO_FORTRAN_ENV USE MOD_WRITE_BINARY_ARRAYS USE MOD_READ_BINARY_ARRAYS -USE mpi IMPLICIT NONE SAVE +include 'mpif.h' integer, parameter :: MAX_LAENDERECK=16 integer, parameter :: MAX_NEIGHBOR_PARTITIONS=32 @@ -186,4 +186,4 @@ subroutine READ_T_PARTIT(partit, unit, iostat, iomsg) read(unit, iostat=iostat, iomsg=iomsg) partit%pe_status end subroutine READ_T_PARTIT -end module MOD_PARTIT \ No newline at end of file +end module MOD_PARTIT From f6e04d6b4ec47ae149fd62f1d1e296fe8031619a Mon Sep 17 00:00:00 2001 From: ackerlar Date: Tue, 8 Oct 2024 16:24:36 +0200 Subject: [PATCH 08/14] fix segfault, really... --- src/icb_step.F90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/icb_step.F90 b/src/icb_step.F90 index 74b38d293..4a133ccfe 100644 --- a/src/icb_step.F90 +++ b/src/icb_step.F90 @@ -464,8 +464,7 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib_single,lengt 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 -if( local_idx_of(iceberg_elem) > 0 ) then - +if((local_idx_of(iceberg_elem)>0) .and. (local_idx_of(iceberg_elem)<=partit%myDim_elem2D) ) then if( elem2D_nodes(1,local_idx_of(iceberg_elem)) <= partit%myDim_nod2D ) then i_have_element=.true. From d3aca8b029dbea5bf0cb162df8b96fe56dbf0842 Mon Sep 17 00:00:00 2001 From: Sebastian Beyer Date: Mon, 16 Sep 2024 13:31:21 +0200 Subject: [PATCH 09/14] update mpif.h to use mpi --- src/MOD_PARTIT.F90 | 4 ++-- src/fortran_utils.F90 | 2 +- src/gen_modules_partitioning.F90 | 3 +-- src/ifs_interface/mpp_io.F90 | 3 +-- src/io_restart.F90 | 6 ++---- src/temp/MOD_PARTIT.F90 | 4 ++-- 6 files changed, 9 insertions(+), 13 deletions(-) diff --git a/src/MOD_PARTIT.F90 b/src/MOD_PARTIT.F90 index 6603b8922..2e8330a4b 100644 --- a/src/MOD_PARTIT.F90 +++ b/src/MOD_PARTIT.F90 @@ -5,12 +5,12 @@ module MOD_PARTIT USE, intrinsic :: ISO_FORTRAN_ENV, only : int32 USE MOD_WRITE_BINARY_ARRAYS USE MOD_READ_BINARY_ARRAYS +USE mpi #if defined(_OPENMP) USE OMP_LIB #endif IMPLICIT NONE SAVE -include 'mpif.h' integer, parameter :: MAX_LAENDERECK=16 integer, parameter :: MAX_NEIGHBOR_PARTITIONS=32 @@ -217,4 +217,4 @@ subroutine READ_T_PARTIT(partit, unit, iostat, iomsg) read(unit, iostat=iostat, iomsg=iomsg) partit%pe_status end subroutine READ_T_PARTIT -end module MOD_PARTIT +end module MOD_PARTIT \ No newline at end of file diff --git a/src/fortran_utils.F90 b/src/fortran_utils.F90 index 35efdb742..86a93d226 100644 --- a/src/fortran_utils.F90 +++ b/src/fortran_utils.F90 @@ -1,5 +1,6 @@ ! synopsis: basic Fortran utilities, no MPI, dependencies only to INTRINSIC modules module fortran_utils + use mpi implicit none contains @@ -48,7 +49,6 @@ function mpirank_to_txt(mpicomm) result(txt) integer mype integer npes integer mpierr - include 'mpif.h' call MPI_Comm_Rank(mpicomm, mype, mpierr) call MPI_Comm_Size(mpicomm, npes, mpierr) diff --git a/src/gen_modules_partitioning.F90 b/src/gen_modules_partitioning.F90 index 04ecd94d3..2e701e427 100644 --- a/src/gen_modules_partitioning.F90 +++ b/src/gen_modules_partitioning.F90 @@ -113,7 +113,7 @@ subroutine par_ex(COMM, mype, abort) ! finalizes MPI #ifndef __oasis if (present(abort)) then if (mype==0) write(*,*) 'Run finished unexpectedly!' - call MPI_ABORT(MPI_COMM_WORLD, 1 ) + call MPI_ABORT(MPI_COMM_WORLD, 1, error) else ! TODO: this is where fesom standalone, ifsinterface etc get to !1. there no abort actually even when model calls abort, and barrier may hang @@ -580,4 +580,3 @@ subroutine status_check(partit) call par_ex(partit%MPI_COMM_FESOM, partit%mype, 1) endif end subroutine status_check - diff --git a/src/ifs_interface/mpp_io.F90 b/src/ifs_interface/mpp_io.F90 index 6c57995e3..d7838cd32 100644 --- a/src/ifs_interface/mpp_io.F90 +++ b/src/ifs_interface/mpp_io.F90 @@ -6,6 +6,7 @@ !----------------------------------------------------- MODULE mpp_io + USE mpi #if defined(__MULTIO) USE iom, only : iom_enable_multio, iom_initialize, iom_init_server, iom_finalize #endif @@ -30,7 +31,6 @@ MODULE mpp_io SUBROUTINE mpp_io_init( iicomm, lio, irequired, iprovided, lmpi1 ) - INCLUDE "mpif.h" INTEGER, INTENT(INOUT) :: iicomm LOGICAL, INTENT(INOUT) :: lio INTEGER, INTENT(INOUT) :: irequired, iprovided @@ -126,7 +126,6 @@ SUBROUTINE mpp_io_init_2( iicomm ) INTEGER :: icode, ierr, icolor, iicommx, iicommm, iicommo INTEGER :: ji,inum LOGICAL :: lcompp - INCLUDE "mpif.h" ! Construct multio server communicator diff --git a/src/io_restart.F90 b/src/io_restart.F90 index 9077e27af..fd1121278 100644 --- a/src/io_restart.F90 +++ b/src/io_restart.F90 @@ -12,6 +12,7 @@ MODULE io_RESTART use MOD_PARTIT use MOD_PARSUP use fortran_utils + use mpi #if defined(__recom) use recom_glovar use recom_config @@ -771,7 +772,6 @@ subroutine read_all_raw_restarts(mpicomm, mype) integer fileunit integer status integer mpierr - include 'mpif.h' if(mype == RAW_RESTART_METADATA_RANK) then ! read metadata info for the raw restart @@ -860,7 +860,6 @@ subroutine finalize_restart() ! !_______________________________________________________________________________ subroutine read_restart(path, filegroup, mpicomm, mype) - include 'mpif.h' character(len=*), intent(in) :: path type(restart_file_group), intent(inout) :: filegroup integer, intent(in) :: mpicomm @@ -1004,7 +1003,6 @@ function is_due(unit, frequency, istep) result(d) ! integer mype ! integer npes ! integer mpierr -! include 'mpif.h' ! ! call MPI_Comm_Rank(mpicomm, mype, mpierr) ! call MPI_Comm_Size(mpicomm, npes, mpierr) @@ -1012,4 +1010,4 @@ function is_due(unit, frequency, istep) result(d) ! end function !!PS --> move this function also to fortran_utils.F90 -end module +end module \ No newline at end of file diff --git a/src/temp/MOD_PARTIT.F90 b/src/temp/MOD_PARTIT.F90 index bd3b7dec2..818b46541 100644 --- a/src/temp/MOD_PARTIT.F90 +++ b/src/temp/MOD_PARTIT.F90 @@ -5,9 +5,9 @@ module MOD_PARTIT USE, intrinsic :: ISO_FORTRAN_ENV USE MOD_WRITE_BINARY_ARRAYS USE MOD_READ_BINARY_ARRAYS +USE mpi IMPLICIT NONE SAVE -include 'mpif.h' integer, parameter :: MAX_LAENDERECK=16 integer, parameter :: MAX_NEIGHBOR_PARTITIONS=32 @@ -186,4 +186,4 @@ subroutine READ_T_PARTIT(partit, unit, iostat, iomsg) read(unit, iostat=iostat, iomsg=iomsg) partit%pe_status end subroutine READ_T_PARTIT -end module MOD_PARTIT +end module MOD_PARTIT \ No newline at end of file From ed87d983366b8caa8f9c81ca6bec72ac8750ccc3 Mon Sep 17 00:00:00 2001 From: Sebastian Beyer Date: Tue, 8 Oct 2024 21:23:40 +0200 Subject: [PATCH 10/14] add use mpi to cpl_driver --- src/cpl_driver.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/cpl_driver.F90 b/src/cpl_driver.F90 index 0c2ad2bc1..8aa733bdc 100644 --- a/src/cpl_driver.F90 +++ b/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 ! @@ -1005,4 +1006,4 @@ end module cpl_driver #else module cpl_driver end module cpl_driver -#endif +#endif \ No newline at end of file From 74457684a3064960eabb4a49149f8187f3747a5b Mon Sep 17 00:00:00 2001 From: ackerlar Date: Wed, 9 Oct 2024 13:19:06 +0200 Subject: [PATCH 11/14] small fix --- src/icb_step.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/icb_step.F90 b/src/icb_step.F90 index 4a133ccfe..a7da3baa8 100644 --- a/src/icb_step.F90 +++ b/src/icb_step.F90 @@ -464,7 +464,7 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib_single,lengt 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 -if((local_idx_of(iceberg_elem)>0) .and. (local_idx_of(iceberg_elem)<=partit%myDim_elem2D) ) then +if((local_idx_of(iceberg_elem)>0) .and. (local_idx_of(iceberg_elem)<=partit%myDim_elem2D+partit%eDim_elem2D) ) then if( elem2D_nodes(1,local_idx_of(iceberg_elem)) <= partit%myDim_nod2D ) then i_have_element=.true. From 9be15d9f9d01cfc0e1e6338005d919a9ee4f876a Mon Sep 17 00:00:00 2001 From: ackerlar Date: Thu, 10 Oct 2024 13:20:55 +0200 Subject: [PATCH 12/14] change commRank type to integer --- src/cpl_driver.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/cpl_driver.F90 b/src/cpl_driver.F90 index 8aa733bdc..9e62800de 100644 --- a/src/cpl_driver.F90 +++ b/src/cpl_driver.F90 @@ -57,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 @@ -774,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 ??? @@ -1006,4 +1006,4 @@ end module cpl_driver #else module cpl_driver end module cpl_driver -#endif \ No newline at end of file +#endif From fcf69bbffe8ab3623ff824d64e379f0d14583807 Mon Sep 17 00:00:00 2001 From: ackerlar Date: Sun, 13 Oct 2024 10:50:27 +0200 Subject: [PATCH 13/14] fix silly misstake... --- src/icb_coupling.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/icb_coupling.F90 b/src/icb_coupling.F90 index 35b29ab9f..7f8f6b7e1 100644 --- a/src/icb_coupling.F90 +++ b/src/icb_coupling.F90 @@ -121,7 +121,7 @@ subroutine prepare_icb2fesom(mesh, partit, ib,i_have_element,localelement,depth_ end if dz = abs( lev_low - lev_up ) if( (abs(lev_low)>=abs(depth_ib)) .and. (abs(lev_up) Date: Sun, 13 Oct 2024 10:50:45 +0200 Subject: [PATCH 14/14] write out 3D icb heat flux --- src/io_meandata.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/io_meandata.F90 b/src/io_meandata.F90 index afc58f87f..5b2fc8876 100644 --- a/src/io_meandata.F90 +++ b/src/io_meandata.F90 @@ -1049,7 +1049,7 @@ subroutine ini_mean_io(ice, dynamics, tracers, partit, mesh) call def_stream(nod2D, myDim_nod2D, 'ibfwbv', 'basal iceberg melting', 'm/s', ibfwbv(:), 1, 'm', i_real4, partit, mesh) call def_stream(nod2D, myDim_nod2D, 'ibfwl', 'lateral iceberg melting', 'm/s', ibfwl(:), 1, 'm', i_real4, partit, mesh) call def_stream(nod2D, myDim_nod2D, 'ibfwe', 'iceberg erosion', 'm/s', ibfwe(:), 1, 'm', i_real4, partit, mesh) - call def_stream(nod2D, myDim_nod2D, 'ibhf', 'heat flux from iceberg melting', 'W/m2', ibhf(:), 1, 'm', i_real4, partit, mesh) + call def_stream((/nl,nod2D/), (/nl,myDim_nod2D/), 'ibhf', 'heat flux from iceberg melting', 'm/s', ibhf_n(:,:), 1, 'm', i_real4, partit, mesh) end if !------------------------------------------ !_______________________________________________________________________________