diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index d4224a14cf..71a2e2cb51 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -143,21 +143,22 @@ module MOM_diag_mediator type(remapping_CS) :: remap_cs type(regridding_CS) :: regrid_cs ! Nominal interface locations for Z remapping - real, dimension(:), allocatable :: zi_nom + real, dimension(:), allocatable :: zi_remap ! Nominal layer locations for Z remapping - real, dimension(:), allocatable :: zl_nom + real, dimension(:), allocatable :: zl_remap ! Number of z levels used for remapping integer :: nz_remap - ! Define z star on u, v, T grids in terms of H. These are used for diagnostic remapping. - real, dimension(:,:,:), allocatable :: h_uz - real, dimension(:,:,:), allocatable :: h_vz - real, dimension(:,:,:), allocatable :: h_Tz + ! Define z star on u, v, T grids, these are the interface positions + real, dimension(:,:,:), allocatable :: zi_u + real, dimension(:,:,:), allocatable :: zi_v + real, dimension(:,:,:), allocatable :: zi_T + ! Keep track of which remapping is needed for diagnostic output logical :: do_z_remapping_on_u, do_z_remapping_on_v, do_z_remapping_on_T logical :: remapping_initialized - ! Pointer to H and G for Z remapping + ! Pointer to H and G for remapping real, dimension(:,:,:), pointer :: h => null() type(ocean_grid_type), pointer :: G => null() @@ -341,74 +342,79 @@ end subroutine set_axes_info subroutine update_diag_target_grids(G, diag_cs) ! Build target grids for diagnostic remapping + type(ocean_grid_type), intent(inout) :: G type(diag_ctrl), intent(inout) :: diag_cs + ! Arguments: ! (inout) G - ocean grid structure. ! (inout) diag_cs - structure used to regulate diagnostic output. - real, dimension(G%nk) :: h_src - real, dimension(G%iscB:G%iecB, G%jsc:G%jec, diag_cs%nz_remap) :: z_ui - real, dimension(G%jsc:G%jec, G%iscB:G%iecB, diag_cs%nz_remap) :: z_vi - real, dimension(G%jsc:G%jec, G%jsc:G%jec, diag_cs%nz_remap) :: z_Ti + real, dimension(size(diag_cs%h, 3)) :: h_src real :: depth + integer :: nz_src, nz_dest + integer :: i, j, k + + nz_dest = diag_cs%nz_remap + nz_src = size(diag_cs%h, 3) if ((diag_cs%do_z_remapping_on_u .or. diag_cs%do_z_remapping_on_v .or. & diag_cs%do_z_remapping_on_T) .and. (.not. diag_cs%remapping_initialized)) then + call assert(allocated(diag_cs%zi_remap), 'Remapping axis not initialized') + ! Initialise remapping system - call initialize_regridding(diag_cs%nz_remap, 'Z*', 'PPM_IH4', diag_cs%regrid_cs) - call initialize_remapping(G%nk, 'PPM_IH4', diag%remap_cs) + call initialize_regridding(nz_dest, 'Z*', 'PPM_IH4', diag_cs%regrid_cs) + call initialize_remapping(nz_src, 'PPM_IH4', diag_cs%remap_cs) call setRegriddingMinimumThickness(G%Angstrom, diag_cs%regrid_cs) - call setCoordinateResolution(diag_cs%zi_remap(2:) - diag_cs%zi_remap(:diag_cs%nz_remap), diag_cs%regrid_cs) + call setCoordinateResolution(diag_cs%zi_remap(2:) - diag_cs%zi_remap(:nz_dest), diag_cs%regrid_cs) diag_cs%remapping_initialized = .true. endif ! Build z-star grid on u points if (diag_cs%do_z_remapping_on_u) then - if (.not. associated(diag%z_ui)) then - allocate(diag%h_uz(G%jsc:G%jse,G%iscB,G%iecB)) + if (.not. allocated(diag_cs%zi_u)) then + allocate(diag_cs%zi_u(G%iscB:G%iecB,G%jsc:G%jec,nz_dest+1)) endif do j=G%jsc, G%jec do i=G%iscB, G%iecB - h_src(:) = 0.5 * (h(i,j,:) + h(i+1,j,:)) + h_src(:) = 0.5 * (diag_cs%h(i,j,:) + diag_cs%h(i+1,j,:)) depth = 0.5 * (G%bathyT(i,j) + G%bathyT(i+1,j)) - call buildGridZstarColumn(regrid_cs, diag_cs%nz_remap, depth, sum(h_src(:)), z_ui(i, j, :)) - z_ui(i, j, :) = -z_ui(i, j, :) - diag_cs%h_uz(:) = z_ui(2:) - z_ui(:size(z_ui, 3)-1) + call buildGridZstarColumn(diag_cs%regrid_cs, nz_dest, depth, sum(h_src(:)), diag_cs%zi_u(i, j, :)) + diag_cs%zi_u(i, j, :) = -diag_cs%zi_u(i, j, :) enddo enddo endif ! Build z-star grid on u points if (diag_cs%do_z_remapping_on_v) then - if (.not. associated(diag%z_ui)) then - allocate(diag%h_vz(G%jscB:G%jseB,G%isc,G%iec)) + if (.not. allocated(diag_cs%zi_v)) then + allocate(diag_cs%zi_v(G%isc:G%iec,G%jscB:G%jecB,nz_dest+1)) endif + do j=G%jscB, G%jecB do i=G%isc, G%iec h_src(:) = 0.5 * (diag_cs%h(i,j,:) + diag_cs%h(i,j+1,:)) depth = 0.5 * (G%bathyT(i, j) + G%bathyT(i, j+1)) - call buildGridZstarColumn(regrid_cs, diag_cs%nz_remap, depth, total_thickness, diag%z_vi(i, j, :)) - z_vi(i, j, :) = -z_vi(i, j, :) - diag_cs%h_vz(:) = z_vi(2:) - z_vi(:size(z_vi, 3)-1) + call buildGridZstarColumn(diag_cs%regrid_cs, nz_dest, depth, sum(h_src(:)), diag_cs%zi_v(i, j, :)) + diag_cs%zi_v(i, j, :) = -diag_cs%zi_v(i, j, :) enddo enddo endif ! Build z-star grid on T points if (diag_cs%do_z_remapping_on_T) then - if (.not. associated(diag%z_Ti)) then - allocate(diag%h_Tz(G%jsc:G%jse,G%isc,G%iec)) + if (.not. allocated(diag_cs%zi_T)) then + allocate(diag_cs%zi_T(G%isc:G%iec,G%jsc:G%jec,nz_dest+1)) endif + do j=G%jsc, G%jec do i=G%isc, G%iec - call buildGridZstarColumn(regrid_cs, diag_cs%nz_remap, G%bathyT(i, j), sum(h(i, j, :)), diag%z_Ti(i, j, :)) - z_Ti(i, j, :) = -z_Ti(i, j, :) - diag_cs%h_Tz(:) = z_Ti(2:) - z_Ti(:size(z_Ti, 3)-1) + call buildGridZstarColumn(diag_cs%regrid_cs, nz_dest, G%bathyT(i, j), sum(diag_cs%h(i, j, :)), diag_cs%zi_T(i, j, :)) + diag_cs%zi_T(i, j, :) = -diag_cs%zi_T(i, j, :) enddo enddo endif -end subroutine build_diag_target_grids() +end subroutine update_diag_target_grids function check_grid_def(filename, varname) ! Do some basic checks on the vertical grid definition file, variable @@ -723,8 +729,9 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) ! (out) remapped_field - the field argument remapped to z star coordinate real, dimension(diag_cs%nz_remap+1) :: dz - integer :: nz_src - integer :: i, j + real, dimension(diag_cs%nz_remap) :: h_dest + integer :: nz_src, nz_dest + integer :: i, j, k call assert(size(field, 3) == size(diag_cs%h, 3), & 'Remap field and thickness z-axes do not match.') @@ -734,72 +741,80 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) nz_dest = diag_cs%nz_remap if (is_u_axes(diag%remap_axes, diag_cs)) then - call assert(allocated(diag_cs%zi_remap), 'Remapping axis not initialized') - do j=diag_cs%G%jsc, diag_cs%G%jec do i=diag_cs%G%iscB, diag_cs%G%iecB if (associated(diag%mask3d)) then if (diag%mask3d(i,j, 1) == 0.) cycle endif - call dzFromH1H2(nz_src, h(:), nz_dest, diag_cs%h_uz, dz) - call remapping_core(remap_cs, nz_src, h(:), field(:), nz_dest, dz, & - remapped_field(:)) + h_dest(:) = diag_cs%zi_u(i, j, 2:) - diag_cs%zi_u(i, j, :size(diag_cs%zi_u, 3)-1) + call dzFromH1H2(nz_src, diag_cs%h(i, j, :), nz_dest, h_dest(:), dz) + call remapping_core(diag_cs%remap_cs, nz_src, diag_cs%h(i, j, :), field(i, j, :), nz_dest, dz, & + remapped_field(i, j, :)) + + ! Lower levels of the remapped data get squashed to follow bathymetry, + ! their depth does not corrospond to the nominal depth at that level + ! (and the nominal depth is below the bathymetry), so remove + do k=1, nz_dest + if (diag_cs%zi_u(i, j, k) >= diag_cs%G%bathyT(i, j)) then + remapped_field(i, j, k:nz_dest) = diag_cs%missing_value + exit + endif + enddo enddo enddo - elseif (is_v_axes(diag%remap_axes, diag_cs)) then - call assert(allocated(diag_cs%zi_remap), 'Remapping axis not initialized') + elseif (is_v_axes(diag%remap_axes, diag_cs)) then do j=diag_cs%G%jscB, diag_cs%G%jecB do i=diag_cs%G%isc, diag_cs%G%iec if (associated(diag%mask3d)) then if (diag%mask3d(i,j, 1) == 0.) cycle endif - call dzFromH1H2(nz_src, h(:), nz_dest, diag_cs%h_vz, dz) - call remapping_core(remap_cs, nz_src, h(:), field(:), nz_dest, dz, & - remapped_field(:)) - enddo - enddo - elseif (is_B_axes(diag%remap_axes, diag_cs)) then - call assert(allocated(diag_cs%zi_remap), 'Remapping axis not initialized') - - do j=diag_cs%G%jscB, diag_cs%G%jecB - do i=diag_cs%G%iscB, diag_cs%G%iecB - if (associated(diag%mask3d)) then - if (diag%mask3d(i,j, 1) == 0.) cycle - endif - call dzFromH1H2(nz_src, h(:), nz_dest, diag_cs%h_Bz, dz) - call remapping_core(remap_cs, nz_src, h(:), field(:), nz_dest, dz, & - remapped_field(:)) + h_dest(:) = diag_cs%zi_v(i, j, 2:) - diag_cs%zi_v(i, j, :size(diag_cs%zi_v, 3)-1) + call dzFromH1H2(nz_src, diag_cs%h(i, j, :), nz_dest, h_dest(:), dz) + call remapping_core(diag_cs%remap_cs, nz_src, diag_cs%h(i, j, :), field(i, j, :), nz_dest, dz, & + remapped_field(i, j, :)) + do k=1, nz_dest + if (diag_cs%zi_v(i, j, k) >= diag_cs%G%bathyT(i, j)) then + remapped_field(i, j, k:nz_dest) = diag_cs%missing_value + exit + endif + enddo enddo enddo + !elseif (is_B_axes(diag%remap_axes, diag_cs)) then + ! do j=diag_cs%G%jscB, diag_cs%G%jecB + ! do i=diag_cs%G%iscB, diag_cs%G%iecB + ! if (associated(diag%mask3d)) then + ! if (diag%mask3d(i,j, 1) == 0.) cycle + ! endif + ! h_dest(:) = diag_cs%zi_B(i, j, 2:) - diag_cs%zi_B(i, j, :size(diag_cs%zi_B, 3)-1) + ! call dzFromH1H2(nz_src, diag_cs%h(i, j, :), nz_dest, h_dest(:), dz) + ! call remapping_core(diag_cs%remap_cs, nz_src, diag_cs%h(i, j, :), field(i, j, :), nz_dest, dz, & + ! remapped_field(i, j, :)) + ! enddo + ! enddo else - call assert(allocated(diag_cs%zi_remap), 'Remapping axis not initialized') - do j=diag_cs%G%jsc, diag_cs%G%jec do i=diag_cs%G%isc, diag_cs%G%iec if (associated(diag%mask3d)) then if (diag%mask3d(i,j, 1) == 0.) cycle endif - h_dest(:) = zi_remap(2:) - zi_remap(:size(zi_remap)-1) - call dzFromH1H2(nz_src, h(:), nz_dest, diag_cs%h_Tz, dz) - call remapping_core(remap_cs, nz_src, h(:), field(:), nz_dest, dz, & - remapped_field(:)) + h_dest(:) = diag_cs%zi_T(i, j, 2:) - diag_cs%zi_T(i, j, :size(diag_cs%zi_T, 3)-1) + call dzFromH1H2(nz_src, diag_cs%h(i, j, :), nz_dest, h_dest(:), dz) + call remapping_core(diag_cs%remap_cs, nz_src, diag_cs%h(i, j, :), field(i, j, :), nz_dest, dz, & + remapped_field(i, j, :)) + do k=1, nz_dest + if (diag_cs%zi_T(i, j, k) >= diag_cs%G%bathyT(i, j)) then + remapped_field(i, j, k:nz_dest) = diag_cs%missing_value + exit + endif + enddo enddo enddo endif - ! Lower levels of the remapped data get squashed to follow bathymetry, - ! their depth does not corrospond to the nominal depth at that level - ! (and the nominal depth is below the bathymetry), so remove - do k=1, nz_dest - if (zi_remap(k) >= depth) then - remapped_field(k:nz_dest) = missing_value - exit - endif - enddo - end subroutine remap_diag_to_z subroutine post_data_3d_low(diag, field, diag_cs, is_static, mask) @@ -1452,9 +1467,9 @@ subroutine diag_mediator_init(G, param_file, diag_cs, err_msg) ! Keep a pointer to the grid, this is needed for regridding diag_cs%G => G diag_cs%nz_remap = -1 - do_z_remapping_on_u = .false. - do_z_remapping_on_v = .false. - do_z_remapping_on_T = .false. + diag_cs%do_z_remapping_on_u = .false. + diag_cs%do_z_remapping_on_v = .false. + diag_cs%do_z_remapping_on_T = .false. diag_cs%is = G%isc - (G%isd-1) ; diag_cs%ie = G%iec - (G%isd-1) diag_cs%js = G%jsc - (G%jsd-1) ; diag_cs%je = G%jec - (G%jsd-1)