Skip to content

Commit

Permalink
Store iceberg properties in forces
Browse files Browse the repository at this point in the history 
  Added code storing aggregate iceberg properties in the mech_forcing type from
the Ice_ocean_boundary type.  Also reduced the index range over which
rigidity_ice_[uv] is set due to sea-ice to the range that is actually used. All
answers are bitwise identical.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed May 1, 2018
1 parent 6570af3 commit e7ffe07
Show file tree
Hide file tree
Showing 2 changed files with 39 additions and 32 deletions.
31 changes: 20 additions & 11 deletions config_src/coupled_driver/MOM_surface_forcing.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -307,6 +307,12 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, &
fluxes%dt_buoy_accum = 0.0
endif ! endif for allocation and initialization


if (((associated(IOB%ustar_berg) .and. (.not. associated(fluxes%ustar_berg))) &
.or. (associated(IOB%area_berg) .and. (.not. associated(fluxes%area_berg)))) &
.or. (associated(IOB%mass_berg) .and. (.not. associated(fluxes%mass_berg)))) &
call allocate_forcing_type(G, fluxes, iceberg=.true.)

if ((.not.coupler_type_initialized(fluxes%tr_fluxes)) .and. &
coupler_type_initialized(IOB%fluxes)) &
call coupler_type_spawn(IOB%fluxes, fluxes%tr_fluxes, &
Expand Down Expand Up @@ -417,11 +423,6 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, &
if (associated(IOB%calving)) &
fluxes%frunoff(i,j) = IOB%calving(i-i0,j-j0) * G%mask2dT(i,j)

if (((associated(IOB%ustar_berg) .and. (.not. associated(fluxes%ustar_berg))) &
.or. (associated(IOB%area_berg) .and. (.not. associated(fluxes%area_berg)))) &
.or. (associated(IOB%mass_berg) .and. (.not. associated(fluxes%mass_berg)))) &
call allocate_forcing_type(G, fluxes, iceberg=.true.)

if (associated(IOB%ustar_berg)) &
fluxes%ustar_berg(i,j) = IOB%ustar_berg(i-i0,j-j0) * G%mask2dT(i,j)

Expand Down Expand Up @@ -628,6 +629,11 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, CS)
forces%initialized = .true.
endif

if ( (associated(IOB%area_berg) .and. (.not. associated(forces%area_berg))) .or. &
(associated(IOB%mass_berg) .and. (.not. associated(forces%mass_berg))) ) &
call allocate_mech_forcing(G, forces, iceberg=.true.)


! applied surface pressure from atmosphere and cryosphere
if (associated(IOB%p)) then
if (CS%max_p_surf >= 0.0) then
Expand Down Expand Up @@ -657,6 +663,12 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, CS)

! obtain fluxes from IOB; note the staggering of indices
do j=js,je ; do i=is,ie
if (associated(IOB%area_berg)) &
forces%area_berg(i,j) = IOB%area_berg(i-i0,j-j0) * G%mask2dT(i,j)

if (associated(IOB%mass_berg)) &
forces%mass_berg(i,j) = IOB%mass_berg(i-i0,j-j0) * G%mask2dT(i,j)

if (wind_stagger == BGRID_NE) then
if (associated(IOB%u_flux)) taux_at_q(I,J) = IOB%u_flux(i-i0,j-j0) * CS%wind_stress_multiplier
if (associated(IOB%v_flux)) tauy_at_q(I,J) = IOB%v_flux(i-i0,j-j0) * CS%wind_stress_multiplier
Expand Down Expand Up @@ -764,13 +776,10 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, CS)

