From c0496e1a0e8448d2ad6ac63b100a13d5315608bc Mon Sep 17 00:00:00 2001 From: "brandon.reichl" Date: Mon, 8 Nov 2021 16:51:30 -0500 Subject: [PATCH 1/4] Adds option to homogenize forces and fluxes fields - Adds functions to do global averages on U and V grids in MOM_spatial_means - Adds functionality to average all forcing and fluxes fields in MOM_forcing_types - Adds flag to average all forcing and fluxes in MOM.F90 - This new feature is primarily for running in single column like configurations with the coupler, which requires perfectly equal forcing across all cells. --- src/core/MOM.F90 | 11 + src/core/MOM_forcing_type.F90 | 416 ++++++++++++++++++++++++++ src/diagnostics/MOM_spatial_means.F90 | 50 +++- 3 files changed, 476 insertions(+), 1 deletion(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index fd798076fa..00a7a3dedf 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -80,6 +80,7 @@ module MOM use MOM_forcing_type, only : allocate_forcing_type, allocate_mech_forcing use MOM_forcing_type, only : deallocate_mech_forcing, deallocate_forcing_type use MOM_forcing_type, only : rotate_forcing, rotate_mech_forcing +use MOM_forcing_type, only : homogenize_forcing, homogenize_mech_forcing use MOM_grid, only : ocean_grid_type, MOM_grid_init, MOM_grid_end use MOM_grid, only : set_first_direction, rescale_grid_bathymetry use MOM_hor_index, only : hor_index_type, hor_index_init @@ -204,6 +205,7 @@ module MOM type(ocean_grid_type) :: G_in !< Input grid metric type(ocean_grid_type), pointer :: G => NULL() !< Model grid metric logical :: rotate_index = .false. !< True if index map is rotated + logical :: homogenize_forcings = .false. !< True if all inputs are homogenized type(verticalGrid_type), pointer :: & GV => NULL() !< structure containing vertical grid info @@ -568,6 +570,12 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS fluxes => fluxes_in endif + ! Homogenize the forces + if (CS%homogenize_forcings) then + call homogenize_mech_forcing(forces, G) + call homogenize_forcing(fluxes, G) + endif + ! First determine the time step that is consistent with this call and an ! integer fraction of time_interval. if (do_dyn) then @@ -2114,6 +2122,9 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & call callTree_waypoint("MOM parameters read (initialize_MOM)") + call get_param(param_file, "MOM", "HOMOGENIZE_FORCINGS", CS%homogenize_forcings, & + "Homogenize the forces and fluxes.", default=.false.) + ! Grid rotation test call get_param(param_file, "MOM", "ROTATE_INDEX", CS%rotate_index, & "Enable rotation of the horizontal indices.", default=.false., & diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 2429ce9d2d..a155ff1c7a 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -17,6 +17,7 @@ module MOM_forcing_type use MOM_grid, only : ocean_grid_type use MOM_opacity, only : sumSWoverBands, optics_type, extract_optics_slice, optics_nbands use MOM_spatial_means, only : global_area_integral, global_area_mean +use MOM_spatial_means, only : global_area_mean_u, global_area_mean_v use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : surface, thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type @@ -35,6 +36,7 @@ module MOM_forcing_type public set_derived_forcing_fields, copy_back_forcing_fields public set_net_mass_forcing, get_net_mass_forcing public rotate_forcing, rotate_mech_forcing +public homogenize_forcing, homogenize_mech_forcing !> Allocate the fields of a (flux) forcing type, based on either a set of input !! flags for each group of fields, or a pre-allocated reference forcing. @@ -3385,6 +3387,7 @@ subroutine rotate_forcing(fluxes_in, fluxes, turns) if (do_iceberg) then call rotate_array(fluxes_in%ustar_berg, turns, fluxes%ustar_berg) call rotate_array(fluxes_in%area_berg, turns, fluxes%area_berg) + !BGR: pretty sure the following line isn't supposed to be here. call rotate_array(fluxes_in%iceshelf_melt, turns, fluxes%iceshelf_melt) endif @@ -3490,6 +3493,419 @@ subroutine rotate_mech_forcing(forces_in, turns, forces) forces%initialized = forces_in%initialized end subroutine rotate_mech_forcing +!< Homogenize the forcing fields from the input domain +subroutine homogenize_mech_forcing(forces, G) + type(mech_forcing), intent(in) :: forces !< Forcing on the input domain + type(ocean_grid_type), intent(in) :: G !< Grid metric of target forcing + + real :: tx_mean, ty_mean, avg + logical :: do_stress, do_ustar, do_shelf, do_press, do_iceberg + integer :: i, j, is, ie, js, je, isB, ieB, jsB, jeB + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + isB = G%iscB ; ieB = G%iecB ; jsB = G%jscB ; jeB = G%jecB + + call get_mech_forcing_groups(forces, do_stress, do_ustar, do_shelf, & + do_press, do_iceberg) + + if (do_stress) then + + tx_mean = global_area_mean_u(forces%taux, G) + do j=js,je ; do i=isB,ieB + if (G%mask2dCu(I,j) > 0.) forces%taux(I,j) = tx_mean + enddo ; enddo + + ty_mean = global_area_mean_v(forces%tauy, G) + do j=jsB,jeB ; do i=is,ie + if (G%mask2dCv(i,J) > 0.) forces%tauy(i,J) = ty_mean + enddo ; enddo + + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) forces%ustar(i,j) = sqrt(tx_mean**2 + ty_mean**2) + enddo ; enddo + + else + + if (do_ustar) then + avg = global_area_mean(forces%ustar, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) forces%ustar(i,j) = avg + enddo ; enddo + endif + endif + + + + if (do_shelf) then + + avg = global_area_mean_u(forces%rigidity_ice_u, G) + do j=js,je ; do i=isB,ieB + if (G%mask2dCu(I,j) > 0.) forces%rigidity_ice_u(I,j) = avg + enddo ; enddo + + avg = global_area_mean_v(forces%rigidity_ice_v, G) + do j=jsB,jeB ; do i=is,ie + if (G%mask2dCv(i,j) > 0.) forces%rigidity_ice_v(i,J) = avg + enddo ; enddo + + avg = global_area_mean_u(forces%frac_shelf_u, G) + do j=js,je ; do i=isB,ieB + if (G%mask2dCu(i,j) > 0.) forces%frac_shelf_u(I,j) = avg + enddo ; enddo + + avg = global_area_mean_v(forces%frac_shelf_v, G) + do j=jsB,jeB ; do i=is,ie + if (G%mask2dCv(i,j) > 0.) forces%frac_shelf_v(i,J) = avg + enddo ; enddo + + endif + + if (do_press) then + + ! NOTE: p_surf_SSH either points to p_surf or p_surf_full + + avg = global_area_mean(forces%p_surf, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) forces%p_surf(i,j) = avg + enddo ; enddo + + avg = global_area_mean(forces%p_surf_full, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) forces%p_surf_full(i,j) = avg + enddo; enddo + + avg = global_area_mean(forces%net_mass_src, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) forces%net_mass_src(i,j) = avg + enddo ; enddo + + endif + + if (do_iceberg) then + + avg = global_area_mean(forces%area_berg, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) forces%area_berg(i,j) = avg + enddo ; enddo + + avg = global_area_mean(forces%mass_berg, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) forces%mass_berg(i,j) = avg + enddo; enddo + + endif + +end subroutine homogenize_mech_forcing + +!< Homogenize the fluxes +subroutine homogenize_forcing(fluxes, G) + type(forcing), intent(inout) :: fluxes !< Input forcing struct + type(ocean_grid_type), intent(in) :: G !< Grid metric of target forcing + + real :: avg + logical :: do_ustar, do_water, do_heat, do_salt, do_press, do_shelf, & + do_iceberg, do_heat_added, do_buoy + integer :: i, j, is, ie, js, je, isB, ieB, jsB, jeB + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + isB = G%iscB ; ieB = G%iecB ; jsB = G%jscB ; jeB = G%jecB + + call get_forcing_groups(fluxes, do_water, do_heat, do_ustar, do_press, & + do_shelf, do_iceberg, do_salt, do_heat_added, do_buoy) + + if (do_ustar) then + + ! This is wrong!!!! + avg = global_area_mean(fluxes%ustar, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%ustar(i,j) = avg + enddo ; enddo + + avg = global_area_mean(fluxes%ustar_gustless, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%ustar_gustless(i,j) = avg + enddo ; enddo + + endif + + if (do_water) then + + avg = global_area_mean(fluxes%evap, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%evap(i,j) = avg + enddo ; enddo + + avg = global_area_mean(fluxes%lprec, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%lprec(i,j) = avg + enddo ; enddo + + avg = global_area_mean(fluxes%fprec, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%fprec (i,j) = avg + enddo ; enddo + + avg = global_area_mean(fluxes%vprec, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%vprec(i,j) = avg + enddo ; enddo + + avg = global_area_mean(fluxes%lrunoff, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%lrunoff(i,j) = avg + enddo ; enddo + + avg = global_area_mean(fluxes%frunoff, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%frunoff(i,j) = avg + enddo ; enddo + + avg = global_area_mean(fluxes%seaice_melt, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%seaice_melt(i,j) = avg + enddo ; enddo + + avg = global_area_mean(fluxes%netMassOut, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%netMassOut(i,j) = avg + enddo ; enddo + + avg = global_area_mean(fluxes%netMassIn, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%netMassIn(i,j) = avg + enddo ; enddo + + avg = global_area_mean(fluxes%netSalt, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%netSalt(i,j) = avg + enddo ; enddo + + endif + + if (do_heat) then + + avg= global_area_mean(fluxes%seaice_melt_heat, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%seaice_melt_heat(i,j) = avg + enddo ; enddo + + avg = global_area_mean(fluxes%sw, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%sw(i,j) = avg + enddo ; enddo + + avg = global_area_mean(fluxes%lw, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%lw(i,j) = avg + enddo ; enddo + + avg = global_area_mean(fluxes%latent, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%latent(i,j) = avg + enddo ; enddo + + avg = global_area_mean(fluxes%sens, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%sens(i,j) = avg + enddo ; enddo + + avg = global_area_mean(fluxes%latent_evap_diag, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%latent_evap_diag(i,j) = avg + enddo ; enddo + + avg = global_area_mean(fluxes%latent_fprec_diag, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%latent_fprec_diag(i,j) = avg + enddo ; enddo + + avg = global_area_mean(fluxes%latent_frunoff_diag, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%latent_frunoff_diag(i,j) = avg + enddo ; enddo + + endif + + if (do_salt) then + avg = global_area_mean(fluxes%salt_flux, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%salt_flux(i,j) = avg + enddo ; enddo + endif + + if (do_heat .and. do_water) then + + avg = global_area_mean(fluxes%heat_content_cond, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%heat_content_cond(i,j) = avg + enddo ; enddo + + avg =global_area_mean(fluxes%heat_content_icemelt, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%heat_content_icemelt(i,j) = avg + enddo ; enddo + + avg =global_area_mean(fluxes%heat_content_lprec, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%heat_content_lprec(i,j) = avg + enddo ; enddo + + avg =global_area_mean(fluxes%heat_content_fprec, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%heat_content_fprec(i,j) = avg + enddo ; enddo + + avg =global_area_mean(fluxes%heat_content_vprec, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%heat_content_vprec(i,j) = avg + enddo ; enddo + + avg =global_area_mean(fluxes%heat_content_lrunoff, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%heat_content_lrunoff(i,j) = avg + enddo ; enddo + + avg =global_area_mean(fluxes%heat_content_frunoff, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%heat_content_frunoff(i,j) = avg + enddo ; enddo + + avg =global_area_mean(fluxes%heat_content_massout, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%heat_content_massout(i,j) = avg + enddo ; enddo + + avg =global_area_mean(fluxes%heat_content_massin, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%heat_content_massin(i,j) = avg + enddo ; enddo + + endif + + if (do_press) then + avg = global_area_mean(fluxes%p_surf, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%p_surf(i,j) = avg + enddo ; enddo + endif + + if (do_shelf) then + + avg = global_area_mean(fluxes%frac_shelf_h, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%frac_shelf_h(i,j) = avg + enddo ; enddo + + avg = global_area_mean(fluxes%ustar_shelf, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%ustar_shelf(i,j) = avg + enddo ; enddo + + avg = global_area_mean(fluxes%iceshelf_melt, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%iceshelf_melt(i,j) = avg + enddo ; enddo + + endif + + if (do_iceberg) then + + avg = global_area_mean(fluxes%ustar_berg, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%ustar_berg (i,j) = avg + enddo ; enddo + + avg = global_area_mean(fluxes%area_berg, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%area_berg (i,j) = avg + enddo ; enddo + + endif + + if (do_heat_added) then + avg = global_area_mean(fluxes%heat_added, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%heat_added(i,j) = avg + enddo ; enddo + endif + + ! The following fields are handled by drivers rather than control flags. + if (associated(fluxes%sw_vis_dir)) then + avg = global_area_mean(fluxes%sw_vis_dir, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%sw_vis_dir(i,j) = avg + enddo ; enddo + endif + + if (associated(fluxes%sw_vis_dif)) then + avg = global_area_mean(fluxes%sw_vis_dif, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%sw_vis_dif(i,j) = avg + enddo ; enddo + endif + + if (associated(fluxes%sw_nir_dir)) then + avg= global_area_mean(fluxes%sw_nir_dir, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%sw_nir_dir(i,j) = avg + enddo ; enddo + endif + + if (associated(fluxes%sw_nir_dif)) then + avg = global_area_mean(fluxes%sw_nir_dif, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%sw_nir_dif(i,j) = avg + enddo ; enddo + endif + + if (associated(fluxes%salt_flux_in)) then + avg = global_area_mean(fluxes%salt_flux_in, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%salt_flux_in(i,j) = avg + enddo ; enddo + endif + + if (associated(fluxes%salt_flux_added)) then + avg = global_area_mean(fluxes%salt_flux_added, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%salt_flux_added(i,j) = avg + enddo ; enddo + endif + + if (associated(fluxes%p_surf_full)) then + avg = global_area_mean(fluxes%p_surf_full, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%p_surf_full(i,j) = avg + enddo ; enddo + endif + + if (associated(fluxes%buoy)) then + avg = global_area_mean(fluxes%buoy, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%buoy(i,j) = avg + enddo ; enddo + endif + + if (associated(fluxes%TKE_tidal)) then + avg = global_area_mean(fluxes%TKE_tidal, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%TKE_tidal(i,j) = avg + enddo ; enddo + endif + + if (associated(fluxes%ustar_tidal)) then + avg = global_area_mean(fluxes%ustar_tidal, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) fluxes%ustar_tidal(i,j) = avg + enddo ; enddo + endif + + ! TODO: tracer flux homogenization + ! Having a warning causes a lot of errors (each time step). + !if (coupler_type_initialized(fluxes%tr_fluxes)) & + ! call MOM_error(WARNING, "Homogenization of tracer BC fluxes not yet implemented.") + +end subroutine homogenize_forcing + + !> \namespace mom_forcing_type !! !! \section section_fluxes Boundary fluxes diff --git a/src/diagnostics/MOM_spatial_means.F90 b/src/diagnostics/MOM_spatial_means.F90 index ffbdc5f810..f68cfe1565 100644 --- a/src/diagnostics/MOM_spatial_means.F90 +++ b/src/diagnostics/MOM_spatial_means.F90 @@ -17,7 +17,7 @@ module MOM_spatial_means #include public :: global_i_mean, global_j_mean -public :: global_area_mean, global_layer_mean +public :: global_area_mean, global_area_mean_u, global_area_mean_v, global_layer_mean public :: global_area_integral public :: global_volume_mean, global_mass_integral public :: adjust_area_mean_to_zero @@ -47,6 +47,54 @@ function global_area_mean(var, G, scale) end function global_area_mean +!> Return the global area mean of a variable. This uses reproducing sums. +function global_area_mean_v(var, G, scale) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + real, dimension(SZI_(G), SZJB_(G)), intent(in) :: var !< The variable to average + real, optional, intent(in) :: scale !< A rescaling factor for the variable + + real, dimension(SZI_(G), SZJB_(G)) :: tmpForSumming + real :: global_area_mean_v + + real :: scalefac ! An overall scaling factor for the areas and variable. + integer :: i, j, is, ie, js, je, isB, ieB, jsB, jeB + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + isB = G%iscB ; ieB = G%iecB ; jsB = G%jscB ; jeB = G%jecB + + scalefac = G%US%L_to_m**2 ; if (present(scale)) scalefac = G%US%L_to_m**2*scale + + tmpForSumming(:,:) = 0. + do J=jsB,jeB ; do i=is,ie + tmpForSumming(i,J) = var(i,J) * (scalefac * G%areaCv(i,J) * G%mask2dCv(i,J)) + enddo ; enddo + global_area_mean_v = reproducing_sum(tmpForSumming) * G%IareaT_global !Need to add IareaCv? + +end function global_area_mean_v + +!> Return the global area mean of a variable on U grid. This uses reproducing sums. +function global_area_mean_u(var, G, scale) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + real, dimension(SZIB_(G), SZJ_(G)), intent(in) :: var !< The variable to average + real, optional, intent(in) :: scale !< A rescaling factor for the variable + + real, dimension(SZI_(G), SZJ_(G)) :: tmpForSumming + real :: global_area_mean_u + + real :: scalefac ! An overall scaling factor for the areas and variable. + integer :: i, j, is, ie, js, je, isB, ieB, jsB, jeB + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + isB = G%iscB ; ieB = G%iecB ; jsB = G%jscB ; jeB = G%jecB + + scalefac = G%US%L_to_m**2 ; if (present(scale)) scalefac = G%US%L_to_m**2*scale + + tmpForSumming(:,:) = 0. + do j=js,je ; do I=isB,ieB + tmpForSumming(I,j) = var(I,j) * (scalefac * G%areaCu(I,j) * G%mask2dCu(I,j)) + enddo ; enddo + global_area_mean_u = reproducing_sum(tmpForSumming) * G%IareaT_global !Need to add IareaCu? + +end function global_area_mean_u + !> Return the global area integral of a variable, by default using the masked area from the !! grid, but an alternate could be used instead. This uses reproducing sums. function global_area_integral(var, G, scale, area) From 8f5a31727ae095d6624437e53387cec2fbeb7932 Mon Sep 17 00:00:00 2001 From: "brandon.reichl" Date: Mon, 8 Nov 2021 17:49:16 -0500 Subject: [PATCH 2/4] Fixing ustar calculation in homogenize_mech_forcing - Adds in irho0 and sqrt that were missing in homogenize mech forcing. --- src/core/MOM.F90 | 2 +- src/core/MOM_forcing_type.F90 | 12 +++++++++--- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 00a7a3dedf..603473fd87 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -572,7 +572,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS ! Homogenize the forces if (CS%homogenize_forcings) then - call homogenize_mech_forcing(forces, G) + call homogenize_mech_forcing(forces, G, US, GV%Rho0) call homogenize_forcing(fluxes, G) endif diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index a155ff1c7a..38d5b13e08 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -3494,16 +3494,22 @@ subroutine rotate_mech_forcing(forces_in, turns, forces) end subroutine rotate_mech_forcing !< Homogenize the forcing fields from the input domain -subroutine homogenize_mech_forcing(forces, G) - type(mech_forcing), intent(in) :: forces !< Forcing on the input domain +subroutine homogenize_mech_forcing(forces, G, US, Rho0) + type(mech_forcing), intent(in) :: forces !< Forcing on the input domain type(ocean_grid_type), intent(in) :: G !< Grid metric of target forcing + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, intent(in) :: Rho0 !< A reference density of seawater [R ~> kg m-3], + !! as used to calculate ustar. real :: tx_mean, ty_mean, avg + real :: iRho0 logical :: do_stress, do_ustar, do_shelf, do_press, do_iceberg integer :: i, j, is, ie, js, je, isB, ieB, jsB, jeB is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec isB = G%iscB ; ieB = G%iecB ; jsB = G%jscB ; jeB = G%jecB + iRho0 = US%L_to_Z / Rho0 + call get_mech_forcing_groups(forces, do_stress, do_ustar, do_shelf, & do_press, do_iceberg) @@ -3520,7 +3526,7 @@ subroutine homogenize_mech_forcing(forces, G) enddo ; enddo do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) forces%ustar(i,j) = sqrt(tx_mean**2 + ty_mean**2) + if (G%mask2dT(i,j) > 0.) forces%ustar(i,j) = sqrt(US%L_to_Z * sqrt(tx_mean**2 + ty_mean**2)*iRho0) enddo ; enddo else From 533dccc6d25ff9283b43bb3072d440c88eee126c Mon Sep 17 00:00:00 2001 From: "brandon.reichl" Date: Tue, 9 Nov 2021 15:48:30 -0500 Subject: [PATCH 3/4] Updates to homogenize_forcings options. - Correct issues in global_area_mean_u and global_area_mean_v to work with symmetric and rotated grids. - Add options to compute mean ustar or compute ustar from mean tau. - Add subroutines to replace averaging blocks in MOM_forcing_type. --- src/core/MOM.F90 | 18 +- src/core/MOM_forcing_type.F90 | 453 ++++++++------------------ src/diagnostics/MOM_spatial_means.F90 | 36 +- 3 files changed, 164 insertions(+), 343 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 603473fd87..8307db076f 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -80,6 +80,7 @@ module MOM use MOM_forcing_type, only : allocate_forcing_type, allocate_mech_forcing use MOM_forcing_type, only : deallocate_mech_forcing, deallocate_forcing_type use MOM_forcing_type, only : rotate_forcing, rotate_mech_forcing +use MOM_forcing_type, only : copy_common_forcing_fields, set_derived_forcing_fields use MOM_forcing_type, only : homogenize_forcing, homogenize_mech_forcing use MOM_grid, only : ocean_grid_type, MOM_grid_init, MOM_grid_end use MOM_grid, only : set_first_direction, rescale_grid_bathymetry @@ -206,6 +207,7 @@ module MOM type(ocean_grid_type), pointer :: G => NULL() !< Model grid metric logical :: rotate_index = .false. !< True if index map is rotated logical :: homogenize_forcings = .false. !< True if all inputs are homogenized + logical :: update_ustar = .false. !< True to update ustar from homogenized tau type(verticalGrid_type), pointer :: & GV => NULL() !< structure containing vertical grid info @@ -572,8 +574,16 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS ! Homogenize the forces if (CS%homogenize_forcings) then - call homogenize_mech_forcing(forces, G, US, GV%Rho0) + ! Homogenize all forcing and fluxes fields. + call homogenize_mech_forcing(forces, G, US, GV%Rho0, CS%update_ustar) + ! Note the following computes the mean ustar as the mean of ustar rather than + ! ustar of the mean of tau. call homogenize_forcing(fluxes, G) + if (CS%update_ustar) then + ! These calls corrects the ustar values + call copy_common_forcing_fields(forces, fluxes, G) + call set_derived_forcing_fields(forces, fluxes, G, US, GV%Rho0) + endif endif ! First determine the time step that is consistent with this call and an @@ -2123,7 +2133,11 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & call callTree_waypoint("MOM parameters read (initialize_MOM)") call get_param(param_file, "MOM", "HOMOGENIZE_FORCINGS", CS%homogenize_forcings, & - "Homogenize the forces and fluxes.", default=.false.) + "If True, homogenize the forces and fluxes.", default=.false.) + call get_param(param_file, "MOM", "UPDATE_USTAR",CS%update_ustar, & + "If True, update ustar from homogenized tau when using the "//& + "HOMOGENIZE_FORCINGS option. Note that this will not work "//& + "with a non-zero gustiness factor.", default=.false.) ! Grid rotation test call get_param(param_file, "MOM", "ROTATE_INDEX", CS%rotate_index, & diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 38d5b13e08..44d5d8c135 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -184,6 +184,7 @@ module MOM_forcing_type !! average of the gustless wind stress. real :: C_p !< heat capacity of seawater [Q degC-1 ~> J kg-1 degC-1]. !! C_p is is the same value as in thermovar_ptrs_type. + real :: gust_const !< Constant unresolved background gustiness for ustar [R L Z T-1 ~> Pa] ! CFC-related arrays needed in the MOM_CFC_cap module real, pointer, dimension(:,:) :: & @@ -3494,110 +3495,69 @@ subroutine rotate_mech_forcing(forces_in, turns, forces) end subroutine rotate_mech_forcing !< Homogenize the forcing fields from the input domain -subroutine homogenize_mech_forcing(forces, G, US, Rho0) - type(mech_forcing), intent(in) :: forces !< Forcing on the input domain - type(ocean_grid_type), intent(in) :: G !< Grid metric of target forcing - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - real, intent(in) :: Rho0 !< A reference density of seawater [R ~> kg m-3], +subroutine homogenize_mech_forcing(forces, G, US, Rho0, UpdateUstar) + type(mech_forcing), intent(inout) :: forces !< Forcing on the input domain + type(ocean_grid_type), intent(in) :: G !< Grid metric of target forcing + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, intent(in) :: Rho0 !< A reference density of seawater [R ~> kg m-3], !! as used to calculate ustar. + logical, optional, intent(in) :: UpdateUstar !< A logical to determine if Ustar should be directly averaged + !! or updated from mean tau. real :: tx_mean, ty_mean, avg real :: iRho0 - logical :: do_stress, do_ustar, do_shelf, do_press, do_iceberg + logical :: do_stress, do_ustar, do_shelf, do_press, do_iceberg, tau2ustar integer :: i, j, is, ie, js, je, isB, ieB, jsB, jeB is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec isB = G%iscB ; ieB = G%iecB ; jsB = G%jscB ; jeB = G%jecB iRho0 = US%L_to_Z / Rho0 + tau2ustar = .false. + if (present(UpdateUstar)) tau2ustar = UpdateUstar + call get_mech_forcing_groups(forces, do_stress, do_ustar, do_shelf, & do_press, do_iceberg) if (do_stress) then - tx_mean = global_area_mean_u(forces%taux, G) do j=js,je ; do i=isB,ieB if (G%mask2dCu(I,j) > 0.) forces%taux(I,j) = tx_mean enddo ; enddo - ty_mean = global_area_mean_v(forces%tauy, G) do j=jsB,jeB ; do i=is,ie if (G%mask2dCv(i,J) > 0.) forces%tauy(i,J) = ty_mean enddo ; enddo - - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) forces%ustar(i,j) = sqrt(US%L_to_Z * sqrt(tx_mean**2 + ty_mean**2)*iRho0) - enddo ; enddo - - else - - if (do_ustar) then - avg = global_area_mean(forces%ustar, G) + if (tau2ustar) then do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) forces%ustar(i,j) = avg + if (G%mask2dT(i,j) > 0.) forces%ustar(i,j) = sqrt(sqrt(tx_mean**2 + ty_mean**2)*iRho0) enddo ; enddo + else + call homogenize_field_t(forces%ustar, G) + endif + else + if (do_ustar) then + call homogenize_field_t(forces%ustar, G) endif endif - - if (do_shelf) then - - avg = global_area_mean_u(forces%rigidity_ice_u, G) - do j=js,je ; do i=isB,ieB - if (G%mask2dCu(I,j) > 0.) forces%rigidity_ice_u(I,j) = avg - enddo ; enddo - - avg = global_area_mean_v(forces%rigidity_ice_v, G) - do j=jsB,jeB ; do i=is,ie - if (G%mask2dCv(i,j) > 0.) forces%rigidity_ice_v(i,J) = avg - enddo ; enddo - - avg = global_area_mean_u(forces%frac_shelf_u, G) - do j=js,je ; do i=isB,ieB - if (G%mask2dCu(i,j) > 0.) forces%frac_shelf_u(I,j) = avg - enddo ; enddo - - avg = global_area_mean_v(forces%frac_shelf_v, G) - do j=jsB,jeB ; do i=is,ie - if (G%mask2dCv(i,j) > 0.) forces%frac_shelf_v(i,J) = avg - enddo ; enddo - + call homogenize_field_u(forces%rigidity_ice_u, G) + call homogenize_field_v(forces%rigidity_ice_v, G) + call homogenize_field_u(forces%frac_shelf_u, G) + call homogenize_field_v(forces%frac_shelf_v, G) endif if (do_press) then - ! NOTE: p_surf_SSH either points to p_surf or p_surf_full - - avg = global_area_mean(forces%p_surf, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) forces%p_surf(i,j) = avg - enddo ; enddo - - avg = global_area_mean(forces%p_surf_full, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) forces%p_surf_full(i,j) = avg - enddo; enddo - - avg = global_area_mean(forces%net_mass_src, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) forces%net_mass_src(i,j) = avg - enddo ; enddo - + call homogenize_field_t(forces%p_surf, G) + call homogenize_field_t(forces%p_surf_full, G) + call homogenize_field_t(forces%net_mass_src, G) endif if (do_iceberg) then - - avg = global_area_mean(forces%area_berg, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) forces%area_berg(i,j) = avg - enddo ; enddo - - avg = global_area_mean(forces%mass_berg, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) forces%mass_berg(i,j) = avg - enddo; enddo - + call homogenize_field_t(forces%area_berg, G) + call homogenize_field_t(forces%mass_berg, G) endif end subroutine homogenize_mech_forcing @@ -3618,291 +3578,96 @@ subroutine homogenize_forcing(fluxes, G) do_shelf, do_iceberg, do_salt, do_heat_added, do_buoy) if (do_ustar) then - - ! This is wrong!!!! - avg = global_area_mean(fluxes%ustar, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%ustar(i,j) = avg - enddo ; enddo - - avg = global_area_mean(fluxes%ustar_gustless, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%ustar_gustless(i,j) = avg - enddo ; enddo - + call homogenize_field_t(fluxes%ustar, G) + call homogenize_field_t(fluxes%ustar_gustless, G) endif if (do_water) then - - avg = global_area_mean(fluxes%evap, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%evap(i,j) = avg - enddo ; enddo - - avg = global_area_mean(fluxes%lprec, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%lprec(i,j) = avg - enddo ; enddo - - avg = global_area_mean(fluxes%fprec, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%fprec (i,j) = avg - enddo ; enddo - - avg = global_area_mean(fluxes%vprec, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%vprec(i,j) = avg - enddo ; enddo - - avg = global_area_mean(fluxes%lrunoff, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%lrunoff(i,j) = avg - enddo ; enddo - - avg = global_area_mean(fluxes%frunoff, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%frunoff(i,j) = avg - enddo ; enddo - - avg = global_area_mean(fluxes%seaice_melt, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%seaice_melt(i,j) = avg - enddo ; enddo - - avg = global_area_mean(fluxes%netMassOut, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%netMassOut(i,j) = avg - enddo ; enddo - - avg = global_area_mean(fluxes%netMassIn, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%netMassIn(i,j) = avg - enddo ; enddo - - avg = global_area_mean(fluxes%netSalt, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%netSalt(i,j) = avg - enddo ; enddo - + call homogenize_field_t(fluxes%evap, G) + call homogenize_field_t(fluxes%lprec, G) + call homogenize_field_t(fluxes%lprec, G) + call homogenize_field_t(fluxes%fprec, G) + call homogenize_field_t(fluxes%vprec, G) + call homogenize_field_t(fluxes%lrunoff, G) + call homogenize_field_t(fluxes%frunoff, G) + call homogenize_field_t(fluxes%seaice_melt, G) + call homogenize_field_t(fluxes%netMassOut, G) + call homogenize_field_t(fluxes%netMassIn, G) + call homogenize_field_t(fluxes%netSalt, G) endif if (do_heat) then - - avg= global_area_mean(fluxes%seaice_melt_heat, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%seaice_melt_heat(i,j) = avg - enddo ; enddo - - avg = global_area_mean(fluxes%sw, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%sw(i,j) = avg - enddo ; enddo - - avg = global_area_mean(fluxes%lw, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%lw(i,j) = avg - enddo ; enddo - - avg = global_area_mean(fluxes%latent, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%latent(i,j) = avg - enddo ; enddo - - avg = global_area_mean(fluxes%sens, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%sens(i,j) = avg - enddo ; enddo - - avg = global_area_mean(fluxes%latent_evap_diag, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%latent_evap_diag(i,j) = avg - enddo ; enddo - - avg = global_area_mean(fluxes%latent_fprec_diag, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%latent_fprec_diag(i,j) = avg - enddo ; enddo - - avg = global_area_mean(fluxes%latent_frunoff_diag, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%latent_frunoff_diag(i,j) = avg - enddo ; enddo - + call homogenize_field_t(fluxes%seaice_melt_heat, G) + call homogenize_field_t(fluxes%sw, G) + call homogenize_field_t(fluxes%lw, G) + call homogenize_field_t(fluxes%latent, G) + call homogenize_field_t(fluxes%sens, G) + call homogenize_field_t(fluxes%latent_evap_diag, G) + call homogenize_field_t(fluxes%latent_fprec_diag, G) + call homogenize_field_t(fluxes%latent_frunoff_diag, G) endif - if (do_salt) then - avg = global_area_mean(fluxes%salt_flux, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%salt_flux(i,j) = avg - enddo ; enddo - endif + if (do_salt) call homogenize_field_t(fluxes%salt_flux, G) if (do_heat .and. do_water) then - - avg = global_area_mean(fluxes%heat_content_cond, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%heat_content_cond(i,j) = avg - enddo ; enddo - - avg =global_area_mean(fluxes%heat_content_icemelt, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%heat_content_icemelt(i,j) = avg - enddo ; enddo - - avg =global_area_mean(fluxes%heat_content_lprec, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%heat_content_lprec(i,j) = avg - enddo ; enddo - - avg =global_area_mean(fluxes%heat_content_fprec, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%heat_content_fprec(i,j) = avg - enddo ; enddo - - avg =global_area_mean(fluxes%heat_content_vprec, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%heat_content_vprec(i,j) = avg - enddo ; enddo - - avg =global_area_mean(fluxes%heat_content_lrunoff, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%heat_content_lrunoff(i,j) = avg - enddo ; enddo - - avg =global_area_mean(fluxes%heat_content_frunoff, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%heat_content_frunoff(i,j) = avg - enddo ; enddo - - avg =global_area_mean(fluxes%heat_content_massout, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%heat_content_massout(i,j) = avg - enddo ; enddo - - avg =global_area_mean(fluxes%heat_content_massin, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%heat_content_massin(i,j) = avg - enddo ; enddo - + call homogenize_field_t(fluxes%heat_content_cond, G) + call homogenize_field_t(fluxes%heat_content_icemelt, G) + call homogenize_field_t(fluxes%heat_content_lprec, G) + call homogenize_field_t(fluxes%heat_content_fprec, G) + call homogenize_field_t(fluxes%heat_content_vprec, G) + call homogenize_field_t(fluxes%heat_content_lrunoff, G) + call homogenize_field_t(fluxes%heat_content_frunoff, G) + call homogenize_field_t(fluxes%heat_content_massout, G) + call homogenize_field_t(fluxes%heat_content_massin, G) endif - if (do_press) then - avg = global_area_mean(fluxes%p_surf, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%p_surf(i,j) = avg - enddo ; enddo - endif + if (do_press) call homogenize_field_t(fluxes%p_surf, G) if (do_shelf) then - - avg = global_area_mean(fluxes%frac_shelf_h, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%frac_shelf_h(i,j) = avg - enddo ; enddo - - avg = global_area_mean(fluxes%ustar_shelf, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%ustar_shelf(i,j) = avg - enddo ; enddo - - avg = global_area_mean(fluxes%iceshelf_melt, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%iceshelf_melt(i,j) = avg - enddo ; enddo - + call homogenize_field_t(fluxes%frac_shelf_h, G) + call homogenize_field_t(fluxes%ustar_shelf, G) + call homogenize_field_t(fluxes%iceshelf_melt, G) endif if (do_iceberg) then - - avg = global_area_mean(fluxes%ustar_berg, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%ustar_berg (i,j) = avg - enddo ; enddo - - avg = global_area_mean(fluxes%area_berg, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%area_berg (i,j) = avg - enddo ; enddo - + call homogenize_field_t(fluxes%ustar_berg, G) + call homogenize_field_t(fluxes%area_berg, G) endif if (do_heat_added) then - avg = global_area_mean(fluxes%heat_added, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%heat_added(i,j) = avg - enddo ; enddo + call homogenize_field_t(fluxes%heat_added, G) endif ! The following fields are handled by drivers rather than control flags. - if (associated(fluxes%sw_vis_dir)) then - avg = global_area_mean(fluxes%sw_vis_dir, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%sw_vis_dir(i,j) = avg - enddo ; enddo - endif + if (associated(fluxes%sw_vis_dir)) & + call homogenize_field_t(fluxes%sw_vis_dir, G) - if (associated(fluxes%sw_vis_dif)) then - avg = global_area_mean(fluxes%sw_vis_dif, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%sw_vis_dif(i,j) = avg - enddo ; enddo - endif + if (associated(fluxes%sw_vis_dif)) & + call homogenize_field_t(fluxes%sw_vis_dif, G) - if (associated(fluxes%sw_nir_dir)) then - avg= global_area_mean(fluxes%sw_nir_dir, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%sw_nir_dir(i,j) = avg - enddo ; enddo - endif + if (associated(fluxes%sw_nir_dir)) & + call homogenize_field_t(fluxes%sw_nir_dir, G) - if (associated(fluxes%sw_nir_dif)) then - avg = global_area_mean(fluxes%sw_nir_dif, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%sw_nir_dif(i,j) = avg - enddo ; enddo - endif + if (associated(fluxes%sw_nir_dif)) & + call homogenize_field_t(fluxes%sw_nir_dif, G) - if (associated(fluxes%salt_flux_in)) then - avg = global_area_mean(fluxes%salt_flux_in, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%salt_flux_in(i,j) = avg - enddo ; enddo - endif + if (associated(fluxes%salt_flux_in)) & + call homogenize_field_t(fluxes%salt_flux_in, G) - if (associated(fluxes%salt_flux_added)) then - avg = global_area_mean(fluxes%salt_flux_added, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%salt_flux_added(i,j) = avg - enddo ; enddo - endif + if (associated(fluxes%salt_flux_added)) & + call homogenize_field_t(fluxes%salt_flux_added, G) - if (associated(fluxes%p_surf_full)) then - avg = global_area_mean(fluxes%p_surf_full, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%p_surf_full(i,j) = avg - enddo ; enddo - endif + if (associated(fluxes%p_surf_full)) & + call homogenize_field_t(fluxes%p_surf_full, G) - if (associated(fluxes%buoy)) then - avg = global_area_mean(fluxes%buoy, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%buoy(i,j) = avg - enddo ; enddo - endif + if (associated(fluxes%buoy)) & + call homogenize_field_t(fluxes%buoy, G) - if (associated(fluxes%TKE_tidal)) then - avg = global_area_mean(fluxes%TKE_tidal, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%TKE_tidal(i,j) = avg - enddo ; enddo - endif + if (associated(fluxes%TKE_tidal)) & + call homogenize_field_t(fluxes%TKE_tidal, G) - if (associated(fluxes%ustar_tidal)) then - avg = global_area_mean(fluxes%ustar_tidal, G) - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.) fluxes%ustar_tidal(i,j) = avg - enddo ; enddo - endif + if (associated(fluxes%ustar_tidal)) & + call homogenize_field_t(fluxes%ustar_tidal, G) ! TODO: tracer flux homogenization ! Having a warning causes a lot of errors (each time step). @@ -3911,6 +3676,50 @@ subroutine homogenize_forcing(fluxes, G) end subroutine homogenize_forcing +subroutine homogenize_field_t(var, G) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + real, dimension(SZI_(G), SZJ_(G)), intent(inout) :: var !< The variable to homogenize + + real :: avg + integer :: i, j, is, ie, js, je + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + + avg = global_area_mean(var, G) + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.) var(i,j) = avg + enddo ; enddo + +end subroutine homogenize_field_t + +subroutine homogenize_field_v(var, G) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + real, dimension(SZI_(G), SZJB_(G)), intent(inout) :: var !< The variable to homogenize + + real :: avg + integer :: i, j, is, ie, jsB, jeB + is = G%isc ; ie = G%iec ; jsB = G%jscB ; jeB = G%jecB + + avg = global_area_mean_v(var, G) + do J=jsB,jeB ; do i=is,ie + if (G%mask2dCv(i,J) > 0.) var(i,J) = avg + enddo ; enddo + +end subroutine homogenize_field_v + +subroutine homogenize_field_u(var, G) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + real, dimension(SZI_(G), SZJB_(G)), intent(inout) :: var !< The variable to homogenize + + real :: avg + integer :: i, j, isB, ieB, js, je + isB = G%iscB ; ieB = G%iecB ; js = G%jsc ; je = G%jec + + avg = global_area_mean_u(var, G) + do j=js,je ; do I=isB,ieB + if (G%mask2dCu(I,j) > 0.) var(I,j) = avg + enddo ; enddo + +end subroutine homogenize_field_u !> \namespace mom_forcing_type !! diff --git a/src/diagnostics/MOM_spatial_means.F90 b/src/diagnostics/MOM_spatial_means.F90 index f68cfe1565..52e7a43355 100644 --- a/src/diagnostics/MOM_spatial_means.F90 +++ b/src/diagnostics/MOM_spatial_means.F90 @@ -48,50 +48,48 @@ function global_area_mean(var, G, scale) end function global_area_mean !> Return the global area mean of a variable. This uses reproducing sums. -function global_area_mean_v(var, G, scale) +function global_area_mean_v(var, G) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure real, dimension(SZI_(G), SZJB_(G)), intent(in) :: var !< The variable to average - real, optional, intent(in) :: scale !< A rescaling factor for the variable - real, dimension(SZI_(G), SZJB_(G)) :: tmpForSumming + real, dimension(SZI_(G), SZJ_(G)) :: tmpForSumming real :: global_area_mean_v - - real :: scalefac ! An overall scaling factor for the areas and variable. integer :: i, j, is, ie, js, je, isB, ieB, jsB, jeB + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec isB = G%iscB ; ieB = G%iecB ; jsB = G%jscB ; jeB = G%jecB - scalefac = G%US%L_to_m**2 ; if (present(scale)) scalefac = G%US%L_to_m**2*scale - tmpForSumming(:,:) = 0. - do J=jsB,jeB ; do i=is,ie - tmpForSumming(i,J) = var(i,J) * (scalefac * G%areaCv(i,J) * G%mask2dCv(i,J)) + do J=js,je ; do i=is,ie + tmpForSumming(i,j) = G%areaT(i,j)/max(1.e-20,& + G%mask2dCv(i,J)+G%mask2dCv(i,J-1)) * & + (var(i,J) * G%mask2dCv(i,J) + & + var(i,J-1) * G%mask2dCv(i,J-1)) enddo ; enddo - global_area_mean_v = reproducing_sum(tmpForSumming) * G%IareaT_global !Need to add IareaCv? + global_area_mean_v = reproducing_sum(tmpForSumming) * G%IareaT_global end function global_area_mean_v !> Return the global area mean of a variable on U grid. This uses reproducing sums. -function global_area_mean_u(var, G, scale) +function global_area_mean_u(var, G) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure real, dimension(SZIB_(G), SZJ_(G)), intent(in) :: var !< The variable to average - real, optional, intent(in) :: scale !< A rescaling factor for the variable real, dimension(SZI_(G), SZJ_(G)) :: tmpForSumming real :: global_area_mean_u - - real :: scalefac ! An overall scaling factor for the areas and variable. integer :: i, j, is, ie, js, je, isB, ieB, jsB, jeB + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec isB = G%iscB ; ieB = G%iecB ; jsB = G%jscB ; jeB = G%jecB - scalefac = G%US%L_to_m**2 ; if (present(scale)) scalefac = G%US%L_to_m**2*scale - tmpForSumming(:,:) = 0. - do j=js,je ; do I=isB,ieB - tmpForSumming(I,j) = var(I,j) * (scalefac * G%areaCu(I,j) * G%mask2dCu(I,j)) + do j=js,je ; do i=is,ie + tmpForSumming(i,j) = G%areaT(i,j)/max(1.e-20,& + G%mask2dCu(I,j)+G%mask2dCu(I-1,j)) * & + (var(I,j) * G%mask2dCu(I,j) + & + var(I-1,j) * G%mask2dCu(I-1,j)) enddo ; enddo - global_area_mean_u = reproducing_sum(tmpForSumming) * G%IareaT_global !Need to add IareaCu? + global_area_mean_u = reproducing_sum(tmpForSumming) * G%IareaT_global end function global_area_mean_u From 7e3f66aafb10378998eab7286d99be102696a4a7 Mon Sep 17 00:00:00 2001 From: "brandon.reichl" Date: Tue, 9 Nov 2021 17:17:15 -0500 Subject: [PATCH 4/4] Minor formatting updates - Move a division and reformat alignment in MOM_spatial_means.F90. - Remove a unused parameter in MOM_forcing_type.F90 - Reformat misalignment of an "if-block" in MOM_forcing_type.F90 --- src/core/MOM_forcing_type.F90 | 3 +-- src/diagnostics/MOM_spatial_means.F90 | 14 ++++++-------- 2 files changed, 7 insertions(+), 10 deletions(-) diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 44d5d8c135..237bb68cc4 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -184,7 +184,6 @@ module MOM_forcing_type !! average of the gustless wind stress. real :: C_p !< heat capacity of seawater [Q degC-1 ~> J kg-1 degC-1]. !! C_p is is the same value as in thermovar_ptrs_type. - real :: gust_const !< Constant unresolved background gustiness for ustar [R L Z T-1 ~> Pa] ! CFC-related arrays needed in the MOM_CFC_cap module real, pointer, dimension(:,:) :: & @@ -3532,7 +3531,7 @@ subroutine homogenize_mech_forcing(forces, G, US, Rho0, UpdateUstar) do j=js,je ; do i=is,ie if (G%mask2dT(i,j) > 0.) forces%ustar(i,j) = sqrt(sqrt(tx_mean**2 + ty_mean**2)*iRho0) enddo ; enddo - else + else call homogenize_field_t(forces%ustar, G) endif else diff --git a/src/diagnostics/MOM_spatial_means.F90 b/src/diagnostics/MOM_spatial_means.F90 index 52e7a43355..182779437c 100644 --- a/src/diagnostics/MOM_spatial_means.F90 +++ b/src/diagnostics/MOM_spatial_means.F90 @@ -61,10 +61,9 @@ function global_area_mean_v(var, G) tmpForSumming(:,:) = 0. do J=js,je ; do i=is,ie - tmpForSumming(i,j) = G%areaT(i,j)/max(1.e-20,& - G%mask2dCv(i,J)+G%mask2dCv(i,J-1)) * & - (var(i,J) * G%mask2dCv(i,J) + & - var(i,J-1) * G%mask2dCv(i,J-1)) + tmpForSumming(i,j) = G%areaT(i,j) * (var(i,J) * G%mask2dCv(i,J) + & + var(i,J-1) * G%mask2dCv(i,J-1)) & + / max(1.e-20,G%mask2dCv(i,J)+G%mask2dCv(i,J-1)) enddo ; enddo global_area_mean_v = reproducing_sum(tmpForSumming) * G%IareaT_global @@ -84,10 +83,9 @@ function global_area_mean_u(var, G) tmpForSumming(:,:) = 0. do j=js,je ; do i=is,ie - tmpForSumming(i,j) = G%areaT(i,j)/max(1.e-20,& - G%mask2dCu(I,j)+G%mask2dCu(I-1,j)) * & - (var(I,j) * G%mask2dCu(I,j) + & - var(I-1,j) * G%mask2dCu(I-1,j)) + tmpForSumming(i,j) = G%areaT(i,j) * (var(I,j) * G%mask2dCu(I,j) + & + var(I-1,j) * G%mask2dCv(I-1,j)) & + / max(1.e-20,G%mask2dCv(I,j)+G%mask2dCv(I-1,j)) enddo ; enddo global_area_mean_u = reproducing_sum(tmpForSumming) * G%IareaT_global