From 47bed87a7191ac9371f7d6dd93f58955f7883b08 Mon Sep 17 00:00:00 2001 From: Nicholas Hannah Date: Fri, 16 Dec 2016 13:47:33 +1100 Subject: [PATCH 1/2] Allow horizontally_average_diag_field to be called by both native and remapped diagnostics. #334 --- src/framework/MOM_diag_mediator.F90 | 40 +++++++++++++++------ src/framework/MOM_diag_remap.F90 | 55 ++++++++++++----------------- 2 files changed, 52 insertions(+), 43 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 9e6fd47518..75fd384599 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -61,9 +61,7 @@ module MOM_diag_mediator implicit none ; private #define __DO_SAFETY_CHECKS__ -#define RANGE_I(a) lbound(a, 1),ubound(a, 1) -#define RANGE_J(a) lbound(a, 2),ubound(a, 2) -#define RANGE_K(a) lbound(a, 3),ubound(a, 3) +#define IMPLIES(A, B) ((.not. (A)) .or. (B)) public set_axes_info, post_data, register_diag_field, time_type public post_data_1d_k @@ -463,6 +461,7 @@ subroutine define_axes_group(diag_cs, handles, axes, nz, vertical_coordinate_num axes%rank = n axes%handles(:) = handles(:) axes%diag_cs => diag_cs ! A [circular] link back to the diag_cs structure + if (present(x_cell_method)) then if (axes%rank<2) call MOM_error(FATAL, 'define_axes_group: ' // & 'Can not set x_cell_method for rank<2.') @@ -930,25 +929,46 @@ subroutine post_xy_average(diag_cs, diag, field) type(diag_type), intent(in) :: diag !< This diagnostic real, target, intent(in) :: field(:,:,:) !< Diagnostic field type(diag_ctrl), intent(in) :: diag_cs !< Diagnostics mediator control structure + ! Local variable - integer :: nk, i, j, k real, dimension(size(field,3)) :: averaged_field logical :: staggered_in_x, staggered_in_y, used + integer :: nz, remap_nz, coord - if (.not.diag%axes%is_h_point) then - call MOM_error(FATAL, 'post_xy_average: Horizontally averaged diagnostic not implemented yet.') + if (.not. diag_cs%ave_enabled) then + return endif staggered_in_x = diag%axes%is_u_point .or. diag%axes%is_q_point staggered_in_y = diag%axes%is_v_point .or. diag%axes%is_q_point - if (diag_cs%ave_enabled) then - call horizontally_average_diag_field(diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number), & - diag_cs%G, staggered_in_x, staggered_in_y, & + + if (diag%axes%is_native) then + call horizontally_average_diag_field(diag_cs%G, diag_cs%h, & + staggered_in_x, staggered_in_y, & + diag%axes%is_layer, diag%v_extensive, & + diag_cs%missing_value, field, averaged_field) + else + nz = size(field, 3) + coord = diag%axes%vertical_coordinate_number + remap_nz = diag_cs%diag_remap_cs(coord)%nz + + call assert(diag_cs%diag_remap_cs(coord)%initialized, & + 'post_xy_average: remap_cs not initialized.') + + call assert(IMPLIES(diag%axes%is_layer, nz == remap_nz), & + 'post_xy_average: layer field dimension mismatch.') + call assert(IMPLIES(.not. diag%axes%is_layer, nz == remap_nz+1), & + 'post_xy_average: interface field dimension mismatch.') + + call horizontally_average_diag_field(diag_cs%G, diag_cs%diag_remap_cs(coord)%h, & + staggered_in_x, staggered_in_y, & diag%axes%is_layer, diag%v_extensive, & diag_cs%missing_value, field, averaged_field) - used = send_data(diag%fms_xyave_diag_id, averaged_field, diag_cs%time_end, weight=diag_cs%time_int) endif + used = send_data(diag%fms_xyave_diag_id, averaged_field, diag_cs%time_end, & + weight=diag_cs%time_int) + end subroutine post_xy_average subroutine enable_averaging(time_int_in, time_end_in, diag_cs) diff --git a/src/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index 120812e19a..5acbed9c02 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -1,7 +1,7 @@ !> This module is used for runtime remapping of diagnostics to z star, sigma and !! rho vertical coordinates. It defines the diag_remap_ctrl type which !! represents a remapping of diagnostics to a particular vertical coordinate. -!! The module is It is used by the diag mediator module in the following way: +!! The module is used by the diag mediator module in the following way: !! 1) _init() is called to initialise a diag_remap_ctrl instance. !! 2) _configure_axes() is called to read the configuration file and set up the !! vertical coordinate / axes definitions. @@ -511,11 +511,11 @@ subroutine vertically_interpolate_diag_field(remap_cs, G, h, staggered_in_x, sta end subroutine vertically_interpolate_diag_field !> Horizontally average field -subroutine horizontally_average_diag_field(remap_cs, G, staggered_in_x, staggered_in_y, & +subroutine horizontally_average_diag_field(G, h, staggered_in_x, staggered_in_y, & is_layer, is_extensive, & missing_value, field, averaged_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 logical, intent(in) :: staggered_in_x !< True if the x-axis location is at u or q points logical, intent(in) :: staggered_in_y !< True if the y-axis location is at v or q points logical, intent(in) :: is_layer !< True if the z-axis location is at h points @@ -524,27 +524,16 @@ subroutine horizontally_average_diag_field(remap_cs, G, staggered_in_x, staggere real, dimension(:,:,:), intent(in) :: field !< The diagnostic field to be remapped real, dimension(:), intent(inout) :: averaged_field !< Field argument horizontally averaged ! Local variables - real, dimension(remap_cs%nz+1) :: vol_sum, stuff_sum ! nz+1 is needed for interface averages + real, dimension(size(field, 3)) :: vol_sum, stuff_sum ! nz+1 is needed for interface averages real :: v1, v2, total_volume, total_stuff, val - integer :: i, j, k + integer :: i, j, k, nz - call assert(remap_cs%initialized, 'horizontally_average_diag_field: remap_cs not initialized.') - if (is_layer) then - call assert(size(field, 3) == remap_cs%nz, & - 'horizontally_average_diag_field: Field vertical dimension does not match diag remap type.') - call assert(size(averaged_field, 1) == remap_cs%nz, & - 'horizontally_average_diag_field: Averaged field dimension does not match diag remap type.') - else - call assert(size(field, 3) == remap_cs%nz+1, & - 'horizontally_average_diag_field: Field vertical dimension does not match diag remap type.') - call assert(size(averaged_field, 1) == remap_cs%nz+1, & - 'horizontally_average_diag_field: Averaged field dimension does not match diag remap type.') - endif + nz = size(field, 3) if (staggered_in_x .and. .not. staggered_in_y) then if (is_layer) then ! U-points - do k=1,remap_cs%nz + do k=1,nz vol_sum(k) = 0. stuff_sum(k) = 0. if (is_extensive) then @@ -556,15 +545,15 @@ subroutine horizontally_average_diag_field(remap_cs, G, staggered_in_x, staggere enddo ; enddo else ! Intensive do j=G%jsc, G%jec ; do i=G%isc, G%iec - v1 = G%areaCu(I,j) * 0.5 * ( remap_cs%h(i,j,k) + remap_cs%h(i+1,j,k) ) - v2 = G%areaCu(I-1,j) * 0.5 * ( remap_cs%h(i,j,k) + remap_cs%h(i-1,j,k) ) + v1 = G%areaCu(I,j) * 0.5 * ( h(i,j,k) + h(i+1,j,k) ) + v2 = G%areaCu(I-1,j) * 0.5 * ( h(i,j,k) + h(i-1,j,k) ) vol_sum(k) = vol_sum(k) + 0.5 * ( v1 + v2 ) * G%mask2dT(i,j) stuff_sum(k) = stuff_sum(k) + 0.5 * ( v1 * field(I,j,k) + v2 * field(I-1,j,k) ) * G%mask2dT(i,j) enddo ; enddo endif enddo else ! Interface - do k=1,remap_cs%nz+1 + do k=1,nz vol_sum(k) = 0. stuff_sum(k) = 0. do j=G%jsc, G%jec ; do i=G%isc, G%iec @@ -578,7 +567,7 @@ subroutine horizontally_average_diag_field(remap_cs, G, staggered_in_x, staggere elseif (staggered_in_y .and. .not. staggered_in_x) then if (is_layer) then ! V-points - do k=1,remap_cs%nz + do k=1,nz vol_sum(k) = 0. stuff_sum(k) = 0. if (is_extensive) then @@ -590,15 +579,15 @@ subroutine horizontally_average_diag_field(remap_cs, G, staggered_in_x, staggere enddo ; enddo else ! Intensive do j=G%jsc, G%jec ; do i=G%isc, G%iec - v1 = G%areaCv(i,J) * 0.5 * ( remap_cs%h(i,j,k) + remap_cs%h(i,j+1,k) ) - v2 = G%areaCv(i,J-1) * 0.5 * ( remap_cs%h(i,j,k) + remap_cs%h(i,j-1,k) ) + v1 = G%areaCv(i,J) * 0.5 * ( h(i,j,k) + h(i,j+1,k) ) + v2 = G%areaCv(i,J-1) * 0.5 * ( h(i,j,k) + h(i,j-1,k) ) vol_sum(k) = vol_sum(k) + 0.5 * ( v1 + v2 ) * G%mask2dT(i,j) stuff_sum(k) = stuff_sum(k) + 0.5 * ( v1 * field(i,J,k) + v2 * field(i,J-1,k) ) * G%mask2dT(i,j) enddo ; enddo endif enddo else ! Interface - do k=1,remap_cs%nz+1 + do k=1,nz vol_sum(k) = 0. stuff_sum(k) = 0. do j=G%jsc, G%jec ; do i=G%isc, G%iec @@ -612,12 +601,12 @@ subroutine horizontally_average_diag_field(remap_cs, G, staggered_in_x, staggere elseif ((.not. staggered_in_x) .and. (.not. staggered_in_y)) then if (is_layer) then ! H-points - do k=1,remap_cs%nz + do k=1,nz vol_sum(k) = 0. stuff_sum(k) = 0. if (is_extensive) then do j=G%jsc, G%jec ; do i=G%isc, G%iec - if (G%mask2dT(i,j)>0. .and. remap_cs%h(i,j,k)>0.) then + if (G%mask2dT(i,j)>0. .and. h(i,j,k)>0.) then v1 = G%areaT(i,j) vol_sum(k) = vol_sum(k) + v1 stuff_sum(k) = stuff_sum(k) + v1 * field(i,j,k) @@ -625,8 +614,8 @@ subroutine horizontally_average_diag_field(remap_cs, G, staggered_in_x, staggere enddo ; enddo else ! Intensive do j=G%jsc, G%jec ; do i=G%isc, G%iec - if (G%mask2dT(i,j)>0. .and. remap_cs%h(i,j,k)>0.) then - v1 = G%areaT(i,j) * remap_cs%h(i,j,k) + if (G%mask2dT(i,j)>0. .and. h(i,j,k)>0.) then + v1 = G%areaT(i,j) * h(i,j,k) vol_sum(k) = vol_sum(k) + v1 stuff_sum(k) = stuff_sum(k) + v1 * field(i,j,k) endif @@ -634,7 +623,7 @@ subroutine horizontally_average_diag_field(remap_cs, G, staggered_in_x, staggere endif enddo else ! Interface - do k=1,remap_cs%nz+1 + do k=1,nz vol_sum(k) = 0. stuff_sum(k) = 0. do j=G%jsc, G%jec ; do i=G%isc, G%iec @@ -651,10 +640,10 @@ subroutine horizontally_average_diag_field(remap_cs, G, staggered_in_x, staggere call assert(.false., 'horizontally_average_diag_field: Q point averaging is not coded yet.') endif - call sum_across_PEs(vol_sum, remap_cs%nz) - call sum_across_PEs(stuff_sum, remap_cs%nz) + call sum_across_PEs(vol_sum, nz) + call sum_across_PEs(stuff_sum, nz) - do k=1,remap_cs%nz + do k=1,nz if (vol_sum(k)>0.) then averaged_field(k) = stuff_sum(k) / vol_sum(k) else From d0f0e7a70b759eecd7a842b0b34262f08f9897e5 Mon Sep 17 00:00:00 2001 From: Nicholas Hannah Date: Fri, 16 Dec 2016 15:30:25 +1100 Subject: [PATCH 2/2] Fix diagnostic (v)umo_2d post. These diags were being registered but not posted properly. #334 --- src/core/MOM.F90 | 4 ++-- src/framework/MOM_diag_mediator.F90 | 3 --- 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index fa5311b9e5..36add03068 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1215,7 +1215,7 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) umo2d(:,:) = umo2d(:,:) + CS%uhtr(:,:,k) enddo umo2d(:,:) = umo2d(:,:) * ( GV%H_to_kg_m2 / CS%dt_trans ) - call post_data(CS%id_umo, umo2d, CS%diag) + call post_data(CS%id_umo_2d, umo2d, CS%diag) endif if (CS%id_umo > 0) then ! Convert to kg/s. Modifying the array for diagnostics is allowed here since it is set to zero immediately below @@ -1228,7 +1228,7 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) vmo2d(:,:) = vmo2d(:,:) + CS%vhtr(:,:,k) enddo vmo2d(:,:) = vmo2d(:,:) * ( GV%H_to_kg_m2 / CS%dt_trans ) - call post_data(CS%id_vmo, vmo2d, CS%diag) + call post_data(CS%id_vmo_2d, vmo2d, CS%diag) endif if (CS%id_vmo > 0) then ! Convert to kg/s. Modifying the array for diagnostics is allowed here since it is set to zero immediately below diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 75fd384599..db3052c756 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -461,7 +461,6 @@ subroutine define_axes_group(diag_cs, handles, axes, nz, vertical_coordinate_num axes%rank = n axes%handles(:) = handles(:) axes%diag_cs => diag_cs ! A [circular] link back to the diag_cs structure - if (present(x_cell_method)) then if (axes%rank<2) call MOM_error(FATAL, 'define_axes_group: ' // & 'Can not set x_cell_method for rank<2.') @@ -929,7 +928,6 @@ subroutine post_xy_average(diag_cs, diag, field) type(diag_type), intent(in) :: diag !< This diagnostic real, target, intent(in) :: field(:,:,:) !< Diagnostic field type(diag_ctrl), intent(in) :: diag_cs !< Diagnostics mediator control structure - ! Local variable real, dimension(size(field,3)) :: averaged_field logical :: staggered_in_x, staggered_in_y, used @@ -968,7 +966,6 @@ subroutine post_xy_average(diag_cs, diag, field) used = send_data(diag%fms_xyave_diag_id, averaged_field, diag_cs%time_end, & weight=diag_cs%time_int) - end subroutine post_xy_average subroutine enable_averaging(time_int_in, time_end_in, diag_cs)