! sea ice related fields
if (CS%rigid_sea_ice) then
! The commented out code here and in the following lines is the correct
! version, but the incorrect version is being retained temporarily to avoid
! changing answers.
call pass_var(forces%p_surf_full, G%Domain)
call pass_var(forces%p_surf_full, G%Domain, halo=1)
I_GEarth = 1.0 / G%G_Earth
Kv_rho_ice = (CS%kv_sea_ice / CS%density_sea_ice)
do I=isd,ied-1 ; do j=jsd,jed
do I=is-1,ie ; do j=js,je
mass_ice = min(forces%p_surf_full(i,j), forces%p_surf_full(i+1,j)) * I_GEarth
mass_eff = 0.0
if (mass_ice > CS%rigid_sea_ice_mass) then
Expand All @@ -781,7 +790,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, CS)
! a maximum for the second call.
forces%rigidity_ice_u(I,j) = Kv_rho_ice * mass_eff
enddo ; enddo
do i=isd,ied ; do J=jsd,jed-1
do i=is,ie ; do J=js-1,je
mass_ice = min(forces%p_surf_full(i,j), forces%p_surf_full(i,j+1)) * I_GEarth
mass_eff = 0.0
if (mass_ice > CS%rigid_sea_ice_mass) then
Expand Down
40 changes: 19 additions & 21 deletions config_src/coupled_driver/ocean_model_MOM.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,10 +24,12 @@ module ocean_model_mod
use MOM_constants, only : CELSIUS_KELVIN_OFFSET, hlf
use MOM_diag_mediator, only : diag_ctrl, enable_averaging, disable_averaging
use MOM_diag_mediator, only : diag_mediator_close_registration, diag_mediator_end
use MOM_domains, only : pass_vector, AGRID, BGRID_NE, CGRID_NE
use MOM_domains, only : pass_var, pass_vector, AGRID, BGRID_NE, CGRID_NE
use MOM_domains, only : TO_ALL, Omit_Corners
use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe
use MOM_error_handler, only : callTree_enter, callTree_leave
use MOM_file_parser, only : get_param, log_version, close_param_file, param_file_type
use MOM_forcing_type, only : allocate_forcing_type
use MOM_forcing_type, only : forcing, mech_forcing
use MOM_forcing_type, only : forcing_accumulate, copy_common_forcing_fields
use MOM_forcing_type, only : copy_back_forcing_fields, set_net_mass_forcing
Expand Down Expand Up @@ -59,10 +61,8 @@ module ocean_model_mod
use mpp_domains_mod, only : domain2d, mpp_get_layout, mpp_get_global_domain
use mpp_domains_mod, only : mpp_define_domains, mpp_get_compute_domain, mpp_get_data_domain
use atmos_ocean_fluxes_mod, only : aof_set_coupler_flux
use MOM_forcing_type, only : allocate_forcing_type
use fms_mod, only : stdout
use mpp_mod, only : mpp_chksum
use MOM_domains, only : pass_var, pass_vector, TO_ALL, CGRID_NE, BGRID_NE
use MOM_EOS, only : gsw_sp_from_sr, gsw_pt_from_ct
use MOM_wave_interface, only: wave_parameters_CS, MOM_wave_interface_init
use MOM_wave_interface, only: MOM_wave_interface_init_lite, Update_Surface_Waves
Expand Down Expand Up @@ -540,7 +540,6 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
call shelf_calc_flux(OS%sfc_state, OS%forces, OS%fluxes, OS%Time, dt_coupling, OS%Ice_shelf_CSp)
endif
if (OS%icebergs_apply_rigid_boundary) then
!This assumes that the iceshelf and ocean are on the same grid. I hope this is true
call add_berg_flux_to_shelf(OS%grid, OS%forces, OS%fluxes, OS%use_ice_shelf, &
OS%density_iceberg, OS%kv_iceberg, OS%latent_heat_fusion, OS%sfc_state, &
dt_coupling, OS%berg_area_threshold)
Expand All @@ -565,7 +564,6 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
call shelf_calc_flux(OS%sfc_state, OS%forces, OS%flux_tmp, OS%Time, dt_coupling, OS%Ice_shelf_CSp)
endif
if (OS%icebergs_apply_rigid_boundary) then
! This assumes that the iceshelf and ocean are on the same grid. I hope this is true
call add_berg_flux_to_shelf(OS%grid, OS%forces, OS%flux_tmp, OS%use_ice_shelf, OS%density_iceberg, &
OS%kv_iceberg, OS%latent_heat_fusion, OS%sfc_state, dt_coupling, OS%berg_area_threshold)
endif
Expand Down Expand Up @@ -728,43 +726,43 @@ subroutine add_berg_flux_to_shelf(G, forces, fluxes, use_ice_shelf, density_ice,
!the ocean model. This routine is taken from the add_shelf_flux subroutine
!within the ice shelf model.

if (.not.(associated(fluxes%area_berg) .and. associated(fluxes%ustar_berg) .and. &
associated(fluxes%mass_berg) ) ) return
if (.not.(associated(forces%area_berg) .and. associated(forces%mass_berg) ) ) return

if (.not.(associated(forces%frac_shelf_u) .and. associated(forces%frac_shelf_v) .and. &
associated(forces%rigidity_ice_u) .and. associated(forces%rigidity_ice_v)) ) return

if (.not.(associated(fluxes%area_berg) .and. associated(fluxes%ustar_berg) .and. &
associated(fluxes%mass_berg) ) ) return
if (.not.(associated(fluxes%frac_shelf_h) .and. associated(fluxes%ustar_shelf)) ) return

! This section sets or augments the values of fields in forces.
if (.not. use_ice_shelf) then
forces%frac_shelf_u(:,:) = 0.
forces%frac_shelf_v(:,:) = 0.
forces%rigidity_ice_u(:,:) = 0.
forces%rigidity_ice_v(:,:) = 0.
forces%frac_shelf_u(:,:) = 0.0 ; forces%frac_shelf_v(:,:) = 0.0
forces%rigidity_ice_u(:,:) = 0.0 ; forces%rigidity_ice_v(:,:) = 0.0
endif

call pass_var(forces%area_berg, G%domain, TO_ALL+Omit_corners, halo=1, complete=.false.)
call pass_var(forces%mass_berg, G%domain, TO_ALL+Omit_corners, halo=1, complete=.true.)
kv_rho_ice = kv_ice / density_ice
do j=jsd,jed ; do I=isd,ied-1
forces%frac_shelf_u(I,j) = 0.0
do j=js,je ; do I=is-1,ie
if ((G%areaT(i,j) + G%areaT(i+1,j) > 0.0)) & ! .and. (G%dxdy_u(I,j) > 0.0)) &
forces%frac_shelf_u(I,j) = forces%frac_shelf_u(I,j) + &
(((fluxes%area_berg(i,j)*G%areaT(i,j)) + &
(fluxes%area_berg(i+1,j)*G%areaT(i+1,j))) / &
(((forces%area_berg(i,j)*G%areaT(i,j)) + &
(forces%area_berg(i+1,j)*G%areaT(i+1,j))) / &
(G%areaT(i,j) + G%areaT(i+1,j)) )
forces%rigidity_ice_u(I,j) = forces%rigidity_ice_u(I,j) + kv_rho_ice * &
min(fluxes%mass_berg(i,j), fluxes%mass_berg(i+1,j))
min(forces%mass_berg(i,j), forces%mass_berg(i+1,j))
enddo ; enddo
do J=jsd,jed-1 ; do i=isd,ied
forces%frac_shelf_v(i,J) = 0.0
do J=js-1,je ; do i=is,ie
if ((G%areaT(i,j) + G%areaT(i,j+1) > 0.0)) & ! .and. (G%dxdy_v(i,J) > 0.0)) &
forces%frac_shelf_v(i,J) = forces%frac_shelf_v(i,J) + &
(((fluxes%area_berg(i,j)*G%areaT(i,j)) + &
(fluxes%area_berg(i,j+1)*G%areaT(i,j+1))) / &
(((forces%area_berg(i,j)*G%areaT(i,j)) + &
(forces%area_berg(i,j+1)*G%areaT(i,j+1))) / &
(G%areaT(i,j) + G%areaT(i,j+1)) )
forces%rigidity_ice_v(i,J) = forces%rigidity_ice_v(i,J) + kv_rho_ice * &
min(fluxes%mass_berg(i,j), fluxes%mass_berg(i,j+1))
min(forces%mass_berg(i,j), forces%mass_berg(i,j+1))
enddo ; enddo
!### This halo update may be unnecessary. Test it. -RWH
call pass_vector(forces%frac_shelf_u, forces%frac_shelf_v, G%domain, TO_ALL, CGRID_NE)

! The remaining code sets or augments the values of fields in fluxes.
Expand Down

0 comments on commit e7ffe07

Please sign in to comment.