Skip to content

User/bgr/homogenize iob #1480

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
breichl wants to merge 2 commits into from
Closed

User/bgr/homogenize iob #1480

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
8c09c14
Adds option to homogenize IOB fields in FMS_cap/ocean_model_MOM
Aug 26, 2021
40368f9
Merge branch 'dev/gfdl' of https://github.com/NOAA-GFDL/MOM6 into HEAD
Aug 26, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
214 changes: 205 additions & 9 deletions config_src/drivers/FMS_cap/ocean_model_MOM.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@ module ocean_model_mod
use MOM_diag_mediator, only : diag_mediator_close_registration, diag_mediator_end
use MOM_domains, only : MOM_domain_type, domain2d, clone_MOM_domain, get_domain_extent
use MOM_domains, only : pass_var, pass_vector, AGRID, BGRID_NE, CGRID_NE, TO_ALL, Omit_Corners
use MOM_domains, only : sum_across_PEs
use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_pe
use MOM_error_handler, only : callTree_enter, callTree_leave
use MOM_EOS, only : gsw_sp_from_sr, gsw_pt_from_ct
Expand Down Expand Up @@ -164,7 +165,8 @@ module ocean_model_mod
!! fields necessary to integrate only the tracer advection
!! and diffusion equation read in from files stored from
!! a previous integration of the prognostic model.

logical :: use_homogenized_IOB !< Creates a homogenized version of the IOB structure
!! that yields homogeneous forcing/fluxes.
logical :: single_step_call !< If true, advance the state of MOM with a single
!! step including both dynamics and thermodynamics.
!! If false, the two phases are advanced with
Expand Down Expand Up @@ -358,6 +360,10 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, wind_stagger, gas
use_melt_pot=.false.
endif

call get_param(param_file, mdl, "HOMOGENIZE_IOB_FORCINGS", OS%use_homogenized_IOB, &
"If true, creates a copy of all IOB fluxes and averages them "//&
"horizontally across the entire domain.", default=.false.)

!allocate(OS%sfc_state)
call allocate_surface_state(OS%sfc_state, OS%grid, use_temperature, do_integrals=.true., &
gas_fields_ocn=gas_fields_ocn, use_meltpot=use_melt_pot)
Expand Down Expand Up @@ -446,6 +452,8 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda
real, optional, intent(in) :: cycle_length !< The duration of a coupled time stepping cycle [s].

! Local variables
type(ice_ocean_boundary_type) :: Ice_ocean_boundary_homogenized ! A copy of the structure containing the
! forcing fields coming from the ice and atmosphere, that is horizontally averaged.
type(time_type) :: Time_seg_start ! Stores the dynamic or thermodynamic ocean model time at the
! start of this call to allow step_MOM to temporarily change the time
! as seen by internal modules.
Expand Down Expand Up @@ -504,9 +512,17 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda
! Translate Ice_ocean_boundary into fluxes and forces.
call get_domain_extent(Ocean_sfc%Domain, index_bnds(1), index_bnds(2), index_bnds(3), index_bnds(4))

if (OS%use_homogenized_IOB) &
call homogenize_IOB(Ice_ocean_boundary, Ice_ocean_boundary_homogenized, OS%grid, index_bnds)
Copy link
Copy Markdown
Contributor

@nikizadehgfdl nikizadehgfdl Sep 3, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@breichl , what is the reason for doing this on an extra copy of IOB? I think if we homogenize the fields on the original Ice_ocean_boundary (rather than the new Ice_ocean_boundary_homogenized ) everything would work fine, including the (generic) tracer fluxes. The problem with Ice_ocean_boundary_homogenized is that it does not have the hidden member "fluxes" allocated.

Copy link
Copy Markdown
Contributor

@nikizadehgfdl nikizadehgfdl Sep 3, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

One issue with modifying IOB itself is that it is intent(in) in update_ocean_model, I had to make it intent(inout). But after that all seem to have worked fine, including the homogenization and the generic tracers.

Copy link
Copy Markdown
Collaborator Author

@breichl breichl Sep 3, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay, thanks Niki. I made the copy to avoid modifying this structure since it is intent(in). If we are comfortable making it intent(inout) then this solution works for me. Do the generic tracers fields also need to be homogenized? I'm not sure how difficult that would be or if it is necessary here.

Copy link
Copy Markdown
Collaborator

@Hallberg-NOAA Hallberg-NOAA Sep 8, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am very uncomfortable making the IOB type intent(inout). The intent is there for a reason, as it helps to document which fields are coming from where, and it does so in a way that the compiler can take advantage of for performance.


