diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 5c7f79fe32..80cc36c577 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -660,6 +660,10 @@ subroutine step_MOM(forces_in, fluxes_in, 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 ceb782ce4b..f9e0bf0fbd 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -333,6 +333,9 @@ module MOM_diag_mediator !> Number of checksum-only diagnostics integer :: num_chksum_diags + real, dimension(:,:,:), allocatable :: h_begin !< Layer thicknesses at the beginning of the timestep used + !! for remapping of extensive variables + end type diag_ctrl !>@{ CPU clocks @@ -456,6 +459,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 @@ -1475,14 +1482,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, & @@ -1502,10 +1511,12 @@ 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%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) 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 @@ -3064,6 +3075,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 @@ -3192,7 +3204,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] @@ -3200,11 +3212,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() ! The layer thickneses for diagnostics [H ~> m or kg m-2] 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 @@ -3223,6 +3241,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 @@ -3230,10 +3257,19 @@ 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 + 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) + 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 cadd74950a..fecf17b84d 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,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) +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 @@ -280,6 +281,7 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state) real, dimension(:, :, :), intent(in) :: T !< New temperatures [degC] real, dimension(:, :, :), intent(in) :: S !< New salinities [ppt] 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 ! Interface positions [H ~> m or kg m-2] @@ -305,7 +307,6 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state) ! 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)) remap_cs%initialized = .true. endif @@ -314,7 +315,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 @@ -339,7 +340,9 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state) ! GV%Z_to_H*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) + do k = 1,nz + h_target(i,j,k) = zInterfaces(k) - zInterfaces(k+1) + enddo enddo ; enddo end subroutine diag_remap_update @@ -485,11 +488,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_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 current thicknesses + 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 real, dimension(:,:,:), pointer :: mask !< A mask for the field @@ -526,7 +530,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 @@ -541,7 +545,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 @@ -554,7 +558,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