From ed51bcd6678cb6b07e031656d7955504a7e0d667 Mon Sep 17 00:00:00 2001 From: Andrew Shao Date: Fri, 6 Mar 2020 12:39:40 -0500 Subject: [PATCH] Modify remapping of vertically extensive diagnostics Vertically extensive diagnostics like advective transports, diffusive fluxes, etc. were being remapped onto a diagnostic grid prior to the physical process that chagne the tracer and/or layer thicknesses. However this means that the effective operator used for each component of a tracer budget was not the same, and so budgets could never be closed cell-by-cell in a diagnostic coordinate. This commit fixes one part of the inconsistency by setting the target diagnostic grid for all vertically extensive quantities to be one constructed at the very beginning of the timestep. However, problems with closing the budget should still be expected 1. Conservation of column integrals cannot be expected because the target grid can be smaller than the source grid (layer thicknesses at the native grid at the current point in the algorithm). This leads to some fluxes being 'thrown away' due to the reintegrate algorithm. 2. To truly have a consistent operator for all terms of a budget, the source grid for all vertically extensive grids should also be the same --- src/core/MOM.F90 | 4 +++ src/framework/MOM_diag_mediator.F90 | 49 ++++++++++++++++++++++------- src/framework/MOM_diag_remap.F90 | 23 ++++++++------ 3 files changed, 54 insertions(+), 22 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index adb916298f..9e57fc2844 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -622,6 +622,10 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_int_in, CS, & if (showCallTree) call callTree_enter("DT cycles (step_MOM) n=",n) + ! Update the vertically extensive diagnostic grids so that they are + ! referenced to the beginning timestep + call diag_update_remap_grids(CS%diag, update_intensive = .false., update_extensive = .true. ) + !=========================================================================== ! This is the first place where the diabatic processes and remapping could occur. if (CS%diabatic_first .and. (CS%t_dyn_rel_adv==0.0) .and. do_thermo) then ! do thermodynamics. diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 1fc012b7b9..a7fe44a93a 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -1473,14 +1473,16 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask, alt_h) logical :: staggered_in_x, staggered_in_y real, dimension(:,:,:), pointer :: h_diag => NULL() + if (id_clock_diag_mediator>0) call cpu_clock_begin(id_clock_diag_mediator) + + ! For intensive variables only, we can choose to use a different diagnostic grid + ! to map to if (present(alt_h)) then h_diag => alt_h else h_diag => diag_cs%h endif - if (id_clock_diag_mediator>0) call cpu_clock_begin(id_clock_diag_mediator) - ! Iterate over list of diag 'variants', e.g. CMOR aliases, different vertical ! grids, and post each. call assert(diag_field_id < diag_cs%next_free_diag_id, & @@ -1500,10 +1502,11 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask, alt_h) if (id_clock_diag_remap>0) call cpu_clock_begin(id_clock_diag_remap) allocate(remapped_field(size(field,1), size(field,2), diag%axes%nz)) - call vertically_reintegrate_diag_field( & - diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number), & - diag_cs%G, h_diag, staggered_in_x, staggered_in_y, & - diag%axes%mask3d, diag_cs%missing_value, field, remapped_field) + call vertically_reintegrate_diag_field( & + diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number), diag_cs%G, & + diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number)%h_extensive, & + staggered_in_x, staggered_in_y, diag%axes%mask3d, diag_cs%missing_value, & + field, remapped_field) if (id_clock_diag_remap>0) call cpu_clock_end(id_clock_diag_remap) if (associated(diag%axes%mask3d)) then ! Since 3d masks do not vary in the vertical, just use as much as is @@ -3191,7 +3194,7 @@ subroutine diag_set_state_ptrs(h, T, S, eqn_of_state, diag_cs) !> Build/update vertical grids for diagnostic remapping. !! \note The target grids need to be updated whenever sea surface !! height changes. -subroutine diag_update_remap_grids(diag_cs, alt_h, alt_T, alt_S) +subroutine diag_update_remap_grids(diag_cs, alt_h, alt_T, alt_S, update_intensive, update_extensive ) type(diag_ctrl), intent(inout) :: diag_cs !< Diagnostics control structure real, target, optional, intent(in ) :: alt_h(:,:,:) !< Used if remapped grids should be something other than !! the current thicknesses [H ~> m or kg m-2] @@ -3199,11 +3202,17 @@ subroutine diag_update_remap_grids(diag_cs, alt_h, alt_T, alt_S) !! the current temperatures real, target, optional, intent(in ) :: alt_S(:,:,:) !< Used if remapped grids should be something other than !! the current salinity + logical, optional, intent(in ) :: update_intensive !< If true (default), update the grids used for + !! intensive diagnostics + logical, optional, intent(in ) :: update_extensive !< If true (not default), update the grids used for + !! intensive diagnostics ! Local variables integer :: i real, dimension(:,:,:), pointer :: h_diag => NULL() real, dimension(:,:,:), pointer :: T_diag => NULL(), S_diag => NULL() + logical :: update_intensive_local, update_extensive_local + ! Set values based on optional input arguments if (present(alt_h)) then h_diag => alt_h else @@ -3222,6 +3231,15 @@ subroutine diag_update_remap_grids(diag_cs, alt_h, alt_T, alt_S) S_diag => diag_CS%S endif + ! Defaults here are based on wanting to update intensive quantities frequently as soon as the model state changes. + ! Conversely, for extensive quantities, in an effort to close budgets and to be consistent with the total time + ! tendency, we construct the diagnostic grid at the beginning of the baroclinic timestep and remap all extensive + ! quantities to the same grid + update_intensive_local = .true. + if (present(update_intensive)) update_intensive_local = update_intensive + update_extensive_local = .false. + if (present(update_extensive)) update_extensive_local = update_extensive + if (id_clock_diag_grid_updates>0) call cpu_clock_begin(id_clock_diag_grid_updates) if (diag_cs%diag_grid_overridden) then @@ -3229,11 +3247,18 @@ subroutine diag_update_remap_grids(diag_cs, alt_h, alt_T, alt_S) "diagnostic structure have been overridden") endif - do i=1, diag_cs%num_diag_coords - call diag_remap_update(diag_cs%diag_remap_cs(i), & - diag_cs%G, diag_cs%GV, diag_cs%US, h_diag, T_diag, S_diag, & - diag_cs%eqn_of_state) - enddo + if (update_intensive_local) then + do i=1, diag_cs%num_diag_coords + call diag_remap_update(diag_cs%diag_remap_cs(i), diag_cs%G, diag_cs%GV, diag_cs%US, h_diag, T_diag, S_diag, & + diag_cs%eqn_of_state, diag_cs%diag_remap_cs(i)%h) + enddo + endif + if (update_extensive_local) then + do i=1, diag_cs%num_diag_coords + call diag_remap_update(diag_cs%diag_remap_cs(i), diag_cs%G, diag_cs%GV, diag_cs%US, h_diag, T_diag, S_diag, & + diag_cs%eqn_of_state, diag_cs%diag_remap_cs(i)%h_extensive) + enddo + endif #if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__) ! Keep a copy of H - used to check whether grids are up-to-date diff --git a/src/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index b61c10eb7e..81e7187786 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -112,6 +112,7 @@ module MOM_diag_remap type(regridding_CS) :: regrid_cs !< Regridding control structure that defines the coordiantes for this axes integer :: nz = 0 !< Number of vertical levels used for remapping real, dimension(:,:,:), allocatable :: h !< Remap grid thicknesses + real, dimension(:,:,:), allocatable :: h_extensive !< Remap grid thicknesses for extensive variables real, dimension(:), allocatable :: dz !< Nominal layer thicknesses integer :: interface_axes_id = 0 !< Vertical axes id for remapping at interfaces integer :: layer_axes_id = 0 !< Vertical axes id for remapping on layers @@ -271,15 +272,16 @@ function diag_remap_axes_configured(remap_cs) !! height or layer thicknesses changes. In the case of density-based !! coordinates then technically we should also regenerate the !! target grid whenever T/S change. -subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state) - type(diag_remap_ctrl), intent(inout) :: remap_cs !< Diagnostic coordinate control structure - type(ocean_grid_type), pointer :: G !< The ocean's grid type - type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - real, dimension(:, :, :), intent(in) :: h !< New thickness - real, dimension(:, :, :), intent(in) :: T !< New T - real, dimension(:, :, :), intent(in) :: S !< New S +subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state, h_target) + type(diag_remap_ctrl), intent(inout) :: remap_cs !< Diagnostic coordinate control structure + type(ocean_grid_type), pointer :: G !< The ocean's grid type + type(verticalGrid_type), intent(in ) :: GV !< ocean vertical grid structure + type(unit_scale_type), intent(in ) :: US !< A dimensional unit scaling type + real, dimension(:, :, :), intent(in ) :: h !< Thicknesses used to construct new diagnostic grid + real, dimension(:, :, :), intent(in ) :: T !< Temperatures used to construct new diagnostic grid + real, dimension(:, :, :), intent(in ) :: S !< Salinity used to construct new diagnostic grid type(EOS_type), pointer :: eqn_of_state !< A pointer to the equation of state + real, dimension(:, :, :), intent(inout) :: h_target !< Where to store the new diagnostic array ! Local variables real, dimension(remap_cs%nz + 1) :: zInterfaces @@ -306,6 +308,7 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state) call initialize_remapping(remap_cs%remap_cs, 'PPM_IH4', boundary_extrapolation=.false., & answers_2018=remap_cs%answers_2018) allocate(remap_cs%h(G%isd:G%ied,G%jsd:G%jed, nz)) + allocate(remap_cs%h_extensive(G%isd:G%ied,G%jsd:G%jed, nz)) remap_cs%initialized = .true. endif @@ -314,7 +317,7 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state) ! assumption that h, T, S has changed. do j=G%jsc-1, G%jec+1 ; do i=G%isc-1, G%iec+1 if (G%mask2dT(i,j)==0.) then - remap_cs%h(i,j,:) = 0. + h_target(i,j,:) = 0. cycle endif @@ -338,7 +341,7 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state) ! US%Z_to_m*G%bathyT(i,j), sum(h(i,j,:)), zInterfaces) call MOM_error(FATAL,"diag_remap_update: HYCOM1 coordinate not coded for diagnostics yet!") endif - remap_cs%h(i,j,:) = zInterfaces(1:nz) - zInterfaces(2:nz+1) + h_target(i,j,:) = zInterfaces(1:nz) - zInterfaces(2:nz+1) enddo ; enddo end subroutine diag_remap_update