From ed51bcd6678cb6b07e031656d7955504a7e0d667 Mon Sep 17 00:00:00 2001 From: Andrew Shao Date: Fri, 6 Mar 2020 12:39:40 -0500 Subject: [PATCH 1/8] 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 From 9abcb103658921e65873ee2cedb15f668694e9a4 Mon Sep 17 00:00:00 2001 From: Andrew Shao Date: Wed, 11 Mar 2020 14:51:20 -0400 Subject: [PATCH 2/8] Set target diagnostic grid based on boolean The target grid for the diagnostic grid update is now based on a assigning a differrent pointer based on the boolean input argument 'intensive'. This is done so that this update is done in a more 'object-oriented' way --- src/framework/MOM_diag_mediator.F90 | 4 ++-- src/framework/MOM_diag_remap.F90 | 12 ++++++++++-- 2 files changed, 12 insertions(+), 4 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index a7fe44a93a..e6efe12b1b 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -3250,13 +3250,13 @@ subroutine diag_update_remap_grids(diag_cs, alt_h, alt_T, alt_S, update_intensiv 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) + diag_cs%eqn_of_state, intensive=.true.) 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) + diag_cs%eqn_of_state, intensive=.false.) enddo endif diff --git a/src/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index 81e7187786..914a815387 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -272,7 +272,7 @@ 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, h_target) +subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state, intensive) 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 @@ -281,7 +281,9 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state, h_targe 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 + logical, intent(in ) :: intensive !< If true, update the intensive diagnostic array + + real, dimension(:,:,:), pointer :: h_target !< Where to store the new diagnostic array ! Local variables real, dimension(remap_cs%nz + 1) :: zInterfaces @@ -312,6 +314,12 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state, h_targe remap_cs%initialized = .true. endif + if (intensive) then + h_target => remap_cs%h + else + h_target => reamp_cs%h_extensive + endif + ! Calculate remapping thicknesses for different target grids based on ! nominal/target interface locations. This happens for every call on the ! assumption that h, T, S has changed. From 3477bbe45e123e5e889680491049a1c0cb1afcf7 Mon Sep 17 00:00:00 2001 From: Andrew Shao Date: Wed, 11 Mar 2020 15:01:29 -0400 Subject: [PATCH 3/8] Revert "Set target diagnostic grid based on boolean" This reverts commit 9abcb103658921e65873ee2cedb15f668694e9a4. --- src/framework/MOM_diag_mediator.F90 | 4 ++-- src/framework/MOM_diag_remap.F90 | 12 ++---------- 2 files changed, 4 insertions(+), 12 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index e6efe12b1b..a7fe44a93a 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -3250,13 +3250,13 @@ subroutine diag_update_remap_grids(diag_cs, alt_h, alt_T, alt_S, update_intensiv 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, intensive=.true.) + 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, intensive=.false.) + diag_cs%eqn_of_state, diag_cs%diag_remap_cs(i)%h_extensive) enddo endif diff --git a/src/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index 914a815387..81e7187786 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -272,7 +272,7 @@ 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, intensive) +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 @@ -281,9 +281,7 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state, intensi 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 - logical, intent(in ) :: intensive !< If true, update the intensive diagnostic array - - real, dimension(:,:,:), pointer :: h_target !< Where to store the new diagnostic array + real, dimension(:, :, :), intent(inout) :: h_target !< Where to store the new diagnostic array ! Local variables real, dimension(remap_cs%nz + 1) :: zInterfaces @@ -314,12 +312,6 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state, intensi remap_cs%initialized = .true. endif - if (intensive) then - h_target => remap_cs%h - else - h_target => reamp_cs%h_extensive - endif - ! Calculate remapping thicknesses for different target grids based on ! nominal/target interface locations. This happens for every call on the ! assumption that h, T, S has changed. From 46ce21ec7d4c7e8f7fe65906cbfa43c654b4f22e Mon Sep 17 00:00:00 2001 From: Andrew Shao Date: Fri, 13 Mar 2020 10:59:59 -0700 Subject: [PATCH 4/8] Store h_extensive for both native and remapped coordinates The model thicknesses --- src/framework/MOM_diag_mediator.F90 | 5 +++++ src/framework/MOM_diag_remap.F90 | 13 +++++++------ 2 files changed, 12 insertions(+), 6 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index e6efe12b1b..b4b4d9162d 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -332,6 +332,8 @@ module MOM_diag_mediator !> Number of checksum-only diagnostics integer :: num_chksum_diags + real, dimension(:,:,:), allocatable :: h_begin + end type diag_ctrl ! CPU clocks @@ -1504,6 +1506,7 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask, alt_h) 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, & + diag_cs%h_extensive, 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) @@ -3066,6 +3069,7 @@ subroutine diag_mediator_init(G, GV, US, nz, param_file, diag_cs, doc_file_dir) diag_cs%S => null() diag_cs%eqn_of_state => null() + allocate(diag_cs%h_begin(G%isd:G%ied,G%jsd:G%jed,nz)) #if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__) allocate(diag_cs%h_old(G%isd:G%ied,G%jsd:G%jed,nz)) diag_cs%h_old(:,:,:) = 0.0 @@ -3254,6 +3258,7 @@ subroutine diag_update_remap_grids(diag_cs, alt_h, alt_T, alt_S, update_intensiv enddo endif if (update_extensive_local) then + CS%h_begin(:,:,:) = CS%h(:,:,:) 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, intensive=.false.) diff --git a/src/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index 914a815387..8839c1df70 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -317,7 +317,7 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state, intensi if (intensive) then h_target => remap_cs%h else - h_target => reamp_cs%h_extensive + h_target => remap_cs%h_extensive endif ! Calculate remapping thicknesses for different target grids based on @@ -495,11 +495,12 @@ subroutine diag_remap_calc_hmask(remap_cs, G, mask) end subroutine diag_remap_calc_hmask !> Vertically re-grid an already vertically-integrated diagnostic field to alternative vertical grid. -subroutine vertically_reintegrate_diag_field(remap_cs, G, h, staggered_in_x, staggered_in_y, & +subroutine vertically_reintegrate_diag_field(remap_cs, G, h, h_dest, staggered_in_x, staggered_in_y, & mask, missing_value, field, reintegrated_field) type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coodinate control structure type(ocean_grid_type), intent(in) :: G !< Ocean grid structure - real, dimension(:,:,:), intent(in) :: h !< The current thicknesses + real, dimension(:,:,:), intent(in) :: h !< The thicknesses of the source grid + real, dimension(:,:,:), intent(in) :: h_target !< The thicknesses of the target grid logical, intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points logical, intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points real, dimension(:,:,:), pointer :: mask !< A mask for the field @@ -536,7 +537,7 @@ subroutine vertically_reintegrate_diag_field(remap_cs, G, h, staggered_in_x, sta if (mask(I,j,1) == 0.) cycle endif h_src(:) = 0.5 * (h(i_lo,j,:) + h(i_hi,j,:)) - h_dest(:) = 0.5 * (remap_cs%h(i_lo,j,:) + remap_cs%h(i_hi,j,:)) + h_dest(:) = 0.5 * (h_target(i_lo,j,:) + h_target(i_hi,j,:)) call reintegrate_column(nz_src, h_src, field(I1,j,:), & nz_dest, h_dest, 0., reintegrated_field(I1,j,:)) enddo @@ -551,7 +552,7 @@ subroutine vertically_reintegrate_diag_field(remap_cs, G, h, staggered_in_x, sta if (mask(i,J,1) == 0.) cycle endif h_src(:) = 0.5 * (h(i,j_lo,:) + h(i,j_hi,:)) - h_dest(:) = 0.5 * (remap_cs%h(i,j_lo,:) + remap_cs%h(i,j_hi,:)) + h_dest(:) = 0.5 * (h_target(i,j_lo,:) + h_target(i,j_hi,:)) call reintegrate_column(nz_src, h_src, field(i,J1,:), & nz_dest, h_dest, 0., reintegrated_field(i,J1,:)) enddo @@ -564,7 +565,7 @@ subroutine vertically_reintegrate_diag_field(remap_cs, G, h, staggered_in_x, sta if (mask(i,j,1) == 0.) cycle endif h_src(:) = h(i,j,:) - h_dest(:) = remap_cs%h(i,j,:) + h_dest(:) = h_target(i,j,:) call reintegrate_column(nz_src, h_src, field(i,j,:), & nz_dest, h_dest, 0., reintegrated_field(i,j,:)) enddo From e23f03c811db789b8d7346c437b25b368d9ac1e0 Mon Sep 17 00:00:00 2001 From: Andrew Shao Date: Tue, 17 Mar 2020 14:08:07 -0400 Subject: [PATCH 5/8] Fix compilation errors - Rename `h_dest` to `h_target` in routine and in signature - Remove extraneous logic --- src/framework/MOM_diag_mediator.F90 | 4 ++-- src/framework/MOM_diag_remap.F90 | 12 +++--------- 2 files changed, 5 insertions(+), 11 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index f18b63f82d..45376b628f 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -1506,7 +1506,7 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask, alt_h) 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, & - diag_cs%h_extensive, + diag_cs%h_begin, & 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) @@ -3258,7 +3258,7 @@ subroutine diag_update_remap_grids(diag_cs, alt_h, alt_T, alt_S, update_intensiv enddo endif if (update_extensive_local) then - CS%h_begin(:,:,:) = CS%h(:,:,:) + diag_cs%h_begin(:,:,:) = diag_cs%h(:,:,:) 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) diff --git a/src/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index bb72ad1af0..77aa9efb91 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -312,12 +312,6 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state, h_targe remap_cs%initialized = .true. endif - if (intensive) then - h_target => remap_cs%h - else - h_target => remap_cs%h_extensive - endif - ! Calculate remapping thicknesses for different target grids based on ! nominal/target interface locations. This happens for every call on the ! assumption that h, T, S has changed. @@ -493,11 +487,11 @@ subroutine diag_remap_calc_hmask(remap_cs, G, mask) end subroutine diag_remap_calc_hmask !> Vertically re-grid an already vertically-integrated diagnostic field to alternative vertical grid. -subroutine vertically_reintegrate_diag_field(remap_cs, G, h, h_dest, staggered_in_x, staggered_in_y, & +subroutine vertically_reintegrate_diag_field(remap_cs, G, h, h_target, staggered_in_x, staggered_in_y, & mask, missing_value, field, reintegrated_field) type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coodinate control structure - type(ocean_grid_type), intent(in) :: G !< Ocean grid structure - real, dimension(:,:,:), intent(in) :: h !< The thicknesses of the source grid + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + real, dimension(:,:,:), intent(in) :: h !< The thicknesses of the source grid real, dimension(:,:,:), intent(in) :: h_target !< The thicknesses of the target grid logical, intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points logical, intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points From 9a330fab88ac5fccb4ac6513c5713ce0ee821ba0 Mon Sep 17 00:00:00 2001 From: Andrew Shao Date: Tue, 17 Mar 2020 14:56:38 -0400 Subject: [PATCH 6/8] Loop over recreation of diagnostic grid explicitly --- src/framework/MOM_diag_remap.F90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index 77aa9efb91..000a3ce518 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -341,7 +341,9 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state, h_targe ! 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 - h_target(i,j,:) = zInterfaces(1:nz) - zInterfaces(2:nz+1) + do k = 1,nz + h_target(i,j,k) = zInterfaces(k) - zInterfaces(k+1) + enddo enddo ; enddo end subroutine diag_remap_update From 24782fa85bdb631c5290fb2ad267aab6dee0f0c0 Mon Sep 17 00:00:00 2001 From: Andrew Shao Date: Tue, 17 Mar 2020 18:19:00 -0400 Subject: [PATCH 7/8] Move allocation of diagnostic arrays Originally, the diagnostic arrays were being allocated on the first call generating the diagnostic grid. This seems overly clunky because the size of the grids are known before that. This was also causing a segfault in updates to the new routine. The allocate statements for these arrays are now done right after the number of levels is known --- src/framework/MOM_diag_mediator.F90 | 4 ++++ src/framework/MOM_diag_remap.F90 | 2 -- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 45376b628f..1d0d204354 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -456,6 +456,10 @@ subroutine set_axes_info(G, GV, US, param_file, diag_cs, set_vertical) ! For each possible diagnostic coordinate call diag_remap_configure_axes(diag_cs%diag_remap_cs(i), GV, US, param_file) + ! Allocate these arrays since the size of the diagnostic array is now known + allocate(diag_cs%diag_remap_cs(i)%h(G%isd:G%ied,G%jsd:G%jed, diag_cs%diag_remap_cs(i)%nz)) + allocate(diag_cs%diag_remap_cs(i)%h_extensive(G%isd:G%ied,G%jsd:G%jed, diag_cs%diag_remap_cs(i)%nz)) + ! This vertical coordinate has been configured so can be used. if (diag_remap_axes_configured(diag_cs%diag_remap_cs(i))) then diff --git a/src/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index 000a3ce518..e4f8b410d2 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -307,8 +307,6 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state, h_targe ! Initialize remapping and regridding on the first call 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 From c309289b6b6930dc7dc14491be3da4ee91a3c024 Mon Sep 17 00:00:00 2001 From: Andrew Shao Date: Thu, 16 Apr 2020 19:05:45 -0400 Subject: [PATCH 8/8] Add doxygen comments for h_begin in diag_ctrl --- src/framework/MOM_diag_mediator.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 1d0d204354..8ec8349e58 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -332,7 +332,8 @@ module MOM_diag_mediator !> Number of checksum-only diagnostics integer :: num_chksum_diags - real, dimension(:,:,:), allocatable :: h_begin + real, dimension(:,:,:), allocatable :: h_begin !< Layer thicknesses at the beginning of the timestep used + !! for remapping of extensive variables end type diag_ctrl