if (do_dyn) then
call convert_IOB_to_forces(Ice_ocean_boundary, OS%forces, index_bnds, OS%Time_dyn, OS%grid, OS%US, &
OS%forcing_CSp, dt_forcing=dt_coupling, reset_avg=OS%fluxes%fluxes_used)
if (OS%use_homogenized_IOB) then
call convert_IOB_to_forces(Ice_ocean_boundary_homogenized, OS%forces, index_bnds, OS%Time_dyn, OS%grid, OS%US, &
OS%forcing_CSp, dt_forcing=dt_coupling, reset_avg=OS%fluxes%fluxes_used)
else
call convert_IOB_to_forces(Ice_ocean_boundary, OS%forces, index_bnds, OS%Time_dyn, OS%grid, OS%US, &
OS%forcing_CSp, dt_forcing=dt_coupling, reset_avg=OS%fluxes%fluxes_used)
endif
if (OS%use_ice_shelf) &
call add_shelf_forces(OS%grid, OS%US, OS%Ice_shelf_CSp, OS%forces)
if (OS%icebergs_alter_ocean) &
Expand All @@ -516,9 +532,13 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda

if (do_thermo) then
if (OS%fluxes%fluxes_used) then
call convert_IOB_to_fluxes(Ice_ocean_boundary, OS%fluxes, index_bnds, OS%Time, dt_coupling, &
OS%grid, OS%US, OS%forcing_CSp, OS%sfc_state)

if (OS%use_homogenized_IOB) then
call convert_IOB_to_fluxes(Ice_ocean_boundary_homogenized, OS%fluxes, index_bnds, OS%Time, dt_coupling, &
OS%grid, OS%US, OS%forcing_CSp, OS%sfc_state)
else
call convert_IOB_to_fluxes(Ice_ocean_boundary, OS%fluxes, index_bnds, OS%Time, dt_coupling, &
OS%grid, OS%US, OS%forcing_CSp, OS%sfc_state)
endif
! Add ice shelf fluxes
if (OS%use_ice_shelf) &
call shelf_calc_flux(OS%sfc_state, OS%fluxes, OS%Time, dt_coupling, OS%Ice_shelf_CSp)
Expand All @@ -534,9 +554,15 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda
else
! The previous fluxes have not been used yet, so translate the input fluxes
! into a temporary type and then accumulate them in about 20 lines.
OS%flux_tmp%C_p = OS%fluxes%C_p
call convert_IOB_to_fluxes(Ice_ocean_boundary, OS%flux_tmp, index_bnds, OS%Time, dt_coupling, &
OS%grid, OS%US, OS%forcing_CSp, OS%sfc_state)
if (OS%use_homogenized_IOB) then
OS%flux_tmp%C_p = OS%fluxes%C_p
call convert_IOB_to_fluxes(Ice_ocean_boundary_homogenized, OS%flux_tmp, index_bnds, OS%Time, dt_coupling, &
OS%grid, OS%US, OS%forcing_CSp, OS%sfc_state)
else
OS%flux_tmp%C_p = OS%fluxes%C_p
call convert_IOB_to_fluxes(Ice_ocean_boundary, OS%flux_tmp, index_bnds, OS%Time, dt_coupling, &
OS%grid, OS%US, OS%forcing_CSp, OS%sfc_state)
endif

if (OS%use_ice_shelf) &
call shelf_calc_flux(OS%sfc_state, OS%flux_tmp, OS%Time, dt_coupling, OS%Ice_shelf_CSp)
Expand Down Expand Up @@ -1194,4 +1220,174 @@ subroutine ocean_model_get_UV_surf(OS, Ocean, name, array2D, isc, jsc)

end subroutine ocean_model_get_UV_surf

!> Homogenizes all associated IOB fields.
subroutine homogenize_IOB(IOB,IOB_h, G, index_bounds)
type(ice_ocean_boundary_type), &
target, intent(in) :: IOB !< An ice-ocean boundary type with fluxes to drive
!! the ocean in a coupled model
type(ice_ocean_boundary_type), &
target, intent(out) :: IOB_h !< Similar to IOB, but with homogenized fields.
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
integer, dimension(4), intent(in) :: index_bounds !< The i- and j- size of the arrays in IOB.

integer :: is, ie, js, je, i0, j0, isc_bnd, jsc_bnd, iec_bnd, jec_bnd, i, j
real :: tmp_sum, tmp_mean, N

is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
isc_bnd = index_bounds(1) ; iec_bnd = index_bounds(2)
jsc_bnd = index_bounds(3) ; jec_bnd = index_bounds(4)

i0 = is - isc_bnd ; j0 = js - jsc_bnd

if (associated(IOB%u_flux)) then
allocate(IOB_h%u_flux (isc_bnd:iec_bnd,jsc_bnd:jec_bnd) ) ; IOB_h%u_flux = 0.0
call homogenize_field(IOB%u_flux,IOB_h%u_flux, G, index_bounds)
endif

