Skip to content

Commit

Permalink
3D heat flux, fix cell saturation
Browse files Browse the repository at this point in the history
  • Loading branch information
@ackerlar
ackerlar committed Sep 20, 2024
1 parent 59f9b2b commit 361b9b9
Show file tree
Hide file tree
Showing 19 changed files with 1,055 additions and 374 deletions.
12 changes: 1 addition & 11 deletions src/fesom_module.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
! --------------

Expand Down Expand Up @@ -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
Expand Down
5 changes: 4 additions & 1 deletion src/gen_modules_config.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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?
Expand Down
43 changes: 33 additions & 10 deletions src/icb_allocate.F90
View file
Original file line number Diff line number Diff line change
@@ -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
Expand Down Expand Up @@ -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
Expand Down
105 changes: 79 additions & 26 deletions src/icb_coupling.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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


Expand All @@ -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"
Expand All @@ -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)<abs(depth_ib)) ) then
dz = abs( lev_up - depth_ib )
end if

if( depth_ib==0.0 ) then
ibhf_n(j,iceberg_node) = ibhf_n(j,iceberg_node) &
- ((hfbv_flux_ib(ib)+hfl_flux_ib(ib)) &
+ hfe_flux_ib(ib)) / tot_area_nods_in_ib_elem(j)
else
ibhf_n(j,iceberg_node) = ibhf_n(j,iceberg_node) &
- ((hfbv_flux_ib(ib)+hfl_flux_ib(ib)) * (dz / abs(depth_ib)) &
+ hfe_flux_ib(ib) * (dz / abs(height_ib(ib)))) &
/ tot_area_nods_in_ib_elem(j)
end if
end do
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))
ibhf_n(1,iceberg_node) = ibhf_n(1,iceberg_node) - hfe_flux_ib(ib) * ((abs(height_ib(ib))-abs(depth_ib))/abs(height_ib(ib))) / tot_area_nods_in_ib_elem(1)
end if
end do
end if
Expand Down Expand Up @@ -90,16 +149,10 @@ subroutine icb2fesom(mesh, partit, ice)
#include "associate_mesh_def.h"
#include "associate_part_ass.h"
#include "associate_mesh_ass.h"

fresh_wa_flux => ice%flx_fw(:)
net_heat_flux => ice%flx_h(:)

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
Expand Down
22 changes: 22 additions & 0 deletions src/icb_coupling.F90.rej
View file
Original file line number Diff line number Diff line change
@@ -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
30 changes: 19 additions & 11 deletions src/icb_dyn.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down
Loading

0 comments on commit 361b9b9

Please sign in to comment.