From e7ffe07f07743bed95eb97b26b6b41f72efb4612 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 1 May 2018 13:14:54 -0400 Subject: [PATCH] Store iceberg properties in forces 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. --- .../coupled_driver/MOM_surface_forcing.F90 | 31 +++++++++----- config_src/coupled_driver/ocean_model_MOM.F90 | 40 +++++++++---------- 2 files changed, 39 insertions(+), 32 deletions(-) diff --git a/config_src/coupled_driver/MOM_surface_forcing.F90 b/config_src/coupled_driver/MOM_surface_forcing.F90 index 7de381f8c0..f832e76a8d 100644 --- a/config_src/coupled_driver/MOM_surface_forcing.F90 +++ b/config_src/coupled_driver/MOM_surface_forcing.F90 @@ -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, & @@ -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) @@ -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 @@ -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 @@ -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 @@ -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 diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index db59637ca6..5fd45f7ea4 100644 --- a/config_src/coupled_driver/ocean_model_MOM.F90 +++ b/config_src/coupled_driver/ocean_model_MOM.F90 @@ -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 @@ -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 @@ -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) @@ -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 @@ -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.