if (associated(IOB%v_flux)) then
allocate(IOB_h%v_flux (isc_bnd:iec_bnd,jsc_bnd:jec_bnd) ) ; IOB_h%v_flux = 0.0
call homogenize_field(IOB%v_flux,IOB_h%v_flux, G, index_bounds)
endif

if (associated(IOB%t_flux)) then
allocate(IOB_h%t_flux (isc_bnd:iec_bnd,jsc_bnd:jec_bnd) ) ; IOB_h%t_flux = 0.0
call homogenize_field(IOB%t_flux,IOB_h%t_flux, G, index_bounds)
endif

if (associated(IOB%q_flux)) then
allocate(IOB_h%q_flux (isc_bnd:iec_bnd,jsc_bnd:jec_bnd) ) ; IOB_h%q_flux = 0.0
call homogenize_field(IOB%q_flux,IOB_h%q_flux, G, index_bounds)
endif

if (associated(IOB%salt_flux)) then
allocate(IOB_h%salt_flux (isc_bnd:iec_bnd,jsc_bnd:jec_bnd) ) ; IOB_h%salt_flux = 0.0
call homogenize_field(IOB%salt_flux,IOB_h%salt_flux, G, index_bounds)
endif


if (associated(IOB%lw_flux)) then
allocate(IOB_h%lw_flux (isc_bnd:iec_bnd,jsc_bnd:jec_bnd) ) ; IOB_h%lw_flux = 0.0
call homogenize_field(IOB%lw_flux,IOB_h%lw_flux, G, index_bounds)
endif

if (associated(IOB%sw_flux_vis_dir)) then
allocate(IOB_h%sw_flux_vis_dir (isc_bnd:iec_bnd,jsc_bnd:jec_bnd) ) ; IOB_h%sw_flux_vis_dir = 0.0
call homogenize_field(IOB%sw_flux_vis_dir,IOB_h%sw_flux_vis_dir, G, index_bounds)
endif

if (associated(IOB%sw_flux_vis_dif)) then
allocate(IOB_h%sw_flux_vis_dif (isc_bnd:iec_bnd,jsc_bnd:jec_bnd) ) ; IOB_h%sw_flux_vis_dif = 0.0
call homogenize_field(IOB%sw_flux_vis_dif,IOB_h%sw_flux_vis_dif, G, index_bounds)
endif

if (associated(IOB%sw_flux_nir_dir)) then
allocate(IOB_h%sw_flux_nir_dir (isc_bnd:iec_bnd,jsc_bnd:jec_bnd) ) ; IOB_h%sw_flux_nir_dir = 0.0
call homogenize_field(IOB%sw_flux_nir_dir,IOB_h%sw_flux_nir_dir, G, index_bounds)
endif

if (associated(IOB%sw_flux_nir_dif)) then
allocate(IOB_h%sw_flux_nir_dif (isc_bnd:iec_bnd,jsc_bnd:jec_bnd) ) ; IOB_h%sw_flux_nir_dif = 0.0
call homogenize_field(IOB%sw_flux_nir_dif,IOB_h%sw_flux_nir_dif, G, index_bounds)
endif

if (associated(IOB%lprec)) then
allocate(IOB_h%lprec (isc_bnd:iec_bnd,jsc_bnd:jec_bnd) ) ; IOB_h%lprec = 0.0
call homogenize_field(IOB%lprec,IOB_h%lprec, G, index_bounds)
endif

if (associated(IOB%fprec)) then
allocate(IOB_h%fprec (isc_bnd:iec_bnd,jsc_bnd:jec_bnd) ) ; IOB_h%fprec = 0.0
call homogenize_field(IOB%fprec,IOB_h%fprec, G, index_bounds)
endif

if (associated(IOB%runoff)) then
allocate(IOB_h%runoff (isc_bnd:iec_bnd,jsc_bnd:jec_bnd) ) ; IOB_h%runoff = 0.0
call homogenize_field(IOB%runoff,IOB_h%runoff, G, index_bounds)
endif

if (associated(IOB%calving)) then
allocate(IOB_h%calving (isc_bnd:iec_bnd,jsc_bnd:jec_bnd) ) ; IOB_h%calving = 0.0
call homogenize_field(IOB%calving,IOB_h%calving, G, index_bounds)
endif

if (associated(IOB%stress_mag)) then
allocate(IOB_h%stress_mag (isc_bnd:iec_bnd,jsc_bnd:jec_bnd) ) ; IOB_h%stress_mag = 0.0
call homogenize_field(IOB%stress_mag,IOB_h%stress_mag, G, index_bounds)
endif

if (associated(IOB%ustar_berg)) then
allocate(IOB_h%ustar_berg (isc_bnd:iec_bnd,jsc_bnd:jec_bnd) ) ; IOB_h%ustar_berg = 0.0
call homogenize_field(IOB%ustar_berg,IOB_h%ustar_berg, G, index_bounds)
endif

if (associated(IOB%area_berg)) then
allocate(IOB_h%area_berg (isc_bnd:iec_bnd,jsc_bnd:jec_bnd) ) ; IOB_h%area_berg = 0.0
call homogenize_field(IOB%area_berg,IOB_h%area_berg, G, index_bounds)
endif

if (associated(IOB%mass_berg)) then
allocate(IOB_h%mass_berg (isc_bnd:iec_bnd,jsc_bnd:jec_bnd) ) ; IOB_h%mass_berg = 0.0
call homogenize_field(IOB%mass_berg,IOB_h%mass_berg, G, index_bounds)
endif

if (associated(IOB%runoff_hflx)) then
allocate(IOB_h%runoff_hflx (isc_bnd:iec_bnd,jsc_bnd:jec_bnd) ) ; IOB_h%runoff_hflx = 0.0
call homogenize_field(IOB%runoff_hflx,IOB_h%runoff_hflx, G, index_bounds)
endif

if (associated(IOB%calving_hflx)) then
allocate(IOB_h%calving_hflx (isc_bnd:iec_bnd,jsc_bnd:jec_bnd) ) ; IOB_h%calving_hflx = 0.0
call homogenize_field(IOB%calving_hflx,IOB_h%calving_hflx, G, index_bounds)
endif

if (associated(IOB%p)) then
allocate(IOB_h%p (isc_bnd:iec_bnd,jsc_bnd:jec_bnd) ) ; IOB_h%p = 0.0
call homogenize_field(IOB%p,IOB_h%p, G, index_bounds)
endif

if (associated(IOB%mi)) then
allocate(IOB_h%mi (isc_bnd:iec_bnd,jsc_bnd:jec_bnd) ) ; IOB_h%mi = 0.0
call homogenize_field(IOB%mi,IOB_h%mi, G, index_bounds)
endif

if (associated(IOB%ice_rigidity)) then
allocate(IOB_h%ice_rigidity (isc_bnd:iec_bnd,jsc_bnd:jec_bnd) ) ; IOB_h%ice_rigidity = 0.0
call homogenize_field(IOB%ice_rigidity,IOB_h%ice_rigidity, G, index_bounds)
endif

Copy link
Copy Markdown
Contributor

@nikizadehgfdl nikizadehgfdl Sep 1, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@breichl , what about other IOBs? When I run this code change with "generic tracers" turned on the model crashes with

FATAL: extract_coupler_values: The requested boundary condition is not associated.

Image              PC                Routine            Line        Source             
MOM6SIS2           00000000026A4A77  Unknown               Unknown  Unknown
MOM6SIS2           0000000001A421EC  mpp_mod_mp_mpp_er          68  mpp_util_mpi.inc
MOM6SIS2           00000000017767FC  fms_coupler_util_          53  FMS_coupler_util.F90
MOM6SIS2           00000000004E20B6  g_tracer_utils_mp        1429  generic_tracer_utils.F90
MOM6SIS2           00000000011995D0  ocean_model_mod_m         551  ocean_model_MOM.F90
MOM6SIS2           00000000013E3A6F  MAIN__                   1041  coupler_main.F90

Copy link
Copy Markdown
Collaborator

@Hallberg-NOAA Hallberg-NOAA Sep 8, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think that this issue that Niki is encountering suggests that it might be worth revisiting whether it would be simpler to homogenize the fields in the forcing and mech_forcing types, even if we had to add the option to recalculate the dependent fields (e.g., ustar) based on the primary fields (like taux and tauy). We should not be changing the intent on the various types.

end subroutine homogenize_IOB

!> Homogenize fields for homogenize_IOB
subroutine homogenize_field(field_in,field_out,G, index_bounds)
real, dimension(:,:), intent(in) :: field_in !< Input 2d array meant to be homogenized
real, dimension(:,:), intent(out) :: field_out !< Output 2d array that is homogenized
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
integer, dimension(4), intent(in) :: index_bounds !< The i- and j- size of the arrays in IOB.

real :: tmp_sum, N, tmp_mean
integer :: is, ie, js, je, isc_bnd, jsc_bnd
integer :: i0, j0, i, j

is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
isc_bnd = index_bounds(1)
jsc_bnd = index_bounds(3)

i0 = is - isc_bnd ; j0 = js - jsc_bnd

tmp_sum = 0.0 ; N = 0.0
do j=js,je ; do i=is,ie
tmp_sum = tmp_sum + field_in(i-i0,j-j0) * G%mask2dT(i,j)
N = N + G%mask2dT(i,j)
enddo; enddo

call sum_across_PEs(tmp_sum)
call sum_across_PEs(N)

tmp_mean = tmp_sum/N
do j=js,je ; do i=is,ie
field_out(i-i0,j-j0) = tmp_mean
enddo; enddo

end subroutine homogenize_field

end module ocean_model_mod