From af2b22185108d818f4adcf28989967a9cb56148d Mon Sep 17 00:00:00 2001 From: Nicholas Hannah Date: Tue, 9 Jun 2015 15:48:17 -0700 Subject: [PATCH 01/18] Support remapping of variables on u, v grids. #62 --- src/framework/MOM_diag_mediator.F90 | 157 ++++++++++++++++++++-------- 1 file changed, 115 insertions(+), 42 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index a6022ab3de..4e9c96cc9c 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -114,7 +114,7 @@ module MOM_diag_mediator type(axesType) :: axesBi, axesTi, axesCui, axesCvi type(axesType) :: axesB1, axesT1, axesCu1, axesCv1 type(axesType) :: axesZi, axesZL - type(axesType) :: axesTzi, axesTZL + type(axesType) :: axesTzi, axesTZL, axesBZL, axesCuZL, axesCvZL ! Mask arrays for diagnostics real, dimension(:,:), pointer :: mask2dT => null() @@ -297,8 +297,11 @@ subroutine set_axes_info(G, param_file, diag_cs, set_vertical) id_zzi = 1 endif - ! Axis for z remapping + ! Axes for z remapping call defineAxes(diag_cs, (/ id_xh, id_yh, id_zzl /), diag_cs%axesTZL) + call defineAxes(diag_cs, (/ id_xq, id_yq, id_zzL /), diag_cs%axesBZL) + call defineAxes(diag_cs, (/ id_xq, id_yh, id_zzL /), diag_cs%axesCuZL) + call defineAxes(diag_cs, (/ id_xh, id_yq, id_zzL /), diag_cs%axesCvZL) ! Vertical axes for the interfaces and layers call defineAxes(diag_cs, (/ id_zi /), diag_cs%axesZi) @@ -639,7 +642,8 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) type(remapping_CS) :: remap_cs type(regridding_CS) :: regrid_cs integer :: nz_src, i, j, k - real, dimension(diag_cs%nz_remap) :: h_new, h_remap + real, dimension(diag_cs%nz_remap) :: h_dest, h_remap + real, dimension(size(diag_cs%h, 3)) :: h_src real, dimension(diag_cs%nz_remap+1) :: dz, zi_tmp call assert(size(field, 3) == size(diag_cs%h, 3), & @@ -662,16 +666,26 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) if (diag%mask3d(i,j, 1) == 0.) cycle endif + ! Calculate source h depending on grid + if (is_u_axes(diag%remap_axes, diag_cs)) then + h_src(:) = 0.5 * (diag_cs%h(i,j,:) + diag_cs%h(i+1,j,:)) + elseif (is_v_axes(diag%remap_axes, diag_cs)) then + h_src(:) = 0.5 * (diag_cs%h(i,j,:) + diag_cs%h(i,j+1,:)) + else + h_src(:) = diag_cs%h(i, j, :) + endif + print*, 'h: ', h_src + ! Calculate thicknesses for new Z* using nominal thicknesses from ! h_remap, current bathymetry and total thickness. call buildGridZstarColumn(regrid_cs, diag_cs%nz_remap, diag_cs%G%bathyT(i, j), & - sum(diag_cs%h(i, j, :)), zi_tmp) + sum(h_src(:)), zi_tmp) ! Calculate how much thicknesses change between source and dest grids, do ! remapping zi_tmp = -zi_tmp - h_new(:) = zi_tmp(2:) - zi_tmp(:size(zi_tmp)-1) - call dzFromH1H2(nz_src, diag_cs%h(i, j, :), diag_cs%nz_remap, h_new, dz) - call remapping_core(remap_cs, nz_src, diag_cs%h(i, j, :), field(i,j,:), & + h_dest(:) = zi_tmp(2:) - zi_tmp(:size(zi_tmp)-1) + call dzFromH1H2(nz_src, h_src(:), diag_cs%nz_remap, h_dest, dz) + call remapping_core(remap_cs, nz_src, h_src(:), field(i,j,:), & diag_cs%nz_remap, dz, remapped_field(i, j, :)) ! Lower levels of the remapped data get squashed to follow bathymetry, @@ -907,7 +921,7 @@ function register_diag_field(module_name, field_name, axes, init_time, & call alloc_diag_with_id(primary_id, diag_cs, diag) call assert(associated(diag), 'register_diag_field: diag allocation failed') diag%fms_diag_id = fms_id - call set_diag_mask(diag, diag_cs, axes) + call set_diag_mask_and_remap_axes(diag, diag_cs, axes) endif ! Set up the CMOR variation of the diagnostic @@ -942,34 +956,30 @@ function register_diag_field(module_name, field_name, axes, init_time, & ! In the case where there is no primary, it will become the primary. call alloc_diag_with_id(primary_id, diag_cs, cmor_diag) cmor_diag%fms_diag_id = fms_id - call set_diag_mask(cmor_diag, diag_cs, axes) + call set_diag_mask_and_remap_axes(cmor_diag, diag_cs, axes) endif endif - ! Remap to z vertical coordinate, note that only - ! diagnostics on T points and layers (not interfaces) are supported, - ! for other diagnostics the old remapping approach can still be used - if (axes%id == diag_cs%axesTL%id) then - fms_id = register_diag_field_fms(module_name//'_z_new', field_name, diag_cs%axesTZL%handles, & - init_time, long_name=long_name, units=units, missing_value=MOM_missing_value, & - range=range, mask_variant=mask_variant, standard_name=standard_name, & - verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, & - interp_method=interp_method, tile_count=tile_count) - if (fms_id /= DIAG_FIELD_NOT_FOUND) then - ! Check that we have read in necessary remapping information from file - if (.not. allocated(diag_cs%zi_remap)) then - call MOM_error(FATAL, 'register_diag_field: Request to regrid but no '// & - 'destination grid spec provided, see param DIAG_REMAP_Z_GRID_DEF') - endif - - if (primary_id == -1) then - primary_id = get_new_diag_id(diag_cs) - endif - call alloc_diag_with_id(primary_id, diag_cs, z_remap_diag) - z_remap_diag%fms_diag_id = fms_id - z_remap_diag%remap_axes => diag_cs%axesTZL - call set_diag_mask(z_remap_diag, diag_cs, diag_cs%axesTZL) + ! Remap to z vertical coordinate, note that only diagnostics on layers + ! (not interfaces) are supported, + if (get_diag_field_id_fms(module_name//'_z_new', field_name) /= DIAG_FIELD_NOT_FOUND & + .and. is_layer_axes(axes, diag_cs)) then + if (.not. allocated(diag_cs%zi_remap)) then + call MOM_error(FATAL, 'register_diag_field: Request to regrid but no '// & + 'destination grid spec provided, see param DIAG_REMAP_Z_GRID_DEF') + endif + if (primary_id == -1) then + primary_id = get_new_diag_id(diag_cs) endif + call alloc_diag_with_id(primary_id, diag_cs, z_remap_diag) + call set_diag_mask_and_remap_axes(z_remap_diag, diag_cs, axes) + fms_id = register_diag_field_fms(module_name//'_z_new', field_name, & + z_remap_diag%remap_axes%handles, & + init_time, long_name=long_name, units=units, missing_value=MOM_missing_value, & + range=range, mask_variant=mask_variant, standard_name=standard_name, & + verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, & + interp_method=interp_method, tile_count=tile_count) + z_remap_diag%fms_diag_id = fms_id endif ! Document diagnostics in list of available diagnostics @@ -981,7 +991,7 @@ function register_diag_field(module_name, field_name, axes, init_time, & posted_cmor_long_name, posted_cmor_units, & posted_cmor_standard_name) endif - if (axes%id == diag_cs%axesTL%id) then + if (is_layer_axes(axes, diag_cs)) then call log_available_diag(associated(z_remap_diag), module_name//'_z_new', field_name, & long_name, units, standard_name) endif @@ -1460,25 +1470,30 @@ function i2s(a,n_in) i2s = adjustl(i2s) end function i2s -subroutine set_diag_mask(diag, diag_cs, axes) +subroutine set_diag_mask_and_remap_axes(diag, diag_cs, axes) ! Associate a mask with the a diagnostic based on the axes info. - type(diag_ctrl), intent(inout) :: diag_cs + type(diag_ctrl), target, intent(inout) :: diag_cs type(diag_type), pointer, intent(out) :: diag type(axesType), intent(in) :: axes + diag%mask2d => null() + diag%mask3d => null() + diag%remap_axes => null() + if (axes%rank .eq. 3) then - diag%mask2d => null() - diag%mask3d => null() - if ((axes%id .eq. diag_cs%axesTL%id) .or. & - (axes%id .eq. diag_cs%axesTZL%id)) then + if ((axes%id .eq. diag_cs%axesTL%id)) then diag%mask3d => diag_cs%mask3dTL + diag%remap_axes => diag_cs%axesTZL elseif(axes%id .eq. diag_cs%axesBL%id) then diag%mask3d => diag_cs%mask3dBuL + diag%remap_axes => diag_cs%axesBZL elseif(axes%id .eq. diag_cs%axesCuL%id ) then diag%mask3d => diag_cs%mask3dCuL + diag%remap_axes => diag_cs%axesCuZL elseif(axes%id .eq. diag_cs%axesCvL%id) then diag%mask3d => diag_cs%mask3dCvL + diag%remap_axes => diag_cs%axesCvZL elseif(axes%id .eq. diag_cs%axesTi%id) then diag%mask3d => diag_cs%mask3dTi elseif(axes%id .eq. diag_cs%axesBi%id) then @@ -1491,8 +1506,6 @@ subroutine set_diag_mask(diag, diag_cs, axes) !call assert(associated(diag%mask3d), "MOM_diag_mediator.F90: Invalid 3d axes id") elseif(axes%rank .eq. 2) then - diag%mask2d => null() - diag%mask3d => null() if (axes%id .eq. diag_cs%axesT1%id) then diag%mask2d => diag_cs%mask2dT elseif(axes%id .eq. diag_cs%axesB1%id) then @@ -1506,7 +1519,67 @@ subroutine set_diag_mask(diag, diag_cs, axes) !call assert(associated(diag%mask2d), "MOM_diag_mediator.F90: Invalid 2d axes id") endif -end subroutine set_diag_mask +end subroutine set_diag_mask_and_remap_axes + +function is_layer_axes(axes, diag_cs) + + type(axesType), intent(in) :: axes + type(diag_ctrl), intent(in) :: diag_cs + + logical :: is_layer_axes + + is_layer_axes = .false. + + if (axes%id == diag_cs%axesTZL%id .or. & + axes%id == diag_cs%axesBZL%id .or. & + axes%id == diag_cs%axesCuZL%id .or. & + axes%id == diag_cs%axesCvZL%id .or. & + axes%id == diag_cs%axesZi%id .or. & + axes%id == diag_cs%axesZL%id .or. & + axes%id == diag_cs%axesTL%id .or. & + axes%id == diag_cs%axesBL%id .or. & + axes%id == diag_cs%axesCuL%id .or. & + axes%id == diag_cs%axesCvL%id) then + is_layer_axes = .true. + endif + +end function is_layer_axes + +function is_u_axes(axes, diag_cs) + + type(axesType), intent(in) :: axes + type(diag_ctrl), intent(in) :: diag_cs + + logical :: is_u_axes + + is_u_axes = .false. + + if (axes%id == diag_cs%axesCuZL%id .or. & + axes%id == diag_cs%axesCuL%id .or. & + axes%id == diag_cs%axesCui%id .or. & + axes%id == diag_cs%axesCu1%id) then + is_u_axes = .true. + endif + +end function is_u_axes + +function is_v_axes(axes, diag_cs) + + type(axesType), intent(in) :: axes + type(diag_ctrl), intent(in) :: diag_cs + + logical :: is_v_axes + + is_v_axes = .false. + + if (axes%id == diag_cs%axesCuZL%id .or. & + axes%id == diag_cs%axesCuL%id .or. & + axes%id == diag_cs%axesCui%id .or. & + axes%id == diag_cs%axesCu1%id) then + is_v_axes = .true. + endif + +end function is_v_axes ! Allocate a new diagnostic id, it may be necessary to expand the diagnostics ! array. From a1510d14e392bf1e50ce950a03270feda9b7154d Mon Sep 17 00:00:00 2001 From: Nicholas Hannah Date: Wed, 10 Jun 2015 08:35:44 -0700 Subject: [PATCH 02/18] Move diag remapping into per-column subroutine, easier to remap u, v grids. #62 --- src/framework/MOM_diag_mediator.F90 | 89 +++++++++++++++++++---------- 1 file changed, 59 insertions(+), 30 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 4e9c96cc9c..a63a7a2074 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -641,10 +641,9 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) type(remapping_CS) :: remap_cs type(regridding_CS) :: regrid_cs - integer :: nz_src, i, j, k - real, dimension(diag_cs%nz_remap) :: h_dest, h_remap - real, dimension(size(diag_cs%h, 3)) :: h_src - real, dimension(diag_cs%nz_remap+1) :: dz, zi_tmp + integer :: nz_src + integer :: i, j, k, ilow, ihigh, jlow, jhigh + real, dimension(diag_cs%nz_remap) :: h_dest call assert(size(field, 3) == size(diag_cs%h, 3), & 'Remap field and thickness z-axes do not match.') @@ -653,40 +652,44 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) remapped_field = diag_cs%missing_value nz_src = size(field, 3) - ! Nominal thicknesses to remap to - h_remap(:) = diag_cs%zi_remap(2:) - diag_cs%zi_remap(:diag_cs%nz_remap) + ! In itialise remapping and set nominal thicknesses to remap to call initialize_regridding(diag_cs%nz_remap, 'Z*', 'PPM_IH4', regrid_cs) call setRegriddingMinimumThickness(diag_cs%G%Angstrom, regrid_cs) - call setCoordinateResolution(h_remap, regrid_cs) + h_dest(:) = diag_cs%zi_remap(2:) - diag_cs%zi_remap(:diag_cs%nz_remap) + call setCoordinateResolution(h_dest, regrid_cs) call initialize_remapping(nz_src, 'PPM_IH4', remap_cs) - do j=RANGE_J(field) - do i=RANGE_I(field) + print*, 'field: ', size(field), lbound(field, 1), ubound(field, 1), lbound(field, 2), ubound(field, 2) + print*, 'h: ', size(field), lbound(field, 1), ubound(field, 1), lbound(field, 2), ubound(field, 2) + print*, 'G%jscB, G%jecB: ', diag_cs%G%jscB, diag_cs%G%jecB + print*, 'G%isc, G%iec: ', diag_cs%G%isc, diag_cs%G%iec + + ilow = diag_cs%G%isc + ihigh = diag_cs%G%iec + jlow = diag_cs%G%jsc + jhigh = diag_cs%G%jec + + if (is_u_axes(diag%remap_axes, diag_cs)) then + print*, 'u_axes' + ilow = diag_cs%G%iscB + ihigh = diag_cs%G%iecB + elseif (is_v_axes(diag%remap_axes, diag_cs)) then + print*, 'v_axes' + jlow = diag_cs%G%jscB + jhigh = diag_cs%G%jecB + else + print*, 'T_axes' + endif + + do j=jlow, jhigh + do i=ilow, ihigh if (associated(diag%mask3d)) then if (diag%mask3d(i,j, 1) == 0.) cycle endif - ! Calculate source h depending on grid - if (is_u_axes(diag%remap_axes, diag_cs)) then - h_src(:) = 0.5 * (diag_cs%h(i,j,:) + diag_cs%h(i+1,j,:)) - elseif (is_v_axes(diag%remap_axes, diag_cs)) then - h_src(:) = 0.5 * (diag_cs%h(i,j,:) + diag_cs%h(i,j+1,:)) - else - h_src(:) = diag_cs%h(i, j, :) - endif - print*, 'h: ', h_src - - ! Calculate thicknesses for new Z* using nominal thicknesses from - ! h_remap, current bathymetry and total thickness. - call buildGridZstarColumn(regrid_cs, diag_cs%nz_remap, diag_cs%G%bathyT(i, j), & - sum(h_src(:)), zi_tmp) - ! Calculate how much thicknesses change between source and dest grids, do - ! remapping - zi_tmp = -zi_tmp - h_dest(:) = zi_tmp(2:) - zi_tmp(:size(zi_tmp)-1) - call dzFromH1H2(nz_src, h_src(:), diag_cs%nz_remap, h_dest, dz) - call remapping_core(remap_cs, nz_src, h_src(:), field(i,j,:), & - diag_cs%nz_remap, dz, remapped_field(i, j, :)) + call remap_column_to_z(regrid_cs, remap_cs, nz_src, diag_cs%nz_remap, & + diag_cs%h(i, j, :), diag_cs%G%bathyT(i, j), & + field(i, j, :), 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 @@ -702,6 +705,32 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) end subroutine remap_diag_to_z +subroutine remap_column_to_z(regrid_cs, remap_cs, nz_src, nz_dest, h, depth, & + field, remapped_field) + + type(remapping_CS), intent(in) :: remap_cs + type(regridding_CS), intent(in) :: regrid_cs + integer, intent(in) :: nz_src, nz_dest + real, dimension(:), intent(in) :: h, field + real, intent(in) :: depth + real, dimension(:), intent(inout) :: remapped_field + + real, dimension(nz_dest+1) :: dz, zi_tmp + real, dimension(nz_dest) :: h_dest + + ! Calculate thicknesses for new Z* using nominal thicknesses from + ! h_remap, current bathymetry and total thickness. + call buildGridZstarColumn(regrid_cs, nz_dest, depth, sum(h(:)), zi_tmp) + ! Calculate how much thicknesses change between source and dest grids, do + ! remapping + zi_tmp = -zi_tmp + h_dest(:) = zi_tmp(2:) - zi_tmp(:size(zi_tmp)-1) + call dzFromH1H2(nz_src, h(:), nz_dest, h_dest, dz) + call remapping_core(remap_cs, nz_src, h(:), field(:), nz_dest, dz, & + remapped_field(:)) + +end subroutine remap_column_to_z + subroutine post_data_3d_low(diag, field, diag_cs, is_static, mask) type(diag_type), intent(in) :: diag real, intent(in) :: field(:,:,:) From d9297bde2ed6d7fc95c8c1273d18b3a3c2ff79b4 Mon Sep 17 00:00:00 2001 From: Nicholas Hannah Date: Wed, 10 Jun 2015 09:06:49 -0700 Subject: [PATCH 03/18] Set 3d mask and axes separately. #62 --- src/framework/MOM_diag_mediator.F90 | 38 +++++++++++++++++++++-------- 1 file changed, 28 insertions(+), 10 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index a63a7a2074..15c5ecf8bd 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -950,7 +950,7 @@ function register_diag_field(module_name, field_name, axes, init_time, & call alloc_diag_with_id(primary_id, diag_cs, diag) call assert(associated(diag), 'register_diag_field: diag allocation failed') diag%fms_diag_id = fms_id - call set_diag_mask_and_remap_axes(diag, diag_cs, axes) + call set_diag_mask(diag, diag_cs, axes) endif ! Set up the CMOR variation of the diagnostic @@ -985,7 +985,7 @@ function register_diag_field(module_name, field_name, axes, init_time, & ! In the case where there is no primary, it will become the primary. call alloc_diag_with_id(primary_id, diag_cs, cmor_diag) cmor_diag%fms_diag_id = fms_id - call set_diag_mask_and_remap_axes(cmor_diag, diag_cs, axes) + call set_diag_mask(cmor_diag, diag_cs, axes) endif endif @@ -1001,7 +1001,8 @@ function register_diag_field(module_name, field_name, axes, init_time, & primary_id = get_new_diag_id(diag_cs) endif call alloc_diag_with_id(primary_id, diag_cs, z_remap_diag) - call set_diag_mask_and_remap_axes(z_remap_diag, diag_cs, axes) + call set_diag_mask(z_remap_diag, diag_cs, axes) + call set_diag_remap_axes(z_remap_diag, diag_cs, axes) fms_id = register_diag_field_fms(module_name//'_z_new', field_name, & z_remap_diag%remap_axes%handles, & init_time, long_name=long_name, units=units, missing_value=MOM_missing_value, & @@ -1499,7 +1500,29 @@ function i2s(a,n_in) i2s = adjustl(i2s) end function i2s -subroutine set_diag_mask_and_remap_axes(diag, diag_cs, axes) +subroutine set_diag_remap_axes(diag, diag_cs, axes) + ! Associate a remapping axes with the a diagnostic based on the axes info. + type(diag_ctrl), target, intent(inout) :: diag_cs + type(diag_type), pointer, intent(out) :: diag + type(axesType), intent(in) :: axes + + diag%remap_axes => null() + + if (axes%rank .eq. 3) then + if ((axes%id .eq. diag_cs%axesTL%id)) then + diag%remap_axes => diag_cs%axesTZL + elseif(axes%id .eq. diag_cs%axesBL%id) then + diag%remap_axes => diag_cs%axesBZL + elseif(axes%id .eq. diag_cs%axesCuL%id ) then + diag%remap_axes => diag_cs%axesCuZL + elseif(axes%id .eq. diag_cs%axesCvL%id) then + diag%remap_axes => diag_cs%axesCvZL + endif + endif + +end subroutine set_diag_remap_axes + +subroutine set_diag_mask(diag, diag_cs, axes) ! Associate a mask with the a diagnostic based on the axes info. type(diag_ctrl), target, intent(inout) :: diag_cs @@ -1508,21 +1531,16 @@ subroutine set_diag_mask_and_remap_axes(diag, diag_cs, axes) diag%mask2d => null() diag%mask3d => null() - diag%remap_axes => null() if (axes%rank .eq. 3) then if ((axes%id .eq. diag_cs%axesTL%id)) then diag%mask3d => diag_cs%mask3dTL - diag%remap_axes => diag_cs%axesTZL elseif(axes%id .eq. diag_cs%axesBL%id) then diag%mask3d => diag_cs%mask3dBuL - diag%remap_axes => diag_cs%axesBZL elseif(axes%id .eq. diag_cs%axesCuL%id ) then diag%mask3d => diag_cs%mask3dCuL - diag%remap_axes => diag_cs%axesCuZL elseif(axes%id .eq. diag_cs%axesCvL%id) then diag%mask3d => diag_cs%mask3dCvL - diag%remap_axes => diag_cs%axesCvZL elseif(axes%id .eq. diag_cs%axesTi%id) then diag%mask3d => diag_cs%mask3dTi elseif(axes%id .eq. diag_cs%axesBi%id) then @@ -1548,7 +1566,7 @@ subroutine set_diag_mask_and_remap_axes(diag, diag_cs, axes) !call assert(associated(diag%mask2d), "MOM_diag_mediator.F90: Invalid 2d axes id") endif -end subroutine set_diag_mask_and_remap_axes +end subroutine set_diag_mask function is_layer_axes(axes, diag_cs) From d39cf62ce58708a5c8942bed8b855926fb7b8b2d Mon Sep 17 00:00:00 2001 From: Nicholas Hannah Date: Wed, 10 Jun 2015 13:24:25 -0700 Subject: [PATCH 04/18] Check relative and absolute errors when doing conservative remapping. #62 --- src/ALE/MOM_remapping.F90 | 23 +++++++++++++++-------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/src/ALE/MOM_remapping.F90 b/src/ALE/MOM_remapping.F90 index 487d12ce8c..0797cd928e 100644 --- a/src/ALE/MOM_remapping.F90 +++ b/src/ALE/MOM_remapping.F90 @@ -534,15 +534,22 @@ subroutine remapping_core( CS, n0, h0, u0, n1, dx, u1 ) if (k<=n0) then; hTmp = h0(k); else; hTmp = 0.; endif z0 = z0 + hTmp ; z1 = z1 + ( hTmp + ( dx(k+1) - dx(k) ) ) enddo + ! Maximum error based on guess at maximum roundoff if (abs(totalHU2-totalHU0) > (err0+err2)*max(real(n0), real(n1)) .and. (err0+err2)/=0.) then - write(0,*) 'h0=',h0 - write(0,*) 'hf=',h0(1:n1)+dx(2:n1+1)-dx(1:n1) - write(0,*) 'u0=',u0 - write(0,*) 'u1=',u1 - write(0,*) 'total HU0,HUf,f-0=',totalHU0,totalHU2,totalHU2-totalHU0 - write(0,*) 'err0,errF=',err0,err2 - call MOM_error( FATAL, 'MOM_remapping, remapping_core: '//& - 'Total stuff on h0 and hF differ by more than roundoff' ) + ! Maximum relative error + if (abs(totalHU2-totalHU0) / totalHU2 > 1e-6) then + ! Maximum absolute error + if (abs(totalHU2-totalHU0) > 1e-18) then + write(0,*) 'h0=',h0 + write(0,*) 'hf=',h0(1:n1)+dx(2:n1+1)-dx(1:n1) + write(0,*) 'u0=',u0 + write(0,*) 'u1=',u1 + write(0,*) 'total HU0,HUf,f-0=',totalHU0,totalHU2,totalHU2-totalHU0 + write(0,*) 'err0,errF=',err0,err2 + call MOM_error( FATAL, 'MOM_remapping, remapping_core: '//& + 'Total stuff on h0 and hF differ by more than roundoff' ) + endif + endif endif #endif From 759adfdfe1f5f1a9d614b020b7ec1bc59f48aa69 Mon Sep 17 00:00:00 2001 From: Nicholas Hannah Date: Wed, 10 Jun 2015 13:25:49 -0700 Subject: [PATCH 05/18] Estimate h and depth for u, v, B grids when doing diag remapping. #62 --- src/framework/MOM_diag_mediator.F90 | 134 ++++++++++++++++++---------- 1 file changed, 88 insertions(+), 46 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 15c5ecf8bd..9b1fcc0077 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -642,8 +642,10 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) type(remapping_CS) :: remap_cs type(regridding_CS) :: regrid_cs integer :: nz_src - integer :: i, j, k, ilow, ihigh, jlow, jhigh + integer :: i, j real, dimension(diag_cs%nz_remap) :: h_dest + real, dimension(size(field, 3)) :: h_src + real :: depth call assert(size(field, 3) == size(diag_cs%h, 3), & 'Remap field and thickness z-axes do not match.') @@ -659,64 +661,76 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) call setCoordinateResolution(h_dest, regrid_cs) call initialize_remapping(nz_src, 'PPM_IH4', remap_cs) - print*, 'field: ', size(field), lbound(field, 1), ubound(field, 1), lbound(field, 2), ubound(field, 2) - print*, 'h: ', size(field), lbound(field, 1), ubound(field, 1), lbound(field, 2), ubound(field, 2) - print*, 'G%jscB, G%jecB: ', diag_cs%G%jscB, diag_cs%G%jecB - print*, 'G%isc, G%iec: ', diag_cs%G%isc, diag_cs%G%iec - - ilow = diag_cs%G%isc - ihigh = diag_cs%G%iec - jlow = diag_cs%G%jsc - jhigh = diag_cs%G%jec - if (is_u_axes(diag%remap_axes, diag_cs)) then - print*, 'u_axes' - ilow = diag_cs%G%iscB - ihigh = diag_cs%G%iecB + 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 + h_src(:) = 0.5 * (diag_cs%h(i,j,:) + diag_cs%h(i+1,j,:)) + depth = 0.5 * (diag_cs%G%bathyT(i, j) + diag_cs%G%bathyT(i+1, j)) + call remap_column_to_z(regrid_cs, remap_cs, nz_src, diag_cs%nz_remap, & + h_src(:), depth, diag_cs%zi_remap(:), & + field(i, j, :), remapped_field(i, j, :), diag_cs%missing_value) + enddo + enddo + elseif (is_v_axes(diag%remap_axes, diag_cs)) then - print*, 'v_axes' - jlow = diag_cs%G%jscB - jhigh = diag_cs%G%jecB - else - print*, 'T_axes' - endif + 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 + h_src(:) = 0.5 * (diag_cs%h(i,j,:) + diag_cs%h(i,j+1,:)) + depth = 0.5 * (diag_cs%G%bathyT(i, j) + diag_cs%G%bathyT(i, j+1)) + call remap_column_to_z(regrid_cs, remap_cs, nz_src, diag_cs%nz_remap, & + h_src(:), depth, diag_cs%zi_remap(:), & + field(i, j, :), remapped_field(i, j, :), diag_cs%missing_value) + enddo + enddo - do j=jlow, jhigh - do i=ilow, ihigh - if (associated(diag%mask3d)) then - if (diag%mask3d(i,j, 1) == 0.) cycle - endif + 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_src(:) = 0.5 * (diag_cs%h(i,j,:) + diag_cs%h(i+1,j+1,:)) + depth = 0.5 * (diag_cs%G%bathyT(i, j) + diag_cs%G%bathyT(i+1, j+1)) + call remap_column_to_z(regrid_cs, remap_cs, nz_src, diag_cs%nz_remap, & + h_src(:), depth, diag_cs%zi_remap(:), & + field(i, j, :), remapped_field(i, j, :), diag_cs%missing_value) + enddo + enddo - call remap_column_to_z(regrid_cs, remap_cs, nz_src, diag_cs%nz_remap, & - diag_cs%h(i, j, :), diag_cs%G%bathyT(i, j), & - field(i, j, :), 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, diag_cs%nz_remap - if (diag_cs%zi_remap(k) >= diag_cs%G%bathyT(i, j)) then - remapped_field(i, j, k:diag_cs%nz_remap) = diag_cs%missing_value - exit + else + 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 + call remap_column_to_z(regrid_cs, remap_cs, nz_src, diag_cs%nz_remap, & + diag_cs%h(i, j, :), diag_cs%G%bathyT(i, j), diag_cs%zi_remap(:), & + field(i, j, :), remapped_field(i, j, :), diag_cs%missing_value) enddo enddo - enddo + endif end subroutine remap_diag_to_z subroutine remap_column_to_z(regrid_cs, remap_cs, nz_src, nz_dest, h, depth, & - field, remapped_field) + zi_remap, field, remapped_field, missing_value) type(remapping_CS), intent(in) :: remap_cs type(regridding_CS), intent(in) :: regrid_cs integer, intent(in) :: nz_src, nz_dest - real, dimension(:), intent(in) :: h, field - real, intent(in) :: depth + real, dimension(:), intent(in) :: h, field, zi_remap + real, intent(in) :: depth, missing_value real, dimension(:), intent(inout) :: remapped_field real, dimension(nz_dest+1) :: dz, zi_tmp real, dimension(nz_dest) :: h_dest + integer :: k ! Calculate thicknesses for new Z* using nominal thicknesses from ! h_remap, current bathymetry and total thickness. @@ -729,6 +743,16 @@ subroutine remap_column_to_z(regrid_cs, remap_cs, nz_src, nz_dest, h, depth, & call remapping_core(remap_cs, nz_src, h(:), field(:), nz_dest, dz, & remapped_field(:)) + ! 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_column_to_z subroutine post_data_3d_low(diag, field, diag_cs, is_static, mask) @@ -992,7 +1016,7 @@ function register_diag_field(module_name, field_name, axes, init_time, & ! Remap to z vertical coordinate, note that only diagnostics on layers ! (not interfaces) are supported, if (get_diag_field_id_fms(module_name//'_z_new', field_name) /= DIAG_FIELD_NOT_FOUND & - .and. is_layer_axes(axes, diag_cs)) then + .and. is_layer_axes(axes, diag_cs) .and. axes%rank == 3) then if (.not. allocated(diag_cs%zi_remap)) then call MOM_error(FATAL, 'register_diag_field: Request to regrid but no '// & 'destination grid spec provided, see param DIAG_REMAP_Z_GRID_DEF') @@ -1003,6 +1027,8 @@ function register_diag_field(module_name, field_name, axes, init_time, & call alloc_diag_with_id(primary_id, diag_cs, z_remap_diag) call set_diag_mask(z_remap_diag, diag_cs, axes) call set_diag_remap_axes(z_remap_diag, diag_cs, axes) + call assert(associated(z_remap_diag%remap_axes), & + 'register_diag_field: remap axes not set') fms_id = register_diag_field_fms(module_name//'_z_new', field_name, & z_remap_diag%remap_axes%handles, & init_time, long_name=long_name, units=units, missing_value=MOM_missing_value, & @@ -1021,7 +1047,7 @@ function register_diag_field(module_name, field_name, axes, init_time, & posted_cmor_long_name, posted_cmor_units, & posted_cmor_standard_name) endif - if (is_layer_axes(axes, diag_cs)) then + if (is_layer_axes(axes, diag_cs) .and. axes%rank == 3) then call log_available_diag(associated(z_remap_diag), module_name//'_z_new', field_name, & long_name, units, standard_name) endif @@ -1507,7 +1533,6 @@ subroutine set_diag_remap_axes(diag, diag_cs, axes) type(axesType), intent(in) :: axes diag%remap_axes => null() - if (axes%rank .eq. 3) then if ((axes%id .eq. diag_cs%axesTL%id)) then diag%remap_axes => diag_cs%axesTZL @@ -1578,13 +1603,12 @@ function is_layer_axes(axes, diag_cs) is_layer_axes = .false. if (axes%id == diag_cs%axesTZL%id .or. & - axes%id == diag_cs%axesBZL%id .or. & + !axes%id == diag_cs%axesBZL%id .or. & axes%id == diag_cs%axesCuZL%id .or. & axes%id == diag_cs%axesCvZL%id .or. & - axes%id == diag_cs%axesZi%id .or. & axes%id == diag_cs%axesZL%id .or. & axes%id == diag_cs%axesTL%id .or. & - axes%id == diag_cs%axesBL%id .or. & + !axes%id == diag_cs%axesBL%id .or. & axes%id == diag_cs%axesCuL%id .or. & axes%id == diag_cs%axesCvL%id) then is_layer_axes = .true. @@ -1628,6 +1652,24 @@ function is_v_axes(axes, diag_cs) end function is_v_axes +function is_B_axes(axes, diag_cs) + + type(axesType), intent(in) :: axes + type(diag_ctrl), intent(in) :: diag_cs + + logical :: is_B_axes + + is_B_axes = .false. + + if (axes%id == diag_cs%axesBZL%id .or. & + axes%id == diag_cs%axesBL%id .or. & + axes%id == diag_cs%axesBi%id .or. & + axes%id == diag_cs%axesB1%id) then + is_B_axes = .true. + endif + +end function is_B_axes + ! Allocate a new diagnostic id, it may be necessary to expand the diagnostics ! array. function get_new_diag_id(diag_cs) From 85497231a20225ff8e1e1947639739d4a0e480b6 Mon Sep 17 00:00:00 2001 From: Nicholas Hannah Date: Thu, 11 Jun 2015 08:00:17 -0700 Subject: [PATCH 06/18] Tune maximum round-off, relative and absolute errors when doing remapping. --- src/ALE/MOM_remapping.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ALE/MOM_remapping.F90 b/src/ALE/MOM_remapping.F90 index 0797cd928e..704fdad396 100644 --- a/src/ALE/MOM_remapping.F90 +++ b/src/ALE/MOM_remapping.F90 @@ -537,7 +537,7 @@ subroutine remapping_core( CS, n0, h0, u0, n1, dx, u1 ) ! Maximum error based on guess at maximum roundoff if (abs(totalHU2-totalHU0) > (err0+err2)*max(real(n0), real(n1)) .and. (err0+err2)/=0.) then ! Maximum relative error - if (abs(totalHU2-totalHU0) / totalHU2 > 1e-6) then + if (abs(totalHU2-totalHU0) / totalHU2 > 1e-09) then ! Maximum absolute error if (abs(totalHU2-totalHU0) > 1e-18) then write(0,*) 'h0=',h0 @@ -547,7 +547,7 @@ subroutine remapping_core( CS, n0, h0, u0, n1, dx, u1 ) write(0,*) 'total HU0,HUf,f-0=',totalHU0,totalHU2,totalHU2-totalHU0 write(0,*) 'err0,errF=',err0,err2 call MOM_error( FATAL, 'MOM_remapping, remapping_core: '//& - 'Total stuff on h0 and hF differ by more than roundoff' ) + 'Total stuff on h0 and hF differ by more than maximum errors' ) endif endif endif From dfc52bce37edf133fa57ed3e7a7fd80e760dd4bf Mon Sep 17 00:00:00 2001 From: Nicholas Hannah Date: Thu, 11 Jun 2015 08:03:46 -0700 Subject: [PATCH 07/18] Don't support remapping of diagnostics on B axes. #62 --- src/framework/MOM_diag_mediator.F90 | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 9b1fcc0077..cf24fb1ce3 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -674,7 +674,6 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) field(i, j, :), remapped_field(i, j, :), diag_cs%missing_value) enddo enddo - 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 @@ -688,7 +687,6 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) field(i, j, :), remapped_field(i, j, :), diag_cs%missing_value) 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 @@ -702,7 +700,6 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) field(i, j, :), remapped_field(i, j, :), diag_cs%missing_value) enddo enddo - else do j=diag_cs%G%jsc, diag_cs%G%jec do i=diag_cs%G%isc, diag_cs%G%iec @@ -1014,9 +1011,10 @@ function register_diag_field(module_name, field_name, axes, init_time, & endif ! Remap to z vertical coordinate, note that only diagnostics on layers - ! (not interfaces) are supported, + ! (not interfaces) are supported, also B axes are not supported yet if (get_diag_field_id_fms(module_name//'_z_new', field_name) /= DIAG_FIELD_NOT_FOUND & - .and. is_layer_axes(axes, diag_cs) .and. axes%rank == 3) then + .and. is_layer_axes(axes, diag_cs) .and. (.not. is_B_axes(axes, diag_cs)) & + .and. axes%rank == 3) then if (.not. allocated(diag_cs%zi_remap)) then call MOM_error(FATAL, 'register_diag_field: Request to regrid but no '// & 'destination grid spec provided, see param DIAG_REMAP_Z_GRID_DEF') @@ -1047,7 +1045,8 @@ function register_diag_field(module_name, field_name, axes, init_time, & posted_cmor_long_name, posted_cmor_units, & posted_cmor_standard_name) endif - if (is_layer_axes(axes, diag_cs) .and. axes%rank == 3) then + if (is_layer_axes(axes, diag_cs) .and. (.not. is_B_axes(axes, diag_cs)) & + .and. axes%rank == 3) then call log_available_diag(associated(z_remap_diag), module_name//'_z_new', field_name, & long_name, units, standard_name) endif @@ -1603,12 +1602,12 @@ function is_layer_axes(axes, diag_cs) is_layer_axes = .false. if (axes%id == diag_cs%axesTZL%id .or. & - !axes%id == diag_cs%axesBZL%id .or. & + axes%id == diag_cs%axesBZL%id .or. & axes%id == diag_cs%axesCuZL%id .or. & axes%id == diag_cs%axesCvZL%id .or. & axes%id == diag_cs%axesZL%id .or. & axes%id == diag_cs%axesTL%id .or. & - !axes%id == diag_cs%axesBL%id .or. & + axes%id == diag_cs%axesBL%id .or. & axes%id == diag_cs%axesCuL%id .or. & axes%id == diag_cs%axesCvL%id) then is_layer_axes = .true. From d7124653afe20f8ddcba96c43e998703631e85fa Mon Sep 17 00:00:00 2001 From: Nicholas Hannah Date: Mon, 13 Jul 2015 13:06:28 -0700 Subject: [PATCH 08/18] Restructuring of z star remapping to build grids less frequently. #62 --- src/framework/MOM_diag_mediator.F90 | 228 ++++++++++++++++++---------- 1 file changed, 148 insertions(+), 80 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index cf24fb1ce3..5c55ae4de8 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -139,14 +139,26 @@ module MOM_diag_mediator !default missing value to be sent to ALL diagnostics registrations real :: missing_value = 1.0e+20 - ! Interface locations for Z remapping + ! Needed to diagnostic regridding using ALE + type(remapping_CS) :: remap_cs + type(regridding_CS) :: regrid_cs + ! Nominal interface locations for Z remapping real, dimension(:), allocatable :: zi_remap - ! Layer locations for Z remapping + ! Nominal layer locations for Z remapping real, dimension(:), allocatable :: zl_remap ! Number of z levels used for remapping integer :: nz_remap - ! Pointer to H and G for Z remapping + ! 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 remapping real, dimension(:,:,:), pointer :: h => null() type(ocean_grid_type), pointer :: G => null() @@ -327,6 +339,83 @@ subroutine set_axes_info(G, param_file, diag_cs, set_vertical) 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(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(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(: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. 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 * (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(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. 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(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. 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(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 update_diag_target_grids + function check_grid_def(filename, varname) ! Do some basic checks on the vertical grid definition file, variable character(len=*), intent(in) :: filename @@ -639,119 +728,86 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) ! (in) diag - structure used to regulate diagnostic output ! (out) remapped_field - the field argument remapped to z star coordinate - type(remapping_CS) :: remap_cs - type(regridding_CS) :: regrid_cs - integer :: nz_src - integer :: i, j + real, dimension(diag_cs%nz_remap+1) :: dz real, dimension(diag_cs%nz_remap) :: h_dest - real, dimension(size(field, 3)) :: h_src - real :: depth + 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.') - call assert(allocated(diag_cs%zi_remap), 'Remapping axis not initialized') remapped_field = diag_cs%missing_value nz_src = size(field, 3) - - ! In itialise remapping and set nominal thicknesses to remap to - call initialize_regridding(diag_cs%nz_remap, 'Z*', 'PPM_IH4', regrid_cs) - call setRegriddingMinimumThickness(diag_cs%G%Angstrom, regrid_cs) - h_dest(:) = diag_cs%zi_remap(2:) - diag_cs%zi_remap(:diag_cs%nz_remap) - call setCoordinateResolution(h_dest, regrid_cs) - call initialize_remapping(nz_src, 'PPM_IH4', remap_cs) + nz_dest = diag_cs%nz_remap if (is_u_axes(diag%remap_axes, diag_cs)) then + call assert(allocated(diag_cs%zi_u), 'remap_diag_to_z: Z grid on u points not built.') 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 - h_src(:) = 0.5 * (diag_cs%h(i,j,:) + diag_cs%h(i+1,j,:)) - depth = 0.5 * (diag_cs%G%bathyT(i, j) + diag_cs%G%bathyT(i+1, j)) - call remap_column_to_z(regrid_cs, remap_cs, nz_src, diag_cs%nz_remap, & - h_src(:), depth, diag_cs%zi_remap(:), & - field(i, j, :), remapped_field(i, j, :), diag_cs%missing_value) + + 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_v), 'remap_diag_to_z: Z grid on v points not built.') 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 - h_src(:) = 0.5 * (diag_cs%h(i,j,:) + diag_cs%h(i,j+1,:)) - depth = 0.5 * (diag_cs%G%bathyT(i, j) + diag_cs%G%bathyT(i, j+1)) - call remap_column_to_z(regrid_cs, remap_cs, nz_src, diag_cs%nz_remap, & - h_src(:), depth, diag_cs%zi_remap(:), & - field(i, j, :), remapped_field(i, j, :), diag_cs%missing_value) - 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_src(:) = 0.5 * (diag_cs%h(i,j,:) + diag_cs%h(i+1,j+1,:)) - depth = 0.5 * (diag_cs%G%bathyT(i, j) + diag_cs%G%bathyT(i+1, j+1)) - call remap_column_to_z(regrid_cs, remap_cs, nz_src, diag_cs%nz_remap, & - h_src(:), depth, diag_cs%zi_remap(:), & - field(i, j, :), remapped_field(i, j, :), diag_cs%missing_value) + + 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 else + call assert(allocated(diag_cs%zi_v), 'remap_diag_to_z: Z grid on T points not built.') 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 - call remap_column_to_z(regrid_cs, remap_cs, nz_src, diag_cs%nz_remap, & - diag_cs%h(i, j, :), diag_cs%G%bathyT(i, j), diag_cs%zi_remap(:), & - field(i, j, :), remapped_field(i, j, :), diag_cs%missing_value) + 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 end subroutine remap_diag_to_z -subroutine remap_column_to_z(regrid_cs, remap_cs, nz_src, nz_dest, h, depth, & - zi_remap, field, remapped_field, missing_value) - - type(remapping_CS), intent(in) :: remap_cs - type(regridding_CS), intent(in) :: regrid_cs - integer, intent(in) :: nz_src, nz_dest - real, dimension(:), intent(in) :: h, field, zi_remap - real, intent(in) :: depth, missing_value - real, dimension(:), intent(inout) :: remapped_field - - real, dimension(nz_dest+1) :: dz, zi_tmp - real, dimension(nz_dest) :: h_dest - integer :: k - - ! Calculate thicknesses for new Z* using nominal thicknesses from - ! h_remap, current bathymetry and total thickness. - call buildGridZstarColumn(regrid_cs, nz_dest, depth, sum(h(:)), zi_tmp) - ! Calculate how much thicknesses change between source and dest grids, do - ! remapping - zi_tmp = -zi_tmp - h_dest(:) = zi_tmp(2:) - zi_tmp(:size(zi_tmp)-1) - call dzFromH1H2(nz_src, h(:), nz_dest, h_dest, dz) - call remapping_core(remap_cs, nz_src, h(:), field(:), nz_dest, dz, & - remapped_field(:)) - - ! 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_column_to_z - subroutine post_data_3d_low(diag, field, diag_cs, is_static, mask) type(diag_type), intent(in) :: diag real, intent(in) :: field(:,:,:) @@ -1034,6 +1090,15 @@ function register_diag_field(module_name, field_name, axes, init_time, & verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, & interp_method=interp_method, tile_count=tile_count) z_remap_diag%fms_diag_id = fms_id + + if (is_u_axes(axes, diag_cs)) then + diag_cs%do_z_remapping_on_u = .true. + elseif (is_v_axes(axes, diag_cs)) then + diag_cs%do_z_remapping_on_v = .true. + else + diag_cs%do_z_remapping_on_T = .true. + endif + endif ! Document diagnostics in list of available diagnostics @@ -1393,6 +1458,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 + 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) From d4bfb94ec3368086a0ad6aa5ae563869797ebf26 Mon Sep 17 00:00:00 2001 From: Nicholas Hannah Date: Tue, 14 Jul 2015 13:14:40 -0700 Subject: [PATCH 09/18] Code cleanup, correct masking of fields below bathymetry. # 62 --- src/framework/MOM_diag_mediator.F90 | 199 +++++++++++++++------------- 1 file changed, 107 insertions(+), 92 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 5c55ae4de8..fb5415b5b1 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -339,83 +339,6 @@ subroutine set_axes_info(G, param_file, diag_cs, set_vertical) 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(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(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(: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. 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 * (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(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. 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(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. 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(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 update_diag_target_grids - function check_grid_def(filename, varname) ! Do some basic checks on the vertical grid definition file, variable character(len=*), intent(in) :: filename @@ -668,7 +591,7 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask) integer, intent(in) :: diag_field_id real, intent(in) :: field(:,:,:) - type(diag_ctrl), target, intent(in) :: diag_cs + type(diag_ctrl), target, intent(inout) :: diag_cs logical, optional, intent(in) :: is_static real, optional, intent(in) :: mask(:,:,:) @@ -698,6 +621,7 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask) endif allocate(remapped_field(DIM_I(field),DIM_J(field), diag_cs%nz_remap)) + call update_diag_target_grids(diag_cs%G, diag_cs) call remap_diag_to_z(field, diag, diag_cs, remapped_field) if (associated(diag%mask3d)) then ! Since 3d masks do not vary in the vertical, just use as much as is @@ -717,6 +641,8 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask) end subroutine post_data_3d subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) +! Remap diagnostic field to z-star vertical grid. +! This uses grids generated by update_diag_target_grids() real, dimension(:,:,:), intent(in) :: field type(diag_type), intent(in) :: diag @@ -741,7 +667,8 @@ 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_u), 'remap_diag_to_z: Z grid on u points not built.') + call assert(allocated(diag_cs%zi_u), & + 'remap_diag_to_z: Z grid on u points not made.') 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 @@ -750,14 +677,13 @@ subroutine remap_diag_to_z(field, diag, diag_cs, 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, :)) + 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 + ! If their nominal depth is below the bathymetry remove. do k=1, nz_dest - if (diag_cs%zi_u(i, j, k) >= diag_cs%G%bathyT(i, j)) then + if (diag_cs%zi_remap(k) >= diag_cs%G%bathyT(i, j)) then remapped_field(i, j, k:nz_dest) = diag_cs%missing_value exit endif @@ -766,7 +692,8 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) enddo elseif (is_v_axes(diag%remap_axes, diag_cs)) then - call assert(allocated(diag_cs%zi_v), 'remap_diag_to_z: Z grid on v points not built.') + call assert(allocated(diag_cs%zi_v), & + 'remap_diag_to_z: Z grid on v points not made.') 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 @@ -775,10 +702,10 @@ subroutine remap_diag_to_z(field, diag, diag_cs, 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, :)) + 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 + if (diag_cs%zi_remap(k) >= diag_cs%G%bathyT(i, j)) then remapped_field(i, j, k:nz_dest) = diag_cs%missing_value exit endif @@ -786,7 +713,8 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) enddo enddo else - call assert(allocated(diag_cs%zi_v), 'remap_diag_to_z: Z grid on T points not built.') + call assert(allocated(diag_cs%zi_T), & + 'remap_diag_to_z: Z grid on T points not built.') 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 @@ -794,10 +722,10 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) endif 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, :)) + 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 + if (diag_cs%zi_remap(k) >= diag_cs%G%bathyT(i, j)) then remapped_field(i, j, k:nz_dest) = diag_cs%missing_value exit endif @@ -808,6 +736,92 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) end subroutine remap_diag_to_z +subroutine update_diag_target_grids(G, diag_cs) +! Build target grids for diagnostic remapping. +! +! The target grids only need to be made when sea surface +! height changes, not every time a diagnostic is posted/written. + + type(ocean_grid_type), intent(in) :: 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(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), & + 'update_diag_target_grids: Remapping axis not initialized') + + ! Initialise remapping system, on the first call + 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(: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. 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 * (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(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. 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(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. 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(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 update_diag_target_grids + + subroutine post_data_3d_low(diag, field, diag_cs, is_static, mask) type(diag_type), intent(in) :: diag real, intent(in) :: field(:,:,:) @@ -1461,6 +1475,7 @@ subroutine diag_mediator_init(G, param_file, diag_cs, err_msg) 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%remapping_initialized = .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) From ecccf84e0eb1d9f7759687e85c196210266574c1 Mon Sep 17 00:00:00 2001 From: Nicholas Hannah Date: Tue, 14 Jul 2015 15:36:34 -0700 Subject: [PATCH 10/18] Get h averaging right on u, v grids. #62 --- src/framework/MOM_diag_mediator.F90 | 21 ++++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index fb5415b5b1..f8c64fc4a8 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -68,6 +68,7 @@ module MOM_diag_mediator public register_scalar_field public defineAxes, diag_masks_set public diag_set_thickness_ptr +public update_diag_target_grids interface post_data module procedure post_data_3d, post_data_2d, post_data_0d @@ -656,6 +657,7 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) real, dimension(diag_cs%nz_remap+1) :: dz real, dimension(diag_cs%nz_remap) :: h_dest + real, dimension(size(diag_cs%h, 3)) :: h_src integer :: nz_src, nz_dest integer :: i, j, k @@ -676,14 +678,16 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) endif 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, :), & + h_src(:) = 0.5 * (diag_cs%h(i,j,:) + diag_cs%h(i+1,j,:)) + call dzFromH1H2(nz_src, h_src(:), nz_dest, h_dest(:), dz) + call remapping_core(diag_cs%remap_cs, nz_src, h_src(:), & field(i, j, :), nz_dest, dz, remapped_field(i, j, :)) ! Lower levels of the remapped data get squashed to follow bathymetry, ! If their nominal depth is below the bathymetry remove. do k=1, nz_dest - if (diag_cs%zi_remap(k) >= diag_cs%G%bathyT(i, j)) then + !if (diag_cs%zi_remap(k) >= diag_cs%G%bathyT(i, j)) then + 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 @@ -701,11 +705,13 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) endif 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, :), & + h_src(:) = 0.5 * (diag_cs%h(i,j,:) + diag_cs%h(i,j+1,:)) + call dzFromH1H2(nz_src, h_src(:), nz_dest, h_dest(:), dz) + call remapping_core(diag_cs%remap_cs, nz_src, h_src(:), & field(i, j, :), nz_dest, dz, remapped_field(i, j, :)) do k=1, nz_dest - if (diag_cs%zi_remap(k) >= diag_cs%G%bathyT(i, j)) then + !if (diag_cs%zi_remap(k) >= diag_cs%G%bathyT(i, j)) then + 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 @@ -725,7 +731,8 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) 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_remap(k) >= diag_cs%G%bathyT(i, j)) then + !if (diag_cs%zi_remap(k) >= diag_cs%G%bathyT(i, j)) then + 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 From 4941877fa2a3a07761236ab42e8c83dec126038c Mon Sep 17 00:00:00 2001 From: Nicholas Hannah Date: Wed, 15 Jul 2015 19:02:50 -0700 Subject: [PATCH 11/18] Added calls to generate/update z-star grid for diagnostic remapping purposes. This grid needs to updated everytime H changes. #62 --- src/core/MOM.F90 | 49 ++++++++---- src/core/MOM_dynamics_split_RK2.F90 | 10 ++- src/framework/MOM_diag_mediator.F90 | 78 +++++++++---------- .../lateral/MOM_mixed_layer_restrat.F90 | 9 +++ .../lateral/MOM_thickness_diffuse.F90 | 5 ++ .../vertical/MOM_bulk_mixed_layer.F90 | 6 +- .../vertical/MOM_diabatic_driver.F90 | 11 ++- 7 files changed, 107 insertions(+), 61 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index caf9e11041..aad0bb2925 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -337,7 +337,8 @@ module MOM use MOM_cpu_clock, only : CLOCK_COMPONENT, CLOCK_SUBCOMPONENT use MOM_cpu_clock, only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE use MOM_coms, only : reproducing_sum -use MOM_diag_mediator, only : diag_mediator_init, enable_averaging, diag_set_thickness_ptr +use MOM_diag_mediator, only : diag_mediator_init, enable_averaging +use MOM_diag_mediator, only : diag_set_thickness_ptr, diag_update_target_grids use MOM_diag_mediator, only : disable_averaging, post_data, safe_alloc_ptr use MOM_diag_mediator, only : register_diag_field, register_static_field use MOM_diag_mediator, only : register_scalar_field @@ -970,6 +971,10 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) call do_group_pass(CS%pass_uv_T_S_h, G%Domain) call cpu_clock_end(id_clock_pass) + ! The diag mediator may need to re-generate target grids for remmapping when + ! total thickness changes. + call diag_update_target_grids(G, h, CS%diag) + if (CS%debug) then call uchksum(u,"Post-dia first u", G, haloshift=2) call vchksum(v,"Post-dia first v", G, haloshift=2) @@ -1032,6 +1037,7 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) else dtth = dt*min(ntstep,n_max-n+1) endif + call enable_averaging(dtth,Time_local+set_time(int(floor(dtth-dt+0.5))), CS%diag) call cpu_clock_begin(id_clock_thick_diff) if (associated(CS%VarMix)) call calc_slope_functions(h, CS%tv, dt, G, CS%VarMix) @@ -1043,6 +1049,11 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) call cpu_clock_end(id_clock_pass) call disable_averaging(CS%diag) if (showCallTree) call callTree_waypoint("finished thickness_diffuse_first (step_MOM)") + + ! The diag mediator may need to re-generate target grids for remmapping when + ! total thickness changes. + call diag_update_target_grids(G, h, CS%diag) + endif endif @@ -1152,11 +1163,9 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) if (CS%useMEKE) call step_forward_MEKE(CS%MEKE, h, CS%VarMix%SN_u, CS%VarMix%SN_v, & CS%visc, dt, G, CS%MEKE_CSp) - call disable_averaging(CS%diag) call cpu_clock_end(id_clock_dynamics) - CS%dt_trans = CS%dt_trans + dt if (thermo_does_span_coupling) then do_advection = (CS%dt_trans + 0.5*dt > dt_therm) @@ -1262,7 +1271,11 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) call hchksum(CS%tv%S,"Post-ALE S", G, haloshift=1) call check_redundant("Post-ALE ", u, v, G) endif - endif + endif + + ! The diag mediator may need to re-generate target grids for remmapping when + ! total thickness changes. + call diag_update_target_grids(G, h, CS%diag) call cpu_clock_begin(id_clock_pass) call do_group_pass(CS%pass_uv_T_S_h, G%Domain) @@ -1920,10 +1933,6 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) endif call callTree_waypoint("state variables allocated (initialize_MOM)") - ! Set up a pointers h within diag mediator control structure, - ! this needs to occur _after_ CS%h has been allocated. - call diag_set_thickness_ptr(CS%h, diag) - ! Set the fields that are needed for bitwise identical restarting ! the time stepping scheme. call restart_init(G, param_file, CS%restart_CSp) @@ -1972,11 +1981,6 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) call cpu_clock_end(id_clock_MOM_init) call callTree_waypoint("returned from MOM_initialize_state() (initialize_MOM)") - ! Initialize the diagnostics mask arrays. - ! This step has to be done after call to MOM_initialize_state - ! and before MOM_diagnostics_init - call diag_masks_set(G, CS%missing, diag) - if (CS%use_ALE_algorithm) then ! For now, this has to follow immediately after MOM_initialize_state because ! the call to initialize_ALE can change CS%h, etc. initialize_ALE should @@ -2000,9 +2004,24 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) endif endif - ! This call sets up the diagnostic axes. - call cpu_clock_begin(id_clock_MOM_init) + ! Initialize the diagnostics mask arrays. + ! This step has to be done after call to MOM_initialize_state + ! and before MOM_diagnostics_init + call diag_masks_set(G, CS%missing, diag) + + ! Set up a pointers h within diag mediator control structure, + ! this needs to occur _after_ CS%h has been allocated. + call diag_set_thickness_ptr(CS%h, diag) + + ! This call sets up the diagnostic axes. These are needed, + ! e.g. to generate the target grids below. call set_axes_info(G, param_file, diag) + + ! The diag mediator may need to (re)generate target grids for remmapping when + ! total thickness changes. + call diag_update_target_grids(G, CS%h, diag) + + call cpu_clock_begin(id_clock_MOM_init) if (CS%use_ALE_algorithm) then call ALE_writeCoordinateFile( CS%ALE_CSp, G, dirs%output_directory ) endif diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 3c602a5612..466ec0d01a 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -78,7 +78,7 @@ module MOM_dynamics_split_RK2 use MOM_diag_mediator, only : diag_mediator_init, enable_averaging use MOM_diag_mediator, only : disable_averaging, post_data, safe_alloc_ptr use MOM_diag_mediator, only : register_diag_field, register_static_field -use MOM_diag_mediator, only : set_diag_mediator_grid, diag_ctrl +use MOM_diag_mediator, only : set_diag_mediator_grid, diag_ctrl, diag_update_target_grids use MOM_domains, only : MOM_domains_init use MOM_domains, only : To_South, To_West, To_All, CGRID_NE, SCALAR_PAIR use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type @@ -701,6 +701,8 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & call cpu_clock_end(id_clock_continuity) if (showCallTree) call callTree_wayPoint("done with continuity (step_MOM_dyn_split_RK2)") + call diag_update_target_grids(G, h, CS%diag) + call cpu_clock_begin(id_clock_pass) call do_group_pass(CS%pass_hp_uv, G%Domain) if (G%nonblocking_updates) then @@ -909,6 +911,10 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & call cpu_clock_end(id_clock_pass) if (showCallTree) call callTree_wayPoint("done with continuity (step_MOM_dyn_split_RK2)") + ! The diag mediator may need to re-generate target grids for remmapping when + ! total thickness changes. + call diag_update_target_grids(G, h, CS%diag) + call cpu_clock_begin(id_clock_pass) if (G%nonblocking_updates) then call start_group_pass(CS%pass_av_uvh, G%Domain) @@ -1269,6 +1275,8 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, param_file, & CS%h_av(:,:,:) = h(:,:,:) endif + call diag_update_target_grids(G, h, CS%diag) + call cpu_clock_begin(id_clock_pass_init) call create_group_pass(pass_av_h_uvh, CS%u_av,CS%v_av, G%Domain) call create_group_pass(pass_av_h_uvh, CS%h_av, G%Domain) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index f8c64fc4a8..cfe9968ea3 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -68,7 +68,7 @@ module MOM_diag_mediator public register_scalar_field public defineAxes, diag_masks_set public diag_set_thickness_ptr -public update_diag_target_grids +public diag_update_target_grids interface post_data module procedure post_data_3d, post_data_2d, post_data_0d @@ -622,7 +622,8 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask) endif allocate(remapped_field(DIM_I(field),DIM_J(field), diag_cs%nz_remap)) - call update_diag_target_grids(diag_cs%G, diag_cs) + !print*, 'post_data_3d: sum(diag_cs%h): ', sum(diag_cs%h) + !call diag_update_target_grids(diag_cs%G, diag_cs%h, diag_cs) call remap_diag_to_z(field, diag, diag_cs, remapped_field) if (associated(diag%mask3d)) then ! Since 3d masks do not vary in the vertical, just use as much as is @@ -643,7 +644,7 @@ end subroutine post_data_3d subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) ! Remap diagnostic field to z-star vertical grid. -! This uses grids generated by update_diag_target_grids() +! This uses grids generated by diag_update_target_grids() real, dimension(:,:,:), intent(in) :: field type(diag_type), intent(in) :: diag @@ -661,16 +662,16 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) integer :: nz_src, nz_dest integer :: i, j, k + call assert(diag_cs%remapping_initialized, & + 'remap_diag_to_z: Remmaping not initialized.') call assert(size(field, 3) == size(diag_cs%h, 3), & - 'Remap field and thickness z-axes do not match.') + 'remap_diag_to_z: Remap field and thickness z-axes do not match.') remapped_field = diag_cs%missing_value nz_src = size(field, 3) nz_dest = diag_cs%nz_remap if (is_u_axes(diag%remap_axes, diag_cs)) then - call assert(allocated(diag_cs%zi_u), & - 'remap_diag_to_z: Z grid on u points not made.') 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 @@ -686,8 +687,7 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) ! Lower levels of the remapped data get squashed to follow bathymetry, ! If their nominal depth is below the bathymetry remove. do k=1, nz_dest - !if (diag_cs%zi_remap(k) >= diag_cs%G%bathyT(i, j)) then - if (diag_cs%zi_u(i, j, k) >= diag_cs%G%bathyT(i, j)) then + if (diag_cs%zi_remap(k) >= diag_cs%G%bathyT(i, j)) then remapped_field(i, j, k:nz_dest) = diag_cs%missing_value exit endif @@ -696,8 +696,6 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) enddo elseif (is_v_axes(diag%remap_axes, diag_cs)) then - call assert(allocated(diag_cs%zi_v), & - 'remap_diag_to_z: Z grid on v points not made.') 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 @@ -710,8 +708,7 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) call remapping_core(diag_cs%remap_cs, nz_src, h_src(:), & field(i, j, :), nz_dest, dz, remapped_field(i, j, :)) do k=1, nz_dest - !if (diag_cs%zi_remap(k) >= diag_cs%G%bathyT(i, j)) then - if (diag_cs%zi_v(i, j, k) >= diag_cs%G%bathyT(i, j)) then + if (diag_cs%zi_remap(k) >= diag_cs%G%bathyT(i, j)) then remapped_field(i, j, k:nz_dest) = diag_cs%missing_value exit endif @@ -719,8 +716,6 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) enddo enddo else - call assert(allocated(diag_cs%zi_T), & - 'remap_diag_to_z: Z grid on T points not built.') 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 @@ -731,8 +726,7 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) 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_remap(k) >= diag_cs%G%bathyT(i, j)) then - if (diag_cs%zi_T(i, j, k) >= diag_cs%G%bathyT(i, j)) then + if (diag_cs%zi_remap(k) >= diag_cs%G%bathyT(i, j)) then remapped_field(i, j, k:nz_dest) = diag_cs%missing_value exit endif @@ -743,29 +737,33 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) end subroutine remap_diag_to_z -subroutine update_diag_target_grids(G, diag_cs) +subroutine diag_update_target_grids(G, h, diag_cs) ! Build target grids for diagnostic remapping. ! ! The target grids only need to be made when sea surface ! height changes, not every time a diagnostic is posted/written. - type(ocean_grid_type), intent(in) :: G - type(diag_ctrl), intent(inout) :: diag_cs + type(ocean_grid_type), intent(in) :: G + real, dimension(:,:,:), intent(in) :: h + type(diag_ctrl), intent(inout) :: diag_cs ! Arguments: -! (inout) G - ocean grid structure. +! (in) G - ocean grid structure. +! (in) h - a pointer to model thickness ! (inout) diag_cs - structure used to regulate diagnostic output. - real, dimension(size(diag_cs%h, 3)) :: h_src + real, dimension(size(h, 3)) :: h_src real :: depth integer :: nz_src, nz_dest integer :: i, j, k + logical :: force nz_dest = diag_cs%nz_remap - nz_src = size(diag_cs%h, 3) + nz_src = size(h, 3) + + !print*, 'diag_update_target_grids sum(h): ', sum(h) - 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 + if (.not. diag_cs%remapping_initialized) then call assert(allocated(diag_cs%zi_remap), & 'update_diag_target_grids: Remapping axis not initialized') @@ -775,17 +773,17 @@ subroutine update_diag_target_grids(G, diag_cs) call setRegriddingMinimumThickness(G%Angstrom, 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. + + allocate(diag_cs%zi_u(G%iscB:G%iecB,G%jsc:G%jec,nz_dest+1)) + allocate(diag_cs%zi_v(G%isc:G%iec,G%jscB:G%jecB,nz_dest+1)) + allocate(diag_cs%zi_T(G%isc:G%iec,G%jsc:G%jec,nz_dest+1)) endif ! Build z-star grid on u points - if (diag_cs%do_z_remapping_on_u) then - 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 + if (diag_cs%do_z_remapping_on_u .or. (.not. diag_cs%remapping_initialized)) then do j=G%jsc, G%jec do i=G%iscB, G%iecB - h_src(:) = 0.5 * (diag_cs%h(i,j,:) + diag_cs%h(i+1,j,:)) + h_src(:) = 0.5 * (h(i,j,:) + h(i+1,j,:)) depth = 0.5 * (G%bathyT(i,j) + G%bathyT(i+1,j)) call buildGridZstarColumn(diag_cs%regrid_cs, nz_dest, depth, & sum(h_src(:)), diag_cs%zi_u(i, j, :)) @@ -795,14 +793,10 @@ subroutine update_diag_target_grids(G, diag_cs) endif ! Build z-star grid on u points - if (diag_cs%do_z_remapping_on_v) then - 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 - + if (diag_cs%do_z_remapping_on_v .or. (.not. diag_cs%remapping_initialized)) then 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,:)) + h_src(:) = 0.5 * (h(i,j,:) + h(i,j+1,:)) depth = 0.5 * (G%bathyT(i, j) + G%bathyT(i, j+1)) call buildGridZstarColumn(diag_cs%regrid_cs, nz_dest, depth, & sum(h_src(:)), diag_cs%zi_v(i, j, :)) @@ -812,22 +806,19 @@ subroutine update_diag_target_grids(G, diag_cs) endif ! Build z-star grid on T points - if (diag_cs%do_z_remapping_on_T) then - 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 - + if (diag_cs%do_z_remapping_on_T .or. (.not. diag_cs%remapping_initialized)) then do j=G%jsc, G%jec do i=G%isc, G%iec call buildGridZstarColumn(diag_cs%regrid_cs, nz_dest, G%bathyT(i, j), & - sum(diag_cs%h(i, j, :)), diag_cs%zi_T(i, j, :)) + sum(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 update_diag_target_grids + diag_cs%remapping_initialized = .true. +end subroutine diag_update_target_grids subroutine post_data_3d_low(diag, field, diag_cs, is_static, mask) type(diag_type), intent(in) :: diag @@ -1478,6 +1469,7 @@ 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%h => null() diag_cs%nz_remap = -1 diag_cs%do_z_remapping_on_u = .false. diag_cs%do_z_remapping_on_v = .false. diff --git a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 index bae71a5933..d1a4b7a5a3 100644 --- a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 +++ b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 @@ -61,6 +61,7 @@ module MOM_mixed_layer_restrat use MOM_diag_mediator, only : post_data, query_averaging_enabled, diag_ctrl use MOM_diag_mediator, only : register_diag_field, safe_alloc_ptr, time_type +use MOM_diag_mediator, only : diag_update_target_grids use MOM_error_handler, only : MOM_error, FATAL, WARNING use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_forcing_type, only : forcing @@ -394,6 +395,10 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, fluxes, dt, G, CS) enddo ; enddo ; enddo !$OMP end parallel + ! The diag mediator may need to re-generate target grids for remmapping when + ! total thickness changes. + call diag_update_target_grids(G, h, CS%diag) + ! Offer fields for averaging. if (query_averaging_enabled(CS%diag)) then if (CS%id_urestrat_time > 0) call post_data(CS%id_urestrat_time, utimescale_diag, CS%diag) @@ -644,6 +649,10 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, fluxes, dt, G, CS) enddo ; enddo ; enddo !$OMP end parallel + ! The diag mediator may need to re-generate target grids for remmapping when + ! total thickness changes. + call diag_update_target_grids(G, h, CS%diag) + ! Offer fields for averaging. if (query_averaging_enabled(CS%diag) .and. & ((CS%id_urestrat_time > 0) .or. (CS%id_vrestrat_time > 0))) then diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index 6ac7db1e1e..f811040257 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -48,6 +48,7 @@ module MOM_thickness_diffuse use MOM_checksums, only : hchksum, uchksum, vchksum use MOM_diag_mediator, only : post_data, query_averaging_enabled, diag_ctrl use MOM_diag_mediator, only : register_diag_field, safe_alloc_ptr, time_type +use MOM_diag_mediator, only : diag_update_target_grids use MOM_error_handler, only : MOM_error, FATAL, WARNING use MOM_EOS, only : calculate_density, calculate_density_derivs use MOM_file_parser, only : get_param, log_version, param_file_type @@ -338,6 +339,10 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, MEKE, VarMix, CDp, CS) enddo ; enddo enddo + ! The diag mediator may need to re-generate target grids for remmapping when + ! total thickness changes. + call diag_update_target_grids(G, h, CS%diag) + if (MEKE_not_null .AND. ASSOCIATED(VarMix)) then if (ASSOCIATED(MEKE%Rd_dx_h) .and. ASSOCIATED(VarMix%Rd_dx_h)) then !$OMP parallel do default(none) shared(is,ie,js,je,MEKE,VarMix) diff --git a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 index da224979ca..2802827a81 100644 --- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 +++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 @@ -55,7 +55,7 @@ module MOM_bulk_mixed_layer use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_alloc -use MOM_diag_mediator, only : time_type, diag_ctrl +use MOM_diag_mediator, only : time_type, diag_ctrl, diag_update_target_grids use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type use MOM_error_handler, only : MOM_error, FATAL, WARNING use MOM_file_parser, only : get_param, log_param, log_version, param_file_type @@ -805,6 +805,10 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, CS, & enddo ! j loop + ! The diag mediator may need to re-generate target grids for remmapping when + ! total thickness changes. + call diag_update_target_grids(G, h_3d, CS%diag) + !$OMP end parallel if (write_diags) then diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 94892f9577..c638a4d2aa 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -72,7 +72,7 @@ module MOM_diabatic_driver use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr -use MOM_diag_mediator, only : diag_ctrl, time_type +use MOM_diag_mediator, only : diag_ctrl, time_type, diag_update_target_grids use MOM_diag_to_Z, only : diag_to_Z_CS, register_Zint_diag, calc_Zint_diags use MOM_diffConvection, only : diffConvection_CS, diffConvection_init use MOM_diffConvection, only : diffConvection_calculate, diffConvection_end @@ -408,6 +408,10 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS) if (CS%debugConservation) call MOM_state_stats('geothermal', u, v, h, tv%T, tv%S, G) endif + ! The diag mediator may need to re-generate target grids for remmapping when + ! total thickness changes. + call diag_update_target_grids(G, h, CS%diag) + ! Set_opacity estimates the optical properties of the water column. ! It will need to be modified later to include information about the ! biological properties and layer thicknesses. @@ -438,6 +442,7 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS) G, CS%bulkmixedlayer_CSp, CS%optics, & CS%aggregate_FW_forcing, dt, last_call=.true.) endif + ! Keep salinity from falling below a small but positive threshold. ! This occurs when the ice model attempts to extract more salt than ! is actually present in the ocean. @@ -902,6 +907,10 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS) if (CS%debugConservation) call MOM_state_stats('regularize_layers', u, v, h, tv%T, tv%S, G) endif + ! The diag mediator may need to re-generate target grids for remmapping when + ! total thickness changes. + call diag_update_target_grids(G, h, CS%diag) + if ((CS%id_Tdif > 0) .or. (CS%id_Tdif_z > 0) .or. & (CS%id_Tadv > 0) .or. (CS%id_Tadv_z > 0)) then do j=js,je ; do i=is,ie From 49cd415e2f56c39a0acdb8592211f5be2a2d91bf Mon Sep 17 00:00:00 2001 From: Nicholas Hannah Date: Thu, 16 Jul 2015 14:14:53 -0700 Subject: [PATCH 12/18] Code to check whether H has changed on diag post(). #62 --- src/framework/MOM_diag_mediator.F90 | 33 +++++++++++++++++++++++------ 1 file changed, 26 insertions(+), 7 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index cfe9968ea3..a4be635a58 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -163,6 +163,10 @@ module MOM_diag_mediator real, dimension(:,:,:), pointer :: h => null() type(ocean_grid_type), pointer :: G => null() + ! Keep a copy of h so that we know whether it has changed. If it has then + ! need to update the target grid for vertical remapping. + real, dimension(:,:,:), allocatable :: h_old + end type diag_ctrl integer :: doc_unit = -1 @@ -622,7 +626,6 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask) endif allocate(remapped_field(DIM_I(field),DIM_J(field), diag_cs%nz_remap)) - !print*, 'post_data_3d: sum(diag_cs%h): ', sum(diag_cs%h) !call diag_update_target_grids(diag_cs%G, diag_cs%h, diag_cs) call remap_diag_to_z(field, diag, diag_cs, remapped_field) if (associated(diag%mask3d)) then @@ -756,13 +759,11 @@ subroutine diag_update_target_grids(G, h, diag_cs) real :: depth integer :: nz_src, nz_dest integer :: i, j, k - logical :: force + logical :: force, h_changed nz_dest = diag_cs%nz_remap nz_src = size(h, 3) - !print*, 'diag_update_target_grids sum(h): ', sum(h) - if (.not. diag_cs%remapping_initialized) then call assert(allocated(diag_cs%zi_remap), & 'update_diag_target_grids: Remapping axis not initialized') @@ -779,8 +780,23 @@ subroutine diag_update_target_grids(G, h, diag_cs) allocate(diag_cs%zi_T(G%isc:G%iec,G%jsc:G%jec,nz_dest+1)) endif + ! See whether H has changed anywhere + h_changed = .true. + !do k=RANGE_K(h) + ! do j=RANGE_J(h) + ! do i=RANGE_I(h) + ! if (diag_cs%h_old(i, j, k) /= h(i, j, k)) then + ! h_changed = .true. + ! exit + ! endif + ! enddo + ! if (h_changed) exit + ! enddo + ! if (h_changed) exit + !enddo + ! Build z-star grid on u points - if (diag_cs%do_z_remapping_on_u .or. (.not. diag_cs%remapping_initialized)) then + if (h_changed .and. (diag_cs%do_z_remapping_on_u .or. (.not. diag_cs%remapping_initialized))) then do j=G%jsc, G%jec do i=G%iscB, G%iecB h_src(:) = 0.5 * (h(i,j,:) + h(i+1,j,:)) @@ -793,7 +809,7 @@ subroutine diag_update_target_grids(G, h, diag_cs) endif ! Build z-star grid on u points - if (diag_cs%do_z_remapping_on_v .or. (.not. diag_cs%remapping_initialized)) then + if (h_changed .and. (diag_cs%do_z_remapping_on_v .or. (.not. diag_cs%remapping_initialized))) then do j=G%jscB, G%jecB do i=G%isc, G%iec h_src(:) = 0.5 * (h(i,j,:) + h(i,j+1,:)) @@ -806,7 +822,7 @@ subroutine diag_update_target_grids(G, h, diag_cs) endif ! Build z-star grid on T points - if (diag_cs%do_z_remapping_on_T .or. (.not. diag_cs%remapping_initialized)) then + if (h_changed .and. (diag_cs%do_z_remapping_on_T .or. (.not. diag_cs%remapping_initialized))) then do j=G%jsc, G%jec do i=G%isc, G%iec call buildGridZstarColumn(diag_cs%regrid_cs, nz_dest, G%bathyT(i, j), & @@ -817,6 +833,7 @@ subroutine diag_update_target_grids(G, h, diag_cs) endif diag_cs%remapping_initialized = .true. + diag_cs%h_old(:,:,:) = h(:,:,:) end subroutine diag_update_target_grids @@ -1475,6 +1492,8 @@ subroutine diag_mediator_init(G, param_file, diag_cs, err_msg) diag_cs%do_z_remapping_on_v = .false. diag_cs%do_z_remapping_on_T = .false. diag_cs%remapping_initialized = .false. + allocate(diag_cs%h_old(G%isd:G%ied,G%jsd:G%jed,G%ks:G%ke)) + diag_cs%h_old(:,:,:) = 0.0 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) From bc8a1850f0477369327daa1703d424dc50f7cf3a Mon Sep 17 00:00:00 2001 From: Nicholas Hannah Date: Thu, 16 Jul 2015 15:44:24 -0700 Subject: [PATCH 13/18] Change interface of diag_update_target_grid - no need to take H as an input. #62 --- src/core/MOM.F90 | 8 +-- src/core/MOM_dynamics_split_RK2.F90 | 8 ++- src/framework/MOM_diag_mediator.F90 | 70 ++++++++++--------- .../lateral/MOM_mixed_layer_restrat.F90 | 4 +- .../lateral/MOM_thickness_diffuse.F90 | 2 +- .../vertical/MOM_bulk_mixed_layer.F90 | 2 +- .../vertical/MOM_diabatic_driver.F90 | 4 +- 7 files changed, 51 insertions(+), 47 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index aad0bb2925..2126535dd1 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -973,7 +973,7 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) ! The diag mediator may need to re-generate target grids for remmapping when ! total thickness changes. - call diag_update_target_grids(G, h, CS%diag) + call diag_update_target_grids(G, CS%diag) if (CS%debug) then call uchksum(u,"Post-dia first u", G, haloshift=2) @@ -1052,7 +1052,7 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) ! The diag mediator may need to re-generate target grids for remmapping when ! total thickness changes. - call diag_update_target_grids(G, h, CS%diag) + call diag_update_target_grids(G, CS%diag) endif endif @@ -1275,7 +1275,7 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) ! The diag mediator may need to re-generate target grids for remmapping when ! total thickness changes. - call diag_update_target_grids(G, h, CS%diag) + call diag_update_target_grids(G, CS%diag) call cpu_clock_begin(id_clock_pass) call do_group_pass(CS%pass_uv_T_S_h, G%Domain) @@ -2019,7 +2019,7 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) ! The diag mediator may need to (re)generate target grids for remmapping when ! total thickness changes. - call diag_update_target_grids(G, CS%h, diag) + call diag_update_target_grids(G, diag) call cpu_clock_begin(id_clock_MOM_init) if (CS%use_ALE_algorithm) then diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 466ec0d01a..996d4bf13f 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -701,7 +701,9 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & call cpu_clock_end(id_clock_continuity) if (showCallTree) call callTree_wayPoint("done with continuity (step_MOM_dyn_split_RK2)") - call diag_update_target_grids(G, h, CS%diag) + ! The diag mediator may need to re-generate target grids for remmapping when + ! total thickness changes. + call diag_update_target_grids(G, CS%diag) call cpu_clock_begin(id_clock_pass) call do_group_pass(CS%pass_hp_uv, G%Domain) @@ -913,7 +915,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & ! The diag mediator may need to re-generate target grids for remmapping when ! total thickness changes. - call diag_update_target_grids(G, h, CS%diag) + call diag_update_target_grids(G, CS%diag) call cpu_clock_begin(id_clock_pass) if (G%nonblocking_updates) then @@ -1275,7 +1277,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, param_file, & CS%h_av(:,:,:) = h(:,:,:) endif - call diag_update_target_grids(G, h, CS%diag) + call diag_update_target_grids(G, CS%diag) call cpu_clock_begin(id_clock_pass_init) call create_group_pass(pass_av_h_uvh, CS%u_av,CS%v_av, G%Domain) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index a4be635a58..c897ccbe5c 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -163,9 +163,11 @@ module MOM_diag_mediator real, dimension(:,:,:), pointer :: h => null() type(ocean_grid_type), pointer :: G => null() +#if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__) ! Keep a copy of h so that we know whether it has changed. If it has then - ! need to update the target grid for vertical remapping. + ! need the target grid for vertical remapping needs to have been updated. real, dimension(:,:,:), allocatable :: h_old +#endif end type diag_ctrl @@ -626,7 +628,6 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask) endif allocate(remapped_field(DIM_I(field),DIM_J(field), diag_cs%nz_remap)) - !call diag_update_target_grids(diag_cs%G, diag_cs%h, diag_cs) call remap_diag_to_z(field, diag, diag_cs, remapped_field) if (associated(diag%mask3d)) then ! Since 3d masks do not vary in the vertical, just use as much as is @@ -674,6 +675,20 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) nz_src = size(field, 3) nz_dest = diag_cs%nz_remap +#if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__) + ! Check that H is up-to-date. + do k=RANGE_K(diag_cs%h) + do j=RANGE_J(diag_cs%h) + do i=RANGE_I(diag_cs%h) + if (diag_cs%h_old(i, j, k) /= diag_cs%h(i, j, k)) then + call MOM_error(FATAL, "remap_diag_to_z: H has changed since "//& + "remapping grids were updated") + endif + enddo + enddo + enddo +#endif + if (is_u_axes(diag%remap_axes, diag_cs)) then do j=diag_cs%G%jsc, diag_cs%G%jec do i=diag_cs%G%iscB, diag_cs%G%iecB @@ -687,8 +702,8 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) call remapping_core(diag_cs%remap_cs, nz_src, h_src(:), & field(i, j, :), nz_dest, dz, remapped_field(i, j, :)) - ! Lower levels of the remapped data get squashed to follow bathymetry, - ! If their nominal depth is below the bathymetry remove. + ! Lower levels of the remapped data get squashed to follow bathymetry. + ! If their nominal depth is below the bathymetry remove them. do k=1, nz_dest if (diag_cs%zi_remap(k) >= diag_cs%G%bathyT(i, j)) then remapped_field(i, j, k:nz_dest) = diag_cs%missing_value @@ -740,29 +755,27 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) end subroutine remap_diag_to_z -subroutine diag_update_target_grids(G, h, diag_cs) -! Build target grids for diagnostic remapping. +subroutine diag_update_target_grids(G, diag_cs) +! Build/update target vertical grids for diagnostic remapping. ! -! The target grids only need to be made when sea surface -! height changes, not every time a diagnostic is posted/written. +! The target grids need to be updated whenever sea surface +! height changes. type(ocean_grid_type), intent(in) :: G - real, dimension(:,:,:), intent(in) :: h type(diag_ctrl), intent(inout) :: diag_cs ! Arguments: ! (in) G - ocean grid structure. -! (in) h - a pointer to model thickness ! (inout) diag_cs - structure used to regulate diagnostic output. - real, dimension(size(h, 3)) :: h_src + real, dimension(size(diag_cs%h, 3)) :: h_src real :: depth integer :: nz_src, nz_dest integer :: i, j, k logical :: force, h_changed nz_dest = diag_cs%nz_remap - nz_src = size(h, 3) + nz_src = size(diag_cs%h, 3) if (.not. diag_cs%remapping_initialized) then call assert(allocated(diag_cs%zi_remap), & @@ -780,26 +793,11 @@ subroutine diag_update_target_grids(G, h, diag_cs) allocate(diag_cs%zi_T(G%isc:G%iec,G%jsc:G%jec,nz_dest+1)) endif - ! See whether H has changed anywhere - h_changed = .true. - !do k=RANGE_K(h) - ! do j=RANGE_J(h) - ! do i=RANGE_I(h) - ! if (diag_cs%h_old(i, j, k) /= h(i, j, k)) then - ! h_changed = .true. - ! exit - ! endif - ! enddo - ! if (h_changed) exit - ! enddo - ! if (h_changed) exit - !enddo - ! Build z-star grid on u points - if (h_changed .and. (diag_cs%do_z_remapping_on_u .or. (.not. diag_cs%remapping_initialized))) then + if (diag_cs%do_z_remapping_on_u .or. (.not. diag_cs%remapping_initialized)) then 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(diag_cs%regrid_cs, nz_dest, depth, & sum(h_src(:)), diag_cs%zi_u(i, j, :)) @@ -809,10 +807,10 @@ subroutine diag_update_target_grids(G, h, diag_cs) endif ! Build z-star grid on u points - if (h_changed .and. (diag_cs%do_z_remapping_on_v .or. (.not. diag_cs%remapping_initialized))) then + if (diag_cs%do_z_remapping_on_v .or. (.not. diag_cs%remapping_initialized)) then do j=G%jscB, G%jecB do i=G%isc, G%iec - h_src(:) = 0.5 * (h(i,j,:) + h(i,j+1,:)) + 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(diag_cs%regrid_cs, nz_dest, depth, & sum(h_src(:)), diag_cs%zi_v(i, j, :)) @@ -822,18 +820,22 @@ subroutine diag_update_target_grids(G, h, diag_cs) endif ! Build z-star grid on T points - if (h_changed .and. (diag_cs%do_z_remapping_on_T .or. (.not. diag_cs%remapping_initialized))) then + if (diag_cs%do_z_remapping_on_T .or. (.not. diag_cs%remapping_initialized)) then do j=G%jsc, G%jec do i=G%isc, G%iec call buildGridZstarColumn(diag_cs%regrid_cs, nz_dest, G%bathyT(i, j), & - sum(h(i, j, :)), diag_cs%zi_T(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 diag_cs%remapping_initialized = .true. - diag_cs%h_old(:,:,:) = h(:,:,:) +#if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__) + ! Keep a copy of H - used to check whether grids are up-to-date + ! when doing remapping. + diag_cs%h_old(:,:,:) = diag_cs%h(:,:,:) +#endif end subroutine diag_update_target_grids diff --git a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 index d1a4b7a5a3..00064b9067 100644 --- a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 +++ b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 @@ -397,7 +397,7 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, fluxes, dt, G, CS) ! The diag mediator may need to re-generate target grids for remmapping when ! total thickness changes. - call diag_update_target_grids(G, h, CS%diag) + call diag_update_target_grids(G, CS%diag) ! Offer fields for averaging. if (query_averaging_enabled(CS%diag)) then @@ -651,7 +651,7 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, fluxes, dt, G, CS) ! The diag mediator may need to re-generate target grids for remmapping when ! total thickness changes. - call diag_update_target_grids(G, h, CS%diag) + call diag_update_target_grids(G, CS%diag) ! Offer fields for averaging. if (query_averaging_enabled(CS%diag) .and. & diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index f811040257..5e6ca27ebc 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -341,7 +341,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, MEKE, VarMix, CDp, CS) ! The diag mediator may need to re-generate target grids for remmapping when ! total thickness changes. - call diag_update_target_grids(G, h, CS%diag) + call diag_update_target_grids(G, CS%diag) if (MEKE_not_null .AND. ASSOCIATED(VarMix)) then if (ASSOCIATED(MEKE%Rd_dx_h) .and. ASSOCIATED(VarMix%Rd_dx_h)) then diff --git a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 index 2802827a81..5cdf994fcb 100644 --- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 +++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 @@ -807,7 +807,7 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, CS, & ! The diag mediator may need to re-generate target grids for remmapping when ! total thickness changes. - call diag_update_target_grids(G, h_3d, CS%diag) + call diag_update_target_grids(G, CS%diag) !$OMP end parallel diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index c638a4d2aa..06ad550324 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -410,7 +410,7 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS) ! The diag mediator may need to re-generate target grids for remmapping when ! total thickness changes. - call diag_update_target_grids(G, h, CS%diag) + call diag_update_target_grids(G, CS%diag) ! Set_opacity estimates the optical properties of the water column. ! It will need to be modified later to include information about the @@ -909,7 +909,7 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS) ! The diag mediator may need to re-generate target grids for remmapping when ! total thickness changes. - call diag_update_target_grids(G, h, CS%diag) + call diag_update_target_grids(G, CS%diag) if ((CS%id_Tdif > 0) .or. (CS%id_Tdif_z > 0) .or. & (CS%id_Tadv > 0) .or. (CS%id_Tadv_z > 0)) then From b690e3bce8458bf51a8901e7211bfbba557f9951 Mon Sep 17 00:00:00 2001 From: Nicholas Hannah Date: Thu, 16 Jul 2015 18:43:13 -0700 Subject: [PATCH 14/18] Clean up comments and interface to diag grid definition update. #62 --- src/core/MOM.F90 | 25 ++++++++++--------- src/core/MOM_continuity.F90 | 5 ++-- src/core/MOM_continuity_PPM.F90 | 9 +++++-- src/core/MOM_dynamics_legacy_split.F90 | 21 ++++++++-------- src/core/MOM_dynamics_split_RK2.F90 | 18 +++---------- src/core/MOM_dynamics_unsplit.F90 | 7 +++--- src/core/MOM_dynamics_unsplit_RK2.F90 | 7 +++--- src/framework/MOM_diag_mediator.F90 | 8 +++--- .../lateral/MOM_mixed_layer_restrat.F90 | 13 +++++----- .../lateral/MOM_thickness_diffuse.F90 | 7 +++--- .../vertical/MOM_bulk_mixed_layer.F90 | 7 +++--- .../vertical/MOM_diabatic_driver.F90 | 12 ++++----- 12 files changed, 72 insertions(+), 67 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 2126535dd1..b3a2ce4593 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -971,9 +971,9 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) call do_group_pass(CS%pass_uv_T_S_h, G%Domain) call cpu_clock_end(id_clock_pass) - ! The diag mediator may need to re-generate target grids for remmapping when - ! total thickness changes. - call diag_update_target_grids(G, CS%diag) + ! Whenever thickness changes let the diag manager know, target grids + ! for vertical remapping may need to be regenerated. + call diag_update_target_grids(CS%diag) if (CS%debug) then call uchksum(u,"Post-dia first u", G, haloshift=2) @@ -1050,9 +1050,9 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) call disable_averaging(CS%diag) if (showCallTree) call callTree_waypoint("finished thickness_diffuse_first (step_MOM)") - ! The diag mediator may need to re-generate target grids for remmapping when - ! total thickness changes. - call diag_update_target_grids(G, CS%diag) + ! Whenever thickness changes let the diag manager know, target grids + ! for vertical remapping may need to be regenerated. + call diag_update_target_grids(CS%diag) endif endif @@ -1273,9 +1273,10 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) endif endif - ! The diag mediator may need to re-generate target grids for remmapping when - ! total thickness changes. - call diag_update_target_grids(G, CS%diag) + ! Whenever thickness changes let the diag manager know, target grids + ! for vertical remapping may need to be regenerated. This needs to + ! happen after the H update and before the next post_data. + call diag_update_target_grids(CS%diag) call cpu_clock_begin(id_clock_pass) call do_group_pass(CS%pass_uv_T_S_h, G%Domain) @@ -2017,9 +2018,9 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) ! e.g. to generate the target grids below. call set_axes_info(G, param_file, diag) - ! The diag mediator may need to (re)generate target grids for remmapping when - ! total thickness changes. - call diag_update_target_grids(G, diag) + ! Whenever thickness changes let the diag manager know, target grids + ! for vertical remapping may need to be regenerated. This needs to + call diag_update_target_grids(diag) call cpu_clock_begin(id_clock_MOM_init) if (CS%use_ALE_algorithm) then diff --git a/src/core/MOM_continuity.F90 b/src/core/MOM_continuity.F90 index b0556eaf06..0f62976a4f 100644 --- a/src/core/MOM_continuity.F90 +++ b/src/core/MOM_continuity.F90 @@ -72,7 +72,7 @@ module MOM_continuity contains -subroutine continuity(u, v, hin, h, uh, vh, dt, G, CS, uhbt, vhbt, OBC, & +subroutine continuity(u, v, hin, h, uh, vh, dt, G, CS, diag_cs, uhbt, vhbt, OBC, & visc_rem_u, visc_rem_v, u_cor, v_cor, & uhbt_aux, vhbt_aux, u_cor_aux, v_cor_aux, BT_cont) real, intent(in), dimension(NIMEMB_,NJMEM_,NKMEM_) :: u @@ -84,6 +84,7 @@ subroutine continuity(u, v, hin, h, uh, vh, dt, G, CS, uhbt, vhbt, OBC, & real, intent(in) :: dt type(ocean_grid_type), intent(inout) :: G type(continuity_CS), pointer :: CS + type(diag_ctrl), intent(inout) :: diag_cs real, intent(in), optional, dimension(NIMEMB_,NJMEM_) :: uhbt real, intent(in), optional, dimension(NIMEM_,NJMEMB_) :: vhbt type(ocean_OBC_type), pointer, optional :: OBC @@ -150,7 +151,7 @@ subroutine continuity(u, v, hin, h, uh, vh, dt, G, CS, uhbt, vhbt, OBC, & " or neither.") if (CS%continuity_scheme == PPM_SCHEME) then - call continuity_PPM(u, v, hin, h, uh, vh, dt, G, CS%PPM_CSp, uhbt, vhbt, OBC, & + call continuity_PPM(u, v, hin, h, uh, vh, dt, G, CS%PPM_CSp, diag_cs, uhbt, vhbt, OBC, & visc_rem_u, visc_rem_v, u_cor, v_cor, & uhbt_aux, vhbt_aux, u_cor_aux, v_cor_aux, BT_cont) else diff --git a/src/core/MOM_continuity_PPM.F90 b/src/core/MOM_continuity_PPM.F90 index 949173d072..552aff2aad 100644 --- a/src/core/MOM_continuity_PPM.F90 +++ b/src/core/MOM_continuity_PPM.F90 @@ -44,7 +44,7 @@ module MOM_continuity_PPM !********+*********+*********+*********+*********+*********+*********+** use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE -use MOM_diag_mediator, only : time_type, diag_ctrl +use MOM_diag_mediator, only : time_type, diag_ctrl, diag_update_target_grids use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type @@ -103,7 +103,7 @@ module MOM_continuity_PPM contains -subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, CS, uhbt, vhbt, OBC, & +subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, CS, diag_cs, uhbt, vhbt, OBC, & visc_rem_u, visc_rem_v, u_cor, v_cor, & uhbt_aux, vhbt_aux, u_cor_aux, v_cor_aux, BT_cont) real, dimension(NIMEMB_,NJMEM_,NKMEM_), intent(in) :: u @@ -115,6 +115,7 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, CS, uhbt, vhbt, OBC, & real, intent(in) :: dt type(ocean_grid_type), intent(inout) :: G type(continuity_PPM_CS), pointer :: CS + type(diag_ctrl), intent(inout) :: diag_cs real, dimension(NIMEMB_,NJMEM_), intent(in), optional :: uhbt real, dimension(NIMEM_,NJMEMB_), intent(in), optional :: vhbt type(ocean_OBC_type), pointer, optional :: OBC @@ -316,6 +317,10 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, CS, uhbt, vhbt, OBC, & endif endif + ! Whenever thickness changes let the diag manager know, target grids + ! for vertical remapping may need to be regenerated. + call diag_update_target_grids(diag_cs) + end subroutine continuity_PPM subroutine zonal_mass_flux(u, h_in, uh, dt, G, CS, LB, uhbt, OBC, & diff --git a/src/core/MOM_dynamics_legacy_split.F90 b/src/core/MOM_dynamics_legacy_split.F90 index 68e4767215..20ddf14572 100644 --- a/src/core/MOM_dynamics_legacy_split.F90 +++ b/src/core/MOM_dynamics_legacy_split.F90 @@ -588,7 +588,7 @@ subroutine step_MOM_dyn_legacy_split(u, v, h, tv, visc, & if (CS%readjust_velocity) then ! Adjust the input velocites so that their transports match uhbt_out & vhbt_out. call continuity(u, v, h, hp, uh_in, vh_in, dt, G, & - CS%continuity_CSp, uhbt_in, vhbt_in, CS%OBC, & + CS%continuity_CSp, CS%diag, uhbt_in, vhbt_in, CS%OBC, & CS%visc_rem_u, CS%visc_rem_v, u_adj, v_adj, & BT_cont=CS%BT_cont) u_init => u_adj ; v_init => v_adj @@ -601,7 +601,7 @@ subroutine step_MOM_dyn_legacy_split(u, v, h, tv, visc, & CS%readjust_velocity = .false. else call continuity(u, v, h, hp, uh_in, vh_in, dt, G, & - CS%continuity_CSp, OBC=CS%OBC, BT_cont=CS%BT_cont) + CS%continuity_CSp, CS%diag, OBC=CS%OBC, BT_cont=CS%BT_cont) !### call continuity(u, v, h, hp, uh_in, vh_in, dt, G, & !### CS%continuity_CSp, OBC=CS%OBC, visc_rem_u=CS%visc_rem_u, & !### visc_rem_v=CS%visc_rem_v, BT_cont=CS%BT_cont) @@ -631,8 +631,9 @@ subroutine step_MOM_dyn_legacy_split(u, v, h, tv, visc, & if (associated(CS%BT_cont) .or. CS%BT_use_layer_fluxes) then call cpu_clock_begin(id_clock_continuity) call continuity(u, v, h, hp, uh_in, vh_in, dt, G, & - CS%continuity_CSp, OBC=CS%OBC, visc_rem_u=CS%visc_rem_u, & - visc_rem_v=CS%visc_rem_v, BT_cont=CS%BT_cont) + CS%continuity_CSp, CS%diag, OBC=CS%OBC, & + visc_rem_u=CS%visc_rem_u, visc_rem_v=CS%visc_rem_v, & + BT_cont=CS%BT_cont) call cpu_clock_end(id_clock_continuity) if (BT_cont_BT_thick) then call cpu_clock_begin(id_clock_pass) @@ -716,8 +717,8 @@ subroutine step_MOM_dyn_legacy_split(u, v, h, tv, visc, & ! hp = h + dt * div . uh call cpu_clock_begin(id_clock_continuity) call continuity(up, vp, h, hp, uh, vh, dt, G, CS%continuity_CSp, & - CS%uhbt, CS%vhbt, CS%OBC, CS%visc_rem_u, CS%visc_rem_v, & - u_av, v_av, BT_cont=CS%BT_cont) + CS%diag, CS%uhbt, CS%vhbt, CS%OBC, CS%visc_rem_u, & + CS%visc_rem_v, u_av, v_av, BT_cont=CS%BT_cont) call cpu_clock_end(id_clock_continuity) call cpu_clock_begin(id_clock_pass) @@ -927,7 +928,7 @@ subroutine step_MOM_dyn_legacy_split(u, v, h, tv, visc, & enddo ; enddo ; enddo ; endif call cpu_clock_begin(id_clock_continuity) call continuity(u, v, h, h, uh, vh, dt, G, & - CS%continuity_CSp, CS%uhbt, CS%vhbt, CS%OBC, & + CS%continuity_CSp, CS%diag, CS%uhbt, CS%vhbt, CS%OBC, & CS%visc_rem_u, CS%visc_rem_v, u_av, v_av, & uhbt_out, vhbt_out, u, v) call cpu_clock_end(id_clock_continuity) @@ -966,7 +967,7 @@ subroutine step_MOM_dyn_legacy_split(u, v, h, tv, visc, & ! u_av and v_av adjusted so their mass transports match uhbt and vhbt. call cpu_clock_begin(id_clock_continuity) call continuity(u, v, h, h, uh, vh, dt, G, & - CS%continuity_CSp, CS%uhbt, CS%vhbt, CS%OBC, & + CS%continuity_CSp, CS%diag, CS%uhbt, CS%vhbt, CS%OBC, & CS%visc_rem_u, CS%visc_rem_v, u_av, v_av) call cpu_clock_end(id_clock_continuity) call cpu_clock_begin(id_clock_pass) @@ -1066,7 +1067,7 @@ subroutine adjustments_dyn_legacy_split(u, v, h, dt, G, CS) if (CS%readjust_BT_trans) then call cpu_clock_begin(id_clock_continuity) call continuity(u, v, h, h_temp, uh_temp, vh_temp, dt, G, & - CS%continuity_CSp, OBC=CS%OBC) + CS%continuity_CSp, CS%diag, OBC=CS%OBC) call cpu_clock_end(id_clock_continuity) !$OMP parallel default(none) shared(is,ie,js,je,nz,CS,uh_temp,vh_temp) !$OMP do @@ -1406,7 +1407,7 @@ subroutine initialize_dyn_legacy_split(u, v, h, uh, vh, eta, Time, G, param_file if (.not. query_initialized(uh,"uh",restart_CS) .or. & .not. query_initialized(vh,"vh",restart_CS)) then h_tmp(:,:,:) = h(:,:,:) - call continuity(u, v, h, h_tmp, uh, vh, dt, G, CS%continuity_CSp, OBC=CS%OBC) + call continuity(u, v, h, h_tmp, uh, vh, dt, G, CS%continuity_CSp, CS%diag, OBC=CS%OBC) call cpu_clock_begin(id_clock_pass_init) call pass_var(h_tmp, G%Domain) call cpu_clock_end(id_clock_pass_init) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 996d4bf13f..39f9d86e3b 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -605,7 +605,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & if (associated(CS%BT_cont) .or. CS%BT_use_layer_fluxes) then call cpu_clock_begin(id_clock_continuity) call continuity(u, v, h, hp, uh_in, vh_in, dt, G, & - CS%continuity_CSp, OBC=CS%OBC, visc_rem_u=CS%visc_rem_u, & + CS%continuity_CSp, CS%diag, OBC=CS%OBC, visc_rem_u=CS%visc_rem_u, & visc_rem_v=CS%visc_rem_v, BT_cont=CS%BT_cont) call cpu_clock_end(id_clock_continuity) if (BT_cont_BT_thick) then @@ -695,16 +695,12 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & ! uh = u_av * h ! hp = h + dt * div . uh call cpu_clock_begin(id_clock_continuity) - call continuity(up, vp, h, hp, uh, vh, dt, G, CS%continuity_CSp, & + call continuity(up, vp, h, hp, uh, vh, dt, G, CS%continuity_CSp, CS%diag, & CS%uhbt, CS%vhbt, CS%OBC, CS%visc_rem_u, CS%visc_rem_v, & u_av, v_av, BT_cont=CS%BT_cont) call cpu_clock_end(id_clock_continuity) if (showCallTree) call callTree_wayPoint("done with continuity (step_MOM_dyn_split_RK2)") - ! The diag mediator may need to re-generate target grids for remmapping when - ! total thickness changes. - call diag_update_target_grids(G, CS%diag) - call cpu_clock_begin(id_clock_pass) call do_group_pass(CS%pass_hp_uv, G%Domain) if (G%nonblocking_updates) then @@ -905,7 +901,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & ! u_av and v_av adjusted so their mass transports match uhbt and vhbt. call cpu_clock_begin(id_clock_continuity) call continuity(u, v, h, h, uh, vh, dt, G, & - CS%continuity_CSp, CS%uhbt, CS%vhbt, CS%OBC, & + CS%continuity_CSp, CS%diag, CS%uhbt, CS%vhbt, CS%OBC, & CS%visc_rem_u, CS%visc_rem_v, u_av, v_av) call cpu_clock_end(id_clock_continuity) call cpu_clock_begin(id_clock_pass) @@ -913,10 +909,6 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & call cpu_clock_end(id_clock_pass) if (showCallTree) call callTree_wayPoint("done with continuity (step_MOM_dyn_split_RK2)") - ! The diag mediator may need to re-generate target grids for remmapping when - ! total thickness changes. - call diag_update_target_grids(G, CS%diag) - call cpu_clock_begin(id_clock_pass) if (G%nonblocking_updates) then call start_group_pass(CS%pass_av_uvh, G%Domain) @@ -1266,7 +1258,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, param_file, & if (.not. query_initialized(uh,"uh",restart_CS) .or. & .not. query_initialized(vh,"vh",restart_CS)) then h_tmp(:,:,:) = h(:,:,:) - call continuity(u, v, h, h_tmp, uh, vh, dt, G, CS%continuity_CSp, OBC=CS%OBC) + call continuity(u, v, h, h_tmp, uh, vh, dt, G, CS%continuity_CSp, CS%diag, OBC=CS%OBC) call cpu_clock_begin(id_clock_pass_init) call create_group_pass(pass_h_tmp, h_tmp, G%Domain) call do_group_pass(pass_h_tmp, G%Domain) @@ -1277,8 +1269,6 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, param_file, & CS%h_av(:,:,:) = h(:,:,:) endif - call diag_update_target_grids(G, CS%diag) - call cpu_clock_begin(id_clock_pass_init) call create_group_pass(pass_av_h_uvh, CS%u_av,CS%v_av, G%Domain) call create_group_pass(pass_av_h_uvh, CS%h_av, G%Domain) diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index dd36f3aa19..9c51432819 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -268,7 +268,8 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, fluxes, & ! uh = u*h ! hp = h + dt/2 div . uh call cpu_clock_begin(id_clock_continuity) - call continuity(u, v, h, hp, uh, vh, dt*0.5, G, CS%continuity_CSp, OBC=CS%OBC) + call continuity(u, v, h, hp, uh, vh, dt*0.5, G, CS%continuity_CSp, CS%diag, & + OBC=CS%OBC) call cpu_clock_end(id_clock_continuity) call cpu_clock_begin(id_clock_pass) call pass_var(hp, G%Domain) @@ -377,7 +378,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, fluxes, & ! h_av = hp + dt/2 div . uh call cpu_clock_begin(id_clock_continuity) call continuity(up, vp, hp, h_av, uh, vh, & - (0.5*dt), G, CS%continuity_CSp, OBC=CS%OBC) + (0.5*dt), G, CS%continuity_CSp, CS%diag, OBC=CS%OBC) call cpu_clock_end(id_clock_continuity) call cpu_clock_begin(id_clock_pass) call pass_var(h_av, G%Domain) @@ -436,7 +437,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, fluxes, & ! h = hp + dt/2 div . uh call cpu_clock_begin(id_clock_continuity) call continuity(upp, vpp, hp, h, uh, vh, & - (dt*0.5), G, CS%continuity_CSp, OBC=CS%OBC) + (dt*0.5), G, CS%continuity_CSp, CS%diag, OBC=CS%OBC) call cpu_clock_end(id_clock_continuity) call cpu_clock_begin(id_clock_pass) call pass_var(h, G%Domain) diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90 index 57edff0a20..87a8b2052f 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/src/core/MOM_dynamics_unsplit_RK2.F90 @@ -282,7 +282,8 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, call cpu_clock_begin(id_clock_continuity) ! This is a duplicate caclulation of the last continuity from the previous step ! and could/should be optimized out. -AJA - call continuity(u_in, v_in, h_in, hp, uh, vh, dt_pred, G, CS%continuity_CSp, OBC=CS%OBC) + call continuity(u_in, v_in, h_in, hp, uh, vh, dt_pred, G, CS%continuity_CSp, & + CS%diag, OBC=CS%OBC) call cpu_clock_end(id_clock_continuity) call cpu_clock_begin(id_clock_pass) call pass_var(hp, G%Domain) @@ -369,7 +370,7 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, ! h_av = h + dt div . uh call cpu_clock_begin(id_clock_continuity) call continuity(up, vp, h_in, hp, uh, vh, & - dt, G, CS%continuity_CSp, OBC=CS%OBC) + dt, G, CS%continuity_CSp, CS%diag, OBC=CS%OBC) call cpu_clock_end(id_clock_continuity) call cpu_clock_begin(id_clock_pass) call pass_var(hp, G%Domain) @@ -426,7 +427,7 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, ! h[n+1] = h[n] + dt div . uh call cpu_clock_begin(id_clock_continuity) call continuity(up, vp, h_in, h_in, uh, vh, & - dt, G, CS%continuity_CSp, OBC=CS%OBC) + dt, G, CS%continuity_CSp, CS%diag, OBC=CS%OBC) call cpu_clock_end(id_clock_continuity) call cpu_clock_begin(id_clock_pass) call pass_var(h_in, G%Domain) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index c897ccbe5c..bedd86ab14 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -58,6 +58,8 @@ module MOM_diag_mediator #define DIM_J(a) lbound(a, 2):ubound(a, 2) #define DIM_K(a) lbound(a, 3):ubound(a, 3) +#define __DO_SAFETY_CHECKS__ + public set_axes_info, post_data, register_diag_field, time_type public post_data_1d_k public safe_alloc_ptr, safe_alloc_alloc @@ -755,20 +757,19 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) end subroutine remap_diag_to_z -subroutine diag_update_target_grids(G, diag_cs) +subroutine diag_update_target_grids(diag_cs) ! Build/update target vertical grids for diagnostic remapping. ! ! The target grids need to be updated whenever sea surface ! height changes. - type(ocean_grid_type), intent(in) :: G type(diag_ctrl), intent(inout) :: diag_cs ! Arguments: -! (in) G - ocean grid structure. ! (inout) diag_cs - structure used to regulate diagnostic output. real, dimension(size(diag_cs%h, 3)) :: h_src + type(ocean_grid_type), pointer :: G real :: depth integer :: nz_src, nz_dest integer :: i, j, k @@ -776,6 +777,7 @@ subroutine diag_update_target_grids(G, diag_cs) nz_dest = diag_cs%nz_remap nz_src = size(diag_cs%h, 3) + G => diag_cs%G if (.not. diag_cs%remapping_initialized) then call assert(allocated(diag_cs%zi_remap), & diff --git a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 index 00064b9067..41286a7a87 100644 --- a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 +++ b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 @@ -395,9 +395,10 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, fluxes, dt, G, CS) enddo ; enddo ; enddo !$OMP end parallel - ! The diag mediator may need to re-generate target grids for remmapping when - ! total thickness changes. - call diag_update_target_grids(G, CS%diag) + ! Whenever thickness changes let the diag manager know, target grids + ! for vertical remapping may need to be regenerated. + ! This needs to happen after the H update and before the next post_data. + call diag_update_target_grids(CS%diag) ! Offer fields for averaging. if (query_averaging_enabled(CS%diag)) then @@ -649,9 +650,9 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, fluxes, dt, G, CS) enddo ; enddo ; enddo !$OMP end parallel - ! The diag mediator may need to re-generate target grids for remmapping when - ! total thickness changes. - call diag_update_target_grids(G, CS%diag) + ! Whenever thickness changes let the diag manager know, target grids + ! for vertical remapping may need to be regenerated. + call diag_update_target_grids(CS%diag) ! Offer fields for averaging. if (query_averaging_enabled(CS%diag) .and. & diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index 5e6ca27ebc..680134381a 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -339,9 +339,10 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, MEKE, VarMix, CDp, CS) enddo ; enddo enddo - ! The diag mediator may need to re-generate target grids for remmapping when - ! total thickness changes. - call diag_update_target_grids(G, CS%diag) + ! Whenever thickness changes let the diag manager know, target grids + ! for vertical remapping may need to be regenerated. + ! This needs to happen after the H update and before the next post_data. + call diag_update_target_grids(CS%diag) if (MEKE_not_null .AND. ASSOCIATED(VarMix)) then if (ASSOCIATED(MEKE%Rd_dx_h) .and. ASSOCIATED(VarMix%Rd_dx_h)) then diff --git a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 index 5cdf994fcb..a1313727b3 100644 --- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 +++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 @@ -805,9 +805,10 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, CS, & enddo ! j loop - ! The diag mediator may need to re-generate target grids for remmapping when - ! total thickness changes. - call diag_update_target_grids(G, CS%diag) + ! Whenever thickness changes let the diag manager know, target grids + ! for vertical remapping may need to be regenerated. + ! This needs to happen after the H update and before the next post_data. + call diag_update_target_grids(CS%diag) !$OMP end parallel diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 06ad550324..157a0119ec 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -408,9 +408,9 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS) if (CS%debugConservation) call MOM_state_stats('geothermal', u, v, h, tv%T, tv%S, G) endif - ! The diag mediator may need to re-generate target grids for remmapping when - ! total thickness changes. - call diag_update_target_grids(G, CS%diag) + ! Whenever thickness changes let the diag manager know, target grids + ! for vertical remapping may need to be regenerated. + call diag_update_target_grids(CS%diag) ! Set_opacity estimates the optical properties of the water column. ! It will need to be modified later to include information about the @@ -907,9 +907,9 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS) if (CS%debugConservation) call MOM_state_stats('regularize_layers', u, v, h, tv%T, tv%S, G) endif - ! The diag mediator may need to re-generate target grids for remmapping when - ! total thickness changes. - call diag_update_target_grids(G, CS%diag) + ! Whenever thickness changes let the diag manager know, target grids + ! for vertical remapping may need to be regenerated. + call diag_update_target_grids(CS%diag) if ((CS%id_Tdif > 0) .or. (CS%id_Tdif_z > 0) .or. & (CS%id_Tadv > 0) .or. (CS%id_Tadv_z > 0)) then From d485907cec2204bcd7cd0f1569dd26c668a58350 Mon Sep 17 00:00:00 2001 From: Nicholas Hannah Date: Fri, 17 Jul 2015 08:38:31 -0700 Subject: [PATCH 15/18] Diag mediator deallocates memory at end(). #62 --- config_src/coupled_driver/ocean_model_MOM.F90 | 5 +-- .../ice_solo_driver/ice_shelf_driver.F90 | 2 +- config_src/solo_driver/MOM_driver.F90 | 2 +- src/framework/MOM_diag_mediator.F90 | 31 ++++++++++++++++--- 4 files changed, 32 insertions(+), 8 deletions(-) diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index 2e317fc09f..136b64e920 100644 --- a/config_src/coupled_driver/ocean_model_MOM.F90 +++ b/config_src/coupled_driver/ocean_model_MOM.F90 @@ -37,7 +37,7 @@ module ocean_model_mod use MOM, only : initialize_MOM, step_MOM, MOM_control_struct, MOM_end use MOM, only : calculate_surface_state use MOM_constants, only : CELSIUS_KELVIN_OFFSET -use MOM_diag_mediator, only : enable_averaging, disable_averaging +use MOM_diag_mediator, only : diag_ctrl, enable_averaging, disable_averaging use MOM_diag_mediator, only : diag_mediator_close_registration, diag_mediator_end use MOM_domains, only : pass_vector, AGRID, BGRID_NE, CGRID_NE use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe @@ -494,6 +494,7 @@ subroutine ocean_model_end(Ocean_sfc, Ocean_state, Time) type(ocean_public_type), intent(inout) :: Ocean_sfc type(ocean_state_type), pointer :: Ocean_state type(time_type), intent(in) :: Time + ! This subroutine terminates the model run, saving the ocean state in a ! restart file and deallocating any data associated with the ocean. @@ -504,7 +505,7 @@ subroutine ocean_model_end(Ocean_sfc, Ocean_state, Time) ! (in) Time - The model time, used for writing restarts. call ocean_model_save_restart(Ocean_state, Time) - call diag_mediator_end(Time) + call diag_mediator_end(Time, Ocean_state%MOM_CSp%diag) call MOM_end(Ocean_state%MOM_CSp) if (Ocean_state%use_ice_shelf) call ice_shelf_end(Ocean_state%Ice_shelf_CSp) end subroutine ocean_model_end diff --git a/config_src/ice_solo_driver/ice_shelf_driver.F90 b/config_src/ice_solo_driver/ice_shelf_driver.F90 index f00071f57b..2e694a2b0a 100644 --- a/config_src/ice_solo_driver/ice_shelf_driver.F90 +++ b/config_src/ice_solo_driver/ice_shelf_driver.F90 @@ -416,7 +416,7 @@ program SHELF_main close(unit) endif - call diag_mediator_end(Time, end_diag_manager=.true.) + call diag_mediator_end(Time, ice_shelf_CSp%diag, end_diag_manager=.true.) call cpu_clock_end(termClock) call io_infra_end ; call MOM_infra_end diff --git a/config_src/solo_driver/MOM_driver.F90 b/config_src/solo_driver/MOM_driver.F90 index 443188fb74..bc17b20091 100644 --- a/config_src/solo_driver/MOM_driver.F90 +++ b/config_src/solo_driver/MOM_driver.F90 @@ -513,7 +513,7 @@ program MOM_main endif call callTree_waypoint("End MOM_main") - call diag_mediator_end(Time, end_diag_manager=.true.) + call diag_mediator_end(Time, MOM_CSp%diag, end_diag_manager=.true.) call cpu_clock_end(termClock) call io_infra_end ; call MOM_infra_end diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index bedd86ab14..ee22320891 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -38,7 +38,8 @@ module MOM_diag_mediator use MOM_string_functions, only : lowercase use MOM_time_manager, only : time_type use MOM_remapping, only : remapping_CS, remapping_core, initialize_remapping, dzFromH1H2 -use MOM_regridding, only : regridding_CS, initialize_regridding, setCoordinateResolution, buildGridZStarColumn, setRegriddingMinimumThickness +use MOM_regridding, only : regridding_CS, initialize_regridding, setCoordinateResolution +use MOM_regridding, only : buildGridZStarColumn, setRegriddingMinimumThickness use diag_manager_mod, only : diag_manager_init, diag_manager_end use diag_manager_mod, only : send_data, diag_axis_init @@ -1596,14 +1597,36 @@ subroutine diag_mediator_close_registration() end subroutine diag_mediator_close_registration -subroutine diag_mediator_end(time, end_diag_manager) - type(time_type), intent(in) :: time - logical, optional, intent(in) :: end_diag_manager !< If true, call diag_manager_end() +subroutine diag_mediator_end(time, diag_cs, end_diag_manager) + type(time_type), intent(in) :: time + type(diag_ctrl), intent(inout) :: diag_cs + logical, optional, intent(in) :: end_diag_manager !< If true, call diag_manager_end() if (doc_unit > -1) then close(doc_unit) ; doc_unit = -3 endif + deallocate(diag_cs%diags) + + if (allocated(diag_cs%zi_remap)) deallocate(diag_cs%zi_remap) + if (allocated(diag_cs%zl_remap)) deallocate(diag_cs%zl_remap) + + if (allocated(diag_cs%zi_u)) deallocate(diag_cs%zi_u) + if (allocated(diag_cs%zi_v)) deallocate(diag_cs%zi_v) + if (allocated(diag_cs%zi_T)) deallocate(diag_cs%zi_T) + deallocate(diag_cs%mask3dTL) + deallocate(diag_cs%mask3dBuL) + deallocate(diag_cs%mask3dCuL) + deallocate(diag_cs%mask3dCvL) + deallocate(diag_cs%mask3dTi) + deallocate(diag_cs%mask3dBui) + deallocate(diag_cs%mask3dCui) + deallocate(diag_cs%mask3dCvi) + +#if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__) + deallocate(diag_cs%h_old) +#endif + if (present(end_diag_manager)) then if (end_diag_manager) call diag_manager_end(time) endif From 8008eea25707354fc2eacb1f99e57ddcb4b5679b Mon Sep 17 00:00:00 2001 From: Nicholas Hannah Date: Wed, 22 Jul 2015 12:59:53 -0700 Subject: [PATCH 16/18] Move calls to update grid for z remapping of diagnostics out of continuity(). It is not necessary to pass the diag mediator structure into calls to continuity(). #62 --- src/core/MOM_continuity.F90 | 5 ++--- src/core/MOM_continuity_PPM.F90 | 9 ++------- src/core/MOM_dynamics_legacy_split.F90 | 24 +++++++++++++++--------- src/core/MOM_dynamics_split_RK2.F90 | 11 +++++++---- src/core/MOM_dynamics_unsplit.F90 | 11 +++++++---- src/core/MOM_dynamics_unsplit_RK2.F90 | 6 +++--- 6 files changed, 36 insertions(+), 30 deletions(-) diff --git a/src/core/MOM_continuity.F90 b/src/core/MOM_continuity.F90 index 0f62976a4f..b0556eaf06 100644 --- a/src/core/MOM_continuity.F90 +++ b/src/core/MOM_continuity.F90 @@ -72,7 +72,7 @@ module MOM_continuity contains -subroutine continuity(u, v, hin, h, uh, vh, dt, G, CS, diag_cs, uhbt, vhbt, OBC, & +subroutine continuity(u, v, hin, h, uh, vh, dt, G, CS, uhbt, vhbt, OBC, & visc_rem_u, visc_rem_v, u_cor, v_cor, & uhbt_aux, vhbt_aux, u_cor_aux, v_cor_aux, BT_cont) real, intent(in), dimension(NIMEMB_,NJMEM_,NKMEM_) :: u @@ -84,7 +84,6 @@ subroutine continuity(u, v, hin, h, uh, vh, dt, G, CS, diag_cs, uhbt, vhbt, OBC, real, intent(in) :: dt type(ocean_grid_type), intent(inout) :: G type(continuity_CS), pointer :: CS - type(diag_ctrl), intent(inout) :: diag_cs real, intent(in), optional, dimension(NIMEMB_,NJMEM_) :: uhbt real, intent(in), optional, dimension(NIMEM_,NJMEMB_) :: vhbt type(ocean_OBC_type), pointer, optional :: OBC @@ -151,7 +150,7 @@ subroutine continuity(u, v, hin, h, uh, vh, dt, G, CS, diag_cs, uhbt, vhbt, OBC, " or neither.") if (CS%continuity_scheme == PPM_SCHEME) then - call continuity_PPM(u, v, hin, h, uh, vh, dt, G, CS%PPM_CSp, diag_cs, uhbt, vhbt, OBC, & + call continuity_PPM(u, v, hin, h, uh, vh, dt, G, CS%PPM_CSp, uhbt, vhbt, OBC, & visc_rem_u, visc_rem_v, u_cor, v_cor, & uhbt_aux, vhbt_aux, u_cor_aux, v_cor_aux, BT_cont) else diff --git a/src/core/MOM_continuity_PPM.F90 b/src/core/MOM_continuity_PPM.F90 index 552aff2aad..949173d072 100644 --- a/src/core/MOM_continuity_PPM.F90 +++ b/src/core/MOM_continuity_PPM.F90 @@ -44,7 +44,7 @@ module MOM_continuity_PPM !********+*********+*********+*********+*********+*********+*********+** use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE -use MOM_diag_mediator, only : time_type, diag_ctrl, diag_update_target_grids +use MOM_diag_mediator, only : time_type, diag_ctrl use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type @@ -103,7 +103,7 @@ module MOM_continuity_PPM contains -subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, CS, diag_cs, uhbt, vhbt, OBC, & +subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, CS, uhbt, vhbt, OBC, & visc_rem_u, visc_rem_v, u_cor, v_cor, & uhbt_aux, vhbt_aux, u_cor_aux, v_cor_aux, BT_cont) real, dimension(NIMEMB_,NJMEM_,NKMEM_), intent(in) :: u @@ -115,7 +115,6 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, CS, diag_cs, uhbt, vhbt, real, intent(in) :: dt type(ocean_grid_type), intent(inout) :: G type(continuity_PPM_CS), pointer :: CS - type(diag_ctrl), intent(inout) :: diag_cs real, dimension(NIMEMB_,NJMEM_), intent(in), optional :: uhbt real, dimension(NIMEM_,NJMEMB_), intent(in), optional :: vhbt type(ocean_OBC_type), pointer, optional :: OBC @@ -317,10 +316,6 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, CS, diag_cs, uhbt, vhbt, endif endif - ! Whenever thickness changes let the diag manager know, target grids - ! for vertical remapping may need to be regenerated. - call diag_update_target_grids(diag_cs) - end subroutine continuity_PPM subroutine zonal_mass_flux(u, h_in, uh, dt, G, CS, LB, uhbt, OBC, & diff --git a/src/core/MOM_dynamics_legacy_split.F90 b/src/core/MOM_dynamics_legacy_split.F90 index 20ddf14572..533c8c9248 100644 --- a/src/core/MOM_dynamics_legacy_split.F90 +++ b/src/core/MOM_dynamics_legacy_split.F90 @@ -78,7 +78,7 @@ module MOM_dynamics_legacy_split use MOM_diag_mediator, only : diag_mediator_init, enable_averaging use MOM_diag_mediator, only : disable_averaging, post_data, safe_alloc_ptr use MOM_diag_mediator, only : register_diag_field, register_static_field -use MOM_diag_mediator, only : set_diag_mediator_grid, diag_ctrl +use MOM_diag_mediator, only : set_diag_mediator_grid, diag_ctrl, diag_update_target_grids use MOM_domains, only : MOM_domains_init, pass_var, pass_vector use MOM_domains, only : pass_var_start, pass_var_complete use MOM_domains, only : pass_vector_start, pass_vector_complete @@ -588,7 +588,7 @@ subroutine step_MOM_dyn_legacy_split(u, v, h, tv, visc, & if (CS%readjust_velocity) then ! Adjust the input velocites so that their transports match uhbt_out & vhbt_out. call continuity(u, v, h, hp, uh_in, vh_in, dt, G, & - CS%continuity_CSp, CS%diag, uhbt_in, vhbt_in, CS%OBC, & + CS%continuity_CSp, uhbt_in, vhbt_in, CS%OBC, & CS%visc_rem_u, CS%visc_rem_v, u_adj, v_adj, & BT_cont=CS%BT_cont) u_init => u_adj ; v_init => v_adj @@ -601,7 +601,7 @@ subroutine step_MOM_dyn_legacy_split(u, v, h, tv, visc, & CS%readjust_velocity = .false. else call continuity(u, v, h, hp, uh_in, vh_in, dt, G, & - CS%continuity_CSp, CS%diag, OBC=CS%OBC, BT_cont=CS%BT_cont) + CS%continuity_CSp, OBC=CS%OBC, BT_cont=CS%BT_cont) !### call continuity(u, v, h, hp, uh_in, vh_in, dt, G, & !### CS%continuity_CSp, OBC=CS%OBC, visc_rem_u=CS%visc_rem_u, & !### visc_rem_v=CS%visc_rem_v, BT_cont=CS%BT_cont) @@ -631,7 +631,7 @@ subroutine step_MOM_dyn_legacy_split(u, v, h, tv, visc, & if (associated(CS%BT_cont) .or. CS%BT_use_layer_fluxes) then call cpu_clock_begin(id_clock_continuity) call continuity(u, v, h, hp, uh_in, vh_in, dt, G, & - CS%continuity_CSp, CS%diag, OBC=CS%OBC, & + CS%continuity_CSp, OBC=CS%OBC, & visc_rem_u=CS%visc_rem_u, visc_rem_v=CS%visc_rem_v, & BT_cont=CS%BT_cont) call cpu_clock_end(id_clock_continuity) @@ -717,7 +717,7 @@ subroutine step_MOM_dyn_legacy_split(u, v, h, tv, visc, & ! hp = h + dt * div . uh call cpu_clock_begin(id_clock_continuity) call continuity(up, vp, h, hp, uh, vh, dt, G, CS%continuity_CSp, & - CS%diag, CS%uhbt, CS%vhbt, CS%OBC, CS%visc_rem_u, & + CS%uhbt, CS%vhbt, CS%OBC, CS%visc_rem_u, & CS%visc_rem_v, u_av, v_av, BT_cont=CS%BT_cont) call cpu_clock_end(id_clock_continuity) @@ -928,10 +928,13 @@ subroutine step_MOM_dyn_legacy_split(u, v, h, tv, visc, & enddo ; enddo ; enddo ; endif call cpu_clock_begin(id_clock_continuity) call continuity(u, v, h, h, uh, vh, dt, G, & - CS%continuity_CSp, CS%diag, CS%uhbt, CS%vhbt, CS%OBC, & + CS%continuity_CSp, CS%uhbt, CS%vhbt, CS%OBC, & CS%visc_rem_u, CS%visc_rem_v, u_av, v_av, & uhbt_out, vhbt_out, u, v) call cpu_clock_end(id_clock_continuity) + ! Whenever thickness changes let the diag manager know, target grids + ! for vertical remapping may need to be regenerated. + call diag_update_target_grids(CS%diag) if (G%nonblocking_updates) then call cpu_clock_begin(id_clock_pass) pid_h = pass_var_start(h, G%Domain) @@ -967,9 +970,12 @@ subroutine step_MOM_dyn_legacy_split(u, v, h, tv, visc, & ! u_av and v_av adjusted so their mass transports match uhbt and vhbt. call cpu_clock_begin(id_clock_continuity) call continuity(u, v, h, h, uh, vh, dt, G, & - CS%continuity_CSp, CS%diag, CS%uhbt, CS%vhbt, CS%OBC, & + CS%continuity_CSp, CS%uhbt, CS%vhbt, CS%OBC, & CS%visc_rem_u, CS%visc_rem_v, u_av, v_av) call cpu_clock_end(id_clock_continuity) + ! Whenever thickness changes let the diag manager know, target grids + ! for vertical remapping may need to be regenerated. + call diag_update_target_grids(CS%diag) call cpu_clock_begin(id_clock_pass) call pass_var(h, G%Domain) call cpu_clock_end(id_clock_pass) @@ -1067,7 +1073,7 @@ subroutine adjustments_dyn_legacy_split(u, v, h, dt, G, CS) if (CS%readjust_BT_trans) then call cpu_clock_begin(id_clock_continuity) call continuity(u, v, h, h_temp, uh_temp, vh_temp, dt, G, & - CS%continuity_CSp, CS%diag, OBC=CS%OBC) + CS%continuity_CSp, OBC=CS%OBC) call cpu_clock_end(id_clock_continuity) !$OMP parallel default(none) shared(is,ie,js,je,nz,CS,uh_temp,vh_temp) !$OMP do @@ -1407,7 +1413,7 @@ subroutine initialize_dyn_legacy_split(u, v, h, uh, vh, eta, Time, G, param_file if (.not. query_initialized(uh,"uh",restart_CS) .or. & .not. query_initialized(vh,"vh",restart_CS)) then h_tmp(:,:,:) = h(:,:,:) - call continuity(u, v, h, h_tmp, uh, vh, dt, G, CS%continuity_CSp, CS%diag, OBC=CS%OBC) + call continuity(u, v, h, h_tmp, uh, vh, dt, G, CS%continuity_CSp, OBC=CS%OBC) call cpu_clock_begin(id_clock_pass_init) call pass_var(h_tmp, G%Domain) call cpu_clock_end(id_clock_pass_init) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 39f9d86e3b..e23c97382f 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -605,7 +605,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & if (associated(CS%BT_cont) .or. CS%BT_use_layer_fluxes) then call cpu_clock_begin(id_clock_continuity) call continuity(u, v, h, hp, uh_in, vh_in, dt, G, & - CS%continuity_CSp, CS%diag, OBC=CS%OBC, visc_rem_u=CS%visc_rem_u, & + CS%continuity_CSp, OBC=CS%OBC, visc_rem_u=CS%visc_rem_u, & visc_rem_v=CS%visc_rem_v, BT_cont=CS%BT_cont) call cpu_clock_end(id_clock_continuity) if (BT_cont_BT_thick) then @@ -695,7 +695,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & ! uh = u_av * h ! hp = h + dt * div . uh call cpu_clock_begin(id_clock_continuity) - call continuity(up, vp, h, hp, uh, vh, dt, G, CS%continuity_CSp, CS%diag, & + call continuity(up, vp, h, hp, uh, vh, dt, G, CS%continuity_CSp, & CS%uhbt, CS%vhbt, CS%OBC, CS%visc_rem_u, CS%visc_rem_v, & u_av, v_av, BT_cont=CS%BT_cont) call cpu_clock_end(id_clock_continuity) @@ -901,9 +901,12 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & ! u_av and v_av adjusted so their mass transports match uhbt and vhbt. call cpu_clock_begin(id_clock_continuity) call continuity(u, v, h, h, uh, vh, dt, G, & - CS%continuity_CSp, CS%diag, CS%uhbt, CS%vhbt, CS%OBC, & + CS%continuity_CSp, CS%uhbt, CS%vhbt, CS%OBC, & CS%visc_rem_u, CS%visc_rem_v, u_av, v_av) call cpu_clock_end(id_clock_continuity) + ! Whenever thickness changes let the diag manager know, target grids + ! for vertical remapping may need to be regenerated. + call diag_update_target_grids(CS%diag) call cpu_clock_begin(id_clock_pass) call do_group_pass(CS%pass_h, G%Domain) call cpu_clock_end(id_clock_pass) @@ -1258,7 +1261,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, param_file, & if (.not. query_initialized(uh,"uh",restart_CS) .or. & .not. query_initialized(vh,"vh",restart_CS)) then h_tmp(:,:,:) = h(:,:,:) - call continuity(u, v, h, h_tmp, uh, vh, dt, G, CS%continuity_CSp, CS%diag, OBC=CS%OBC) + call continuity(u, v, h, h_tmp, uh, vh, dt, G, CS%continuity_CSp, OBC=CS%OBC) call cpu_clock_begin(id_clock_pass_init) call create_group_pass(pass_h_tmp, h_tmp, G%Domain) call do_group_pass(pass_h_tmp, G%Domain) diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index 9c51432819..ff911873f6 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -78,7 +78,7 @@ module MOM_dynamics_unsplit use MOM_diag_mediator, only : diag_mediator_init, enable_averaging use MOM_diag_mediator, only : disable_averaging, post_data, safe_alloc_ptr use MOM_diag_mediator, only : register_diag_field, register_static_field -use MOM_diag_mediator, only : set_diag_mediator_grid, diag_ctrl +use MOM_diag_mediator, only : set_diag_mediator_grid, diag_ctrl, diag_update_target_grids use MOM_domains, only : MOM_domains_init, pass_var, pass_vector use MOM_domains, only : pass_var_start, pass_var_complete use MOM_domains, only : pass_vector_start, pass_vector_complete @@ -268,7 +268,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, fluxes, & ! uh = u*h ! hp = h + dt/2 div . uh call cpu_clock_begin(id_clock_continuity) - call continuity(u, v, h, hp, uh, vh, dt*0.5, G, CS%continuity_CSp, CS%diag, & + call continuity(u, v, h, hp, uh, vh, dt*0.5, G, CS%continuity_CSp, & OBC=CS%OBC) call cpu_clock_end(id_clock_continuity) call cpu_clock_begin(id_clock_pass) @@ -378,7 +378,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, fluxes, & ! h_av = hp + dt/2 div . uh call cpu_clock_begin(id_clock_continuity) call continuity(up, vp, hp, h_av, uh, vh, & - (0.5*dt), G, CS%continuity_CSp, CS%diag, OBC=CS%OBC) + (0.5*dt), G, CS%continuity_CSp, OBC=CS%OBC) call cpu_clock_end(id_clock_continuity) call cpu_clock_begin(id_clock_pass) call pass_var(h_av, G%Domain) @@ -437,8 +437,11 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, fluxes, & ! h = hp + dt/2 div . uh call cpu_clock_begin(id_clock_continuity) call continuity(upp, vpp, hp, h, uh, vh, & - (dt*0.5), G, CS%continuity_CSp, CS%diag, OBC=CS%OBC) + (dt*0.5), G, CS%continuity_CSp, OBC=CS%OBC) call cpu_clock_end(id_clock_continuity) + ! Whenever thickness changes let the diag manager know, target grids + ! for vertical remapping may need to be regenerated. + call diag_update_target_grids(CS%diag) call cpu_clock_begin(id_clock_pass) call pass_var(h, G%Domain) call pass_vector(uh, vh, G%Domain) diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90 index 87a8b2052f..53ca73f4a0 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/src/core/MOM_dynamics_unsplit_RK2.F90 @@ -283,7 +283,7 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, ! This is a duplicate caclulation of the last continuity from the previous step ! and could/should be optimized out. -AJA call continuity(u_in, v_in, h_in, hp, uh, vh, dt_pred, G, CS%continuity_CSp, & - CS%diag, OBC=CS%OBC) + OBC=CS%OBC) call cpu_clock_end(id_clock_continuity) call cpu_clock_begin(id_clock_pass) call pass_var(hp, G%Domain) @@ -370,7 +370,7 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, ! h_av = h + dt div . uh call cpu_clock_begin(id_clock_continuity) call continuity(up, vp, h_in, hp, uh, vh, & - dt, G, CS%continuity_CSp, CS%diag, OBC=CS%OBC) + dt, G, CS%continuity_CSp, OBC=CS%OBC) call cpu_clock_end(id_clock_continuity) call cpu_clock_begin(id_clock_pass) call pass_var(hp, G%Domain) @@ -427,7 +427,7 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, ! h[n+1] = h[n] + dt div . uh call cpu_clock_begin(id_clock_continuity) call continuity(up, vp, h_in, h_in, uh, vh, & - dt, G, CS%continuity_CSp, CS%diag, OBC=CS%OBC) + dt, G, CS%continuity_CSp, OBC=CS%OBC) call cpu_clock_end(id_clock_continuity) call cpu_clock_begin(id_clock_pass) call pass_var(h_in, G%Domain) From 604681e8351371e5a137e92d6cb18a977119e858 Mon Sep 17 00:00:00 2001 From: Nic Hannah Date: Tue, 28 Jul 2015 10:46:37 -0400 Subject: [PATCH 17/18] Fix dummy argument intent. Does not compile on Gaea. #62 --- src/framework/MOM_diag_mediator.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index ee22320891..2386ec0de1 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -601,7 +601,7 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask) integer, intent(in) :: diag_field_id real, intent(in) :: field(:,:,:) - type(diag_ctrl), target, intent(inout) :: diag_cs + type(diag_ctrl), target, intent(in) :: diag_cs logical, optional, intent(in) :: is_static real, optional, intent(in) :: mask(:,:,:) From d0e99df664dedf1eefd597612331d39c5f917593 Mon Sep 17 00:00:00 2001 From: Nic Hannah Date: Tue, 28 Jul 2015 10:49:28 -0400 Subject: [PATCH 18/18] Don't build diagnostic Z remapping grids if grid definition param not provided. #62 --- src/framework/MOM_diag_mediator.F90 | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 2386ec0de1..dd0a901e8b 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -780,6 +780,13 @@ subroutine diag_update_target_grids(diag_cs) nz_src = size(diag_cs%h, 3) G => diag_cs%G + ! The interface positions for z remapping were not provided, so don't do + ! anything. If z remapping of diagnostics is requested then an error will + ! be triggered on post(). See param DIAG_REMAP_Z_GRID_DEF + if (.not. allocated(diag_cs%zi_remap)) then + return + endif + if (.not. diag_cs%remapping_initialized) then call assert(allocated(diag_cs%zi_remap), & 'update_diag_target_grids: Remapping axis not initialized')