mom_diag_remap Module Reference

Detailed Description

provides runtime remapping of diagnostics to z star, sigma and rho vertical coordinates.

The diag_remap_ctrl type represents a remapping of diagnostics to a particular vertical coordinate. The module is used by the diag mediator module in the following way:

1. diag_remap_init() is called to initialize a diag_remap_ctrl instance.
2. diag_remap_configure_axes() is called to read the configuration file and set up the vertical coordinate / axes definitions.
3. diag_remap_get_axes_info() returns information needed for the diag mediator to define new axes for the remapped diagnostics.
4. diag_remap_update() is called periodically (whenever h, T or S change) to either create or update the target remapping grids.
5. diag_remap_do_remap() is called from within a diag post() to do the remapping before the diagnostic is written out.

Data Types
type	diag_remap_ctrl
	Represents remapping of diagnostics to a particular vertical coordinate. More...

Functions/Subroutines
subroutine, public	diag_remap_init (remap_cs, coord_tuple, answers_2018)
	Initialize a diagnostic remapping type with the given vertical coordinate. More...
subroutine, public	diag_remap_end (remap_cs)
	De-init a diagnostic remapping type. Free allocated memory. More...
subroutine, public	diag_remap_diag_registration_closed (remap_cs)
	Inform that all diagnostics have been registered. If _set_active() has not been called on the remapping control structure will be disabled. This saves time in the case that a vertical coordinate was configured but no diagnostics which use the coordinate appeared in the diag_table. More...
subroutine, public	diag_remap_set_active (remap_cs)
	Indicate that this remapping type is actually used by the diag manager. If this is never called then the type will be disabled to save time. See further explanation with diag_remap_registration_closed. More...
subroutine, public	diag_remap_configure_axes (remap_cs, GV, US, param_file)
	Configure the vertical axes for a diagnostic remapping control structure. Reads a configuration parameters to determine coordinate generation. More...
subroutine, public	diag_remap_get_axes_info (remap_cs, nz, id_layer, id_interface)
	Get layer and interface axes ids for this coordinate Needed when defining axes groups. More...
logical function, public	diag_remap_axes_configured (remap_cs)
	Whether or not the axes for this vertical coordinated has been configured. Configuration is complete when diag_remap_configure_axes() has been successfully called. More...
subroutine, public	diag_remap_update (remap_cs, G, GV, US, h, T, S, eqn_of_state, h_target)
	Build/update target vertical grids for diagnostic remapping. More...
subroutine, public	diag_remap_do_remap (remap_cs, G, GV, h, staggered_in_x, staggered_in_y, mask, field, remapped_field)
	Remap diagnostic field to alternative vertical grid. More...
subroutine, public	diag_remap_calc_hmask (remap_cs, G, mask)
	Calculate masks for target grid. More...
subroutine, public	vertically_reintegrate_diag_field (remap_cs, G, h, h_target, staggered_in_x, staggered_in_y, mask, field, reintegrated_field)
	Vertically re-grid an already vertically-integrated diagnostic field to alternative vertical grid. More...
subroutine, public	vertically_interpolate_diag_field (remap_cs, G, h, staggered_in_x, staggered_in_y, mask, field, interpolated_field)
	Vertically interpolate diagnostic field to alternative vertical grid. More...
subroutine, public	horizontally_average_diag_field (G, GV, h, staggered_in_x, staggered_in_y, is_layer, is_extensive, field, averaged_field, averaged_mask)
	Horizontally average field. More...

Function/Subroutine Documentation

diag_remap_axes_configured()

logical function, public mom_diag_remap::diag_remap_axes_configured

(

type(diag_remap_ctrl), intent(in)

remap_cs

)

Whether or not the axes for this vertical coordinated has been configured. Configuration is complete when diag_remap_configure_axes() has been successfully called.

Parameters

[in]

remap_cs

Diagnostic coordinate control structure

Definition at line 262 of file MOM_diag_remap.F90.

  type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coordinate control structure
  logical :: diag_remap_axes_configured
 
  diag_remap_axes_configured = remap_cs%configured
 

diag_remap_calc_hmask()

subroutine, public mom_diag_remap::diag_remap_calc_hmask	(	type(diag_remap_ctrl), intent(in)	remap_cs,
		type(ocean_grid_type), intent(in)	G,
		real, dimension(:,:,:), intent(out)	mask
	)

Calculate masks for target grid.

Parameters

[in]	remap_cs	Diagnostic coodinate control structure
[in]	g	Ocean grid structure
[out]	mask	h-point mask for target grid [nondim]

Definition at line 446 of file MOM_diag_remap.F90.

  type(diag_remap_ctrl),  intent(in)  :: remap_cs !< Diagnostic coodinate control structure
  type(ocean_grid_type),  intent(in)  :: G    !< Ocean grid structure
  real, dimension(:,:,:), intent(out) :: mask !< h-point mask for target grid [nondim]
  ! Local variables
  real, dimension(remap_cs%nz) :: h_dest ! Destination thicknesses [H ~> m or kg m-2]
  integer :: i, j, k
  logical :: mask_vanished_layers
  real :: h_tot      ! Sum of all thicknesses [H ~> m or kg m-2]
  real :: h_err      ! An estimate of a negligible thickness [H ~> m or kg m-2]
 
  call assert(remap_cs%initialized, 'diag_remap_calc_hmask: remap_cs not initialized.')
 
  ! Only z*-like diagnostic coordinates should have a 3d mask
  mask_vanished_layers = (remap_cs%vertical_coord == coordinatemode('ZSTAR'))
  mask(:,:,:) = 0.
 
  do j=g%jsc-1, g%jec+1 ; do i=g%isc-1, g%iec+1
    if (g%mask2dT(i,j)>0.) then
      if (mask_vanished_layers) then
        h_dest(:) = remap_cs%h(i,j,:)
        h_tot = 0.
        h_err = 0.
        do k=1, remap_cs%nz
          h_tot = h_tot + h_dest(k)
          ! This is an overestimate of how thick a vanished layer might be, that
          ! appears due to round-off.
          h_err = h_err + epsilon(h_tot) * h_tot
          ! Mask out vanished layers
          if (h_dest(k)<=8.*h_err) then
            mask(i,j,k) = 0.
          else
            mask(i,j,k) = 1.
          endif
        enddo
      else ! all layers might contain data
        mask(i,j,:) = 1.
      endif
    endif
  enddo ; enddo
 

diag_remap_configure_axes()

subroutine, public mom_diag_remap::diag_remap_configure_axes	(	type(diag_remap_ctrl), intent(inout)	remap_cs,
		type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		type(param_file_type), intent(in)	param_file
	)

Configure the vertical axes for a diagnostic remapping control structure. Reads a configuration parameters to determine coordinate generation.

Parameters

[in,out]	remap_cs	Diag remap control structure
[in]	gv	ocean vertical grid structure
[in]	us	A dimensional unit scaling type
[in]	param_file	Parameter file structure

Definition at line 187 of file MOM_diag_remap.F90.

  type(diag_remap_ctrl),   intent(inout) :: remap_cs !< Diag remap control structure
  type(verticalGrid_type), intent(in)    :: GV !< ocean vertical grid structure
  type(unit_scale_type),   intent(in)    :: US !< A dimensional unit scaling type
  type(param_file_type),   intent(in)    :: param_file !< Parameter file structure
  ! Local variables
  integer :: nzi(4), nzl(4), k
  character(len=200) :: inputdir, string, filename, int_varname, layer_varname
  character(len=40)  :: mod  = "MOM_diag_remap" ! This module's name.
  character(len=8)   :: units, expected_units
  character(len=34)  :: longname, string2
 
  character(len=256) :: err_msg
  logical :: ierr
 
  real, allocatable, dimension(:) :: interfaces, layers
 
  call initialize_regridding(remap_cs%regrid_cs, gv, us, gv%max_depth, param_file, mod, &
           trim(remap_cs%vertical_coord_name), "DIAG_COORD", trim(remap_cs%diag_coord_name))
  call set_regrid_params(remap_cs%regrid_cs, min_thickness=0., integrate_downward_for_e=.false.)
 
  remap_cs%nz = get_regrid_size(remap_cs%regrid_cs)
 
  if (remap_cs%vertical_coord == coordinatemode('SIGMA')) then
    units = 'nondim'
    longname = 'Fraction'
  elseif (remap_cs%vertical_coord == coordinatemode('RHO')) then
    units = 'kg m-3'
    longname = 'Target Potential Density'
  else
    units = 'meters'
    longname = 'Depth'
  endif
 
  ! Make axes objects
  allocate(interfaces(remap_cs%nz+1))
  allocate(layers(remap_cs%nz))
 
  interfaces(:) = getcoordinateinterfaces(remap_cs%regrid_cs)
  layers(:) = 0.5 * ( interfaces(1:remap_cs%nz) + interfaces(2:remap_cs%nz+1) )
 
  remap_cs%interface_axes_id = diag_axis_init(lowercase(trim(remap_cs%diag_coord_name))//'_i', &
                                              interfaces, trim(units), 'z', &
                                              trim(longname)//' at interface', direction=-1)
  remap_cs%layer_axes_id = diag_axis_init(lowercase(trim(remap_cs%diag_coord_name))//'_l', &
                                          layers, trim(units), 'z', &
                                          trim(longname)//' at cell center', direction=-1, &
                                          edges=remap_cs%interface_axes_id)
 
  ! Axes have now been configured.
  remap_cs%configured = .true.
 
  deallocate(interfaces)
  deallocate(layers)
 

diag_remap_diag_registration_closed()

subroutine, public mom_diag_remap::diag_remap_diag_registration_closed

(

type(diag_remap_ctrl), intent(inout)

remap_cs

)

Inform that all diagnostics have been registered. If _set_active() has not been called on the remapping control structure will be disabled. This saves time in the case that a vertical coordinate was configured but no diagnostics which use the coordinate appeared in the diag_table.

Parameters

[in,out]

remap_cs

Diag remapping control structure

Definition at line 166 of file MOM_diag_remap.F90.

  type(diag_remap_ctrl), intent(inout) :: remap_cs !< Diag remapping control structure
 
  if (.not. remap_cs%used) then
    call diag_remap_end(remap_cs)
  endif
 

diag_remap_do_remap()

subroutine, public mom_diag_remap::diag_remap_do_remap	(	type(diag_remap_ctrl), intent(in)	remap_cs,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		real, dimension(:,:,:), intent(in)	h,
		logical, intent(in)	staggered_in_x,
		logical, intent(in)	staggered_in_y,
		real, dimension(:,:,:), pointer	mask,
		real, dimension(:,:,:), intent(in)	field,
		real, dimension(:,:,:), intent(inout)	remapped_field
	)

Remap diagnostic field to alternative vertical grid.

Parameters

[in]	remap_cs	Diagnostic coodinate control structure
[in]	g	Ocean grid structure
[in]	gv	ocean vertical grid structure
[in]	h	The current thicknesses [H ~> m or kg m-2]
[in]	staggered_in_x	True is the x-axis location is at u or q points
[in]	staggered_in_y	True is the y-axis location is at v or q points
	mask	A mask for the field [nondim]
[in]	field	The diagnostic field to be remapped [A]
[in,out]	remapped_field	Field remapped to new coordinate [A]

Definition at line 350 of file MOM_diag_remap.F90.

  type(diag_remap_ctrl),   intent(in) :: remap_cs !< Diagnostic coodinate control structure
  type(ocean_grid_type),   intent(in) :: G  !< Ocean grid structure
  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
  real, dimension(:,:,:),  intent(in) :: h  !< The current thicknesses [H ~> m or kg m-2]
  logical,                 intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points
  logical,                 intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points
  real, dimension(:,:,:),  pointer    :: mask !< A mask for the field [nondim]
  real, dimension(:,:,:),  intent(in) :: field(:,:,:) !< The diagnostic field to be remapped [A]
  real, dimension(:,:,:),  intent(inout) :: remapped_field !< Field remapped to new coordinate [A]
  ! Local variables
  real, dimension(remap_cs%nz) :: h_dest ! Destination thicknesses [H ~> m or kg m-2]
  real, dimension(size(h,3)) :: h_src    ! A column of source thicknesses [H ~> m or kg m-2]
  real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2]
  integer :: nz_src, nz_dest
  integer :: i, j, k                !< Grid index
  integer :: i1, j1                 !< 1-based index
  integer :: i_lo, i_hi, j_lo, j_hi !< (uv->h) interpolation indices
  integer :: shift                  !< Symmetric offset for 1-based indexing
 
  call assert(remap_cs%initialized, 'diag_remap_do_remap: remap_cs not initialized.')
  call assert(size(field, 3) == size(h, 3), &
              'diag_remap_do_remap: Remap field and thickness z-axes do not match.')
 
  if (.not.remap_cs%answers_2018) then
    h_neglect = gv%H_subroundoff ; h_neglect_edge = gv%H_subroundoff
  elseif (gv%Boussinesq) then
    h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
  else
    h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
  endif
 
  nz_src = size(field,3)
  nz_dest = remap_cs%nz
  remapped_field(:,:,:) = 0.
 
  ! Symmetric grid offset under 1-based indexing; see header for details.
  shift = 0; if (g%symmetric) shift = 1
 
  if (staggered_in_x .and. .not. staggered_in_y) then
    ! U-points
    do j=g%jsc, g%jec
      do i=g%iscB, g%iecB
        i1 = i - g%isdB + 1
        i_lo = i1 - shift; i_hi = i_lo + 1
        if (associated(mask)) then
          if (mask(i,j,1) == 0.) cycle
        endif
        h_src(:) = 0.5 * (h(i_lo,j,:) + h(i_hi,j,:))
        h_dest(:) = 0.5 * (remap_cs%h(i_lo,j,:) + remap_cs%h(i_hi,j,:))
        call remapping_core_h(remap_cs%remap_cs, &
                              nz_src, h_src(:), field(i1,j,:), &
                              nz_dest, h_dest(:), remapped_field(i1,j,:), &
                              h_neglect, h_neglect_edge)
      enddo
    enddo
  elseif (staggered_in_y .and. .not. staggered_in_x) then
    ! V-points
    do j=g%jscB, g%jecB
      j1 = j - g%jsdB + 1
      j_lo = j1 - shift; j_hi = j_lo + 1
      do i=g%isc, g%iec
        if (associated(mask)) then
          if (mask(i,j,1) == 0.) cycle
        endif
        h_src(:) = 0.5 * (h(i,j_lo,:) + h(i,j_hi,:))
        h_dest(:) = 0.5 * (remap_cs%h(i,j_lo,:) + remap_cs%h(i,j_hi,:))
        call remapping_core_h(remap_cs%remap_cs, &
                              nz_src, h_src(:), field(i,j1,:), &
                              nz_dest, h_dest(:), remapped_field(i,j1,:), &
                              h_neglect, h_neglect_edge)
      enddo
    enddo
  elseif ((.not. staggered_in_x) .and. (.not. staggered_in_y)) then
    ! H-points
    do j=g%jsc, g%jec
      do i=g%isc, g%iec
        if (associated(mask)) then
          if (mask(i,j,1) == 0.) cycle
        endif
        h_src(:) = h(i,j,:)
        h_dest(:) = remap_cs%h(i,j,:)
        call remapping_core_h(remap_cs%remap_cs, &
                              nz_src, h_src(:), field(i,j,:), &
                              nz_dest, h_dest(:), remapped_field(i,j,:), &
                              h_neglect, h_neglect_edge)
      enddo
    enddo
  else
    call assert(.false., 'diag_remap_do_remap: Unsupported axis combination')
  endif
 

diag_remap_end()

subroutine, public mom_diag_remap::diag_remap_end

(

type(diag_remap_ctrl), intent(inout)

remap_cs

)

De-init a diagnostic remapping type. Free allocated memory.

Parameters

[in,out]

remap_cs

Diag remapping control structure

Definition at line 150 of file MOM_diag_remap.F90.

  type(diag_remap_ctrl), intent(inout) :: remap_cs !< Diag remapping control structure
 
  if (allocated(remap_cs%h)) deallocate(remap_cs%h)
  remap_cs%configured = .false.
  remap_cs%initialized = .false.
  remap_cs%used = .false.
  remap_cs%nz = 0
 

diag_remap_get_axes_info()

subroutine, public mom_diag_remap::diag_remap_get_axes_info	(	type(diag_remap_ctrl), intent(in)	remap_cs,
		integer, intent(out)	nz,
		integer, intent(out)	id_layer,
		integer, intent(out)	id_interface
	)

Get layer and interface axes ids for this coordinate Needed when defining axes groups.

Parameters

[in]	remap_cs	Diagnostic coordinate control structure
[out]	nz	Number of vertical levels for the coordinate
[out]	id_layer	1D-axes id for layer points
[out]	id_interface	1D-axes id for interface points

Definition at line 246 of file MOM_diag_remap.F90.

  type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coordinate control structure
  integer, intent(out) :: nz !< Number of vertical levels for the coordinate
  integer, intent(out) :: id_layer !< 1D-axes id for layer points
  integer, intent(out) :: id_interface !< 1D-axes id for interface points
 
  nz = remap_cs%nz
  id_layer = remap_cs%layer_axes_id
  id_interface = remap_cs%interface_axes_id
 

diag_remap_init()

subroutine, public mom_diag_remap::diag_remap_init	(	type(diag_remap_ctrl), intent(inout)	remap_cs,
		character(len=*), intent(in)	coord_tuple,
		logical, intent(in)	answers_2018
	)

Initialize a diagnostic remapping type with the given vertical coordinate.

Parameters

[in,out]	remap_cs	Diag remapping control structure
[in]	coord_tuple	A string in form of MODULE_SUFFIX PARAMETER_SUFFIX COORDINATE_NAME
[in]	answers_2018	If true, use the order of arithmetic and expressions for remapping that recover the answers from the end of 2018. Otherwise, use more robust forms of the same expressions.

Definition at line 128 of file MOM_diag_remap.F90.

  type(diag_remap_ctrl), intent(inout) :: remap_cs !< Diag remapping control structure
  character(len=*),      intent(in)    :: coord_tuple !< A string in form of
                                                      !! MODULE_SUFFIX PARAMETER_SUFFIX COORDINATE_NAME
  logical,               intent(in)    :: answers_2018 !< If true, use the order of arithmetic and expressions
                                                      !! for remapping that recover the answers from the end of 2018.
                                                      !! Otherwise, use more robust forms of the same expressions.
 
  remap_cs%diag_module_suffix = trim(extractword(coord_tuple, 1))
  remap_cs%diag_coord_name = trim(extractword(coord_tuple, 2))
  remap_cs%vertical_coord_name = trim(extractword(coord_tuple, 3))
  remap_cs%vertical_coord = coordinatemode(remap_cs%vertical_coord_name)
  remap_cs%configured = .false.
  remap_cs%initialized = .false.
  remap_cs%used = .false.
  remap_cs%answers_2018 = answers_2018
  remap_cs%nz = 0
 

diag_remap_set_active()

subroutine, public mom_diag_remap::diag_remap_set_active

(

type(diag_remap_ctrl), intent(inout)

remap_cs

)

Indicate that this remapping type is actually used by the diag manager. If this is never called then the type will be disabled to save time. See further explanation with diag_remap_registration_closed.

Parameters

[in,out]

remap_cs

Diag remapping control structure

Definition at line 178 of file MOM_diag_remap.F90.

  type(diag_remap_ctrl), intent(inout) :: remap_cs !< Diag remapping control structure
 
  remap_cs%used = .true.
 

diag_remap_update()

subroutine, public mom_diag_remap::diag_remap_update	(	type(diag_remap_ctrl), intent(inout)	remap_cs,
		type(ocean_grid_type), pointer	G,
		type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		real, dimension(:,:,:), intent(in)	h,
		real, dimension(:,:,:), intent(in)	T,
		real, dimension(:,:,:), intent(in)	S,
		type(eos_type), pointer	eqn_of_state,
		real, dimension(:,:,:), intent(inout)	h_target
	)

Build/update target vertical grids for diagnostic remapping.

Note

The target grids need to be updated whenever sea surface height or layer thicknesses changes. In the case of density-based coordinates then technically we should also regenerate the target grid whenever T/S change.

Parameters

[in,out]	remap_cs	Diagnostic coordinate control structure
	g	The ocean's grid type
[in]	gv	ocean vertical grid structure
[in]	us	A dimensional unit scaling type
[in]	h	New thickness [H ~> m or kg m-2]
[in]	t	New temperatures [degC]
[in]	s	New salinities [ppt]
	eqn_of_state	A pointer to the equation of state
[in,out]	h_target	The new diagnostic thicknesses [H ~> m or kg m-2]

Definition at line 275 of file MOM_diag_remap.F90.

  type(diag_remap_ctrl),   intent(inout) :: remap_cs !< Diagnostic coordinate control structure
  type(ocean_grid_type),   pointer    :: G  !< The ocean's grid type
  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
  type(unit_scale_type),   intent(in) :: US !< A dimensional unit scaling type
  real, dimension(:,:,:),  intent(in) :: h  !< New thickness [H ~> m or kg m-2]
  real, dimension(:,:,:),  intent(in) :: T  !< New temperatures [degC]
  real, dimension(:,:,:),  intent(in) :: S  !< New salinities [ppt]
  type(EOS_type),          pointer    :: eqn_of_state !< A pointer to the equation of state
  real, dimension(:,:,:),  intent(inout) :: h_target  !< The new diagnostic thicknesses [H ~> m or kg m-2]
 
  ! Local variables
  real, dimension(remap_cs%nz + 1) :: zInterfaces ! Interface positions [H ~> m or kg m-2]
  real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2]
  integer :: i, j, k, nz
 
  ! Note that coordinateMode('LAYER') is never 'configured' so will
  ! always return here.
  if (.not. remap_cs%configured) then
    return
  endif
 
  if (.not.remap_cs%answers_2018) then
    h_neglect = gv%H_subroundoff ; h_neglect_edge = gv%H_subroundoff
  elseif (gv%Boussinesq) then
    h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
  else
    h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
  endif
  nz = remap_cs%nz
 
  if (.not. remap_cs%initialized) then
    ! Initialize remapping and regridding on the first call
    call initialize_remapping(remap_cs%remap_cs, 'PPM_IH4', boundary_extrapolation=.false., &
                              answers_2018=remap_cs%answers_2018)
    remap_cs%initialized = .true.
  endif
 
  ! Calculate remapping thicknesses for different target grids based on
  ! nominal/target interface locations. This happens for every call on the
  ! assumption that h, T, S has changed.
  do j=g%jsc-1, g%jec+1 ; do i=g%isc-1, g%iec+1
    if (g%mask2dT(i,j)==0.) then
      h_target(i,j,:) = 0.
      cycle
    endif
 
    if (remap_cs%vertical_coord == coordinatemode('ZSTAR')) then
      call build_zstar_column(get_zlike_cs(remap_cs%regrid_cs), &
                              gv%Z_to_H*g%bathyT(i,j), sum(h(i,j,:)), &
                              zinterfaces, zscale=gv%Z_to_H)
    elseif (remap_cs%vertical_coord == coordinatemode('SIGMA')) then
      call build_sigma_column(get_sigma_cs(remap_cs%regrid_cs), &
                              gv%Z_to_H*g%bathyT(i,j), sum(h(i,j,:)), zinterfaces)
    elseif (remap_cs%vertical_coord == coordinatemode('RHO')) then
      call build_rho_column(get_rho_cs(remap_cs%regrid_cs), g%ke, &
                            gv%Z_to_H*g%bathyT(i,j), h(i,j,:), t(i,j,:), s(i,j,:), &
                            eqn_of_state, zinterfaces, h_neglect, h_neglect_edge)
    elseif (remap_cs%vertical_coord == coordinatemode('SLIGHT')) then
!     call build_slight_column(remap_cs%regrid_cs,remap_cs%remap_cs, nz, &
!                           GV%Z_to_H*G%bathyT(i,j), sum(h(i,j,:)), zInterfaces)
      call mom_error(fatal,"diag_remap_update: SLIGHT coordinate not coded for diagnostics yet!")
    elseif (remap_cs%vertical_coord == coordinatemode('HYCOM1')) then
!     call build_hycom1_column(remap_cs%regrid_cs, nz, &
!                           GV%Z_to_H*G%bathyT(i,j), sum(h(i,j,:)), zInterfaces)
      call mom_error(fatal,"diag_remap_update: HYCOM1 coordinate not coded for diagnostics yet!")
    endif
    do k = 1,nz
      h_target(i,j,k) = zinterfaces(k) - zinterfaces(k+1)
    enddo
  enddo ; enddo
 

horizontally_average_diag_field()

subroutine, public mom_diag_remap::horizontally_average_diag_field	(	type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		real, dimension(:,:,:), intent(in)	h,
		logical, intent(in)	staggered_in_x,
		logical, intent(in)	staggered_in_y,
		logical, intent(in)	is_layer,
		logical, intent(in)	is_extensive,
		real, dimension(:,:,:), intent(in)	field,
		real, dimension(:), intent(inout)	averaged_field,
		logical, dimension(:), intent(inout)	averaged_mask
	)

Horizontally average field.

Parameters

[in]	g	Ocean grid structure
[in]	gv	The ocean vertical grid structure
[in]	h	The current thicknesses [H ~> m or kg m-2]
[in]	staggered_in_x	True if the x-axis location is at u or q points
[in]	staggered_in_y	True if the y-axis location is at v or q points
[in]	is_layer	True if the z-axis location is at h points
[in]	is_extensive	True if the z-direction is spatially integrated (over layers)
[in]	field	The diagnostic field to be remapped [A]
[in,out]	averaged_field	Field argument horizontally averaged [A]
[in,out]	averaged_mask	Mask for horizontally averaged field [nondim]

Definition at line 652 of file MOM_diag_remap.F90.

  type(ocean_grid_type),  intent(in) :: G !< Ocean grid structure
  type(verticalGrid_type), intent(in) :: GV !< The ocean vertical grid structure
  real, dimension(:,:,:), intent(in) :: h !< The current thicknesses [H ~> m or kg m-2]
  logical,                intent(in) :: staggered_in_x !< True if the x-axis location is at u or q points
  logical,                intent(in) :: staggered_in_y !< True if the y-axis location is at v or q points
  logical,                intent(in) :: is_layer !< True if the z-axis location is at h points
  logical,                intent(in) :: is_extensive !< True if the z-direction is spatially integrated (over layers)
  real, dimension(:,:,:), intent(in) :: field !<  The diagnostic field to be remapped [A]
  real, dimension(:),  intent(inout) :: averaged_field !< Field argument horizontally averaged [A]
  logical, dimension(:), intent(inout) :: averaged_mask  !< Mask for horizontally averaged field [nondim]
 
  ! Local variables
  real, dimension(G%isc:G%iec, G%jsc:G%jec, size(field,3)) :: volume, stuff
  real, dimension(size(field, 3)) :: vol_sum, stuff_sum ! nz+1 is needed for interface averages
  type(EFP_type), dimension(2*size(field,3)) :: sums_EFP ! Sums of volume or stuff by layer
  real :: height  ! An average thickness attributed to an velocity point [H ~> m or kg m-2]
  integer :: i, j, k, nz
  integer :: i1, j1                 !< 1-based index
 
  nz = size(field, 3)
 
  ! TODO: These averages could potentially be modified to use the function in
  !       the MOM_spatial_means module.
  ! NOTE: Reproducible sums must be computed in the original MKS units
 
  if (staggered_in_x .and. .not. staggered_in_y) then
    if (is_layer) then
      ! U-points
      do k=1,nz
        vol_sum(k) = 0.
        stuff_sum(k) = 0.
        if (is_extensive) then
          do j=g%jsc, g%jec ; do i=g%isc, g%iec
            i1 = i - g%isdB + 1
            volume(i,j,k) = (g%US%L_to_m**2 * g%areaCu(i,j)) * g%mask2dCu(i,j)
            stuff(i,j,k) = volume(i,j,k) * field(i1,j,k)
          enddo ; enddo
        else ! Intensive
          do j=g%jsc, g%jec ; do i=g%isc, g%iec
            i1 = i - g%isdB + 1
            height = 0.5 * (h(i,j,k) + h(i+1,j,k))
            volume(i,j,k) = (g%US%L_to_m**2 * g%areaCu(i,j)) &
                * (gv%H_to_m * height) * g%mask2dCu(i,j)
            stuff(i,j,k) = volume(i,j,k) * field(i1,j,k)
          enddo ; enddo
        endif
      enddo
    else ! Interface
      do k=1,nz
        do j=g%jsc, g%jec ; do i=g%isc, g%iec
          i1 = i - g%isdB + 1
          volume(i,j,k) = (g%US%L_to_m**2 * g%areaCu(i,j)) * g%mask2dCu(i,j)
          stuff(i,j,k) = volume(i,j,k) * field(i1,j,k)
        enddo ; enddo
      enddo
    endif
  elseif (staggered_in_y .and. .not. staggered_in_x) then
    if (is_layer) then
      ! V-points
      do k=1,nz
        if (is_extensive) then
          do j=g%jsc, g%jec ; do i=g%isc, g%iec
            j1 = j - g%jsdB + 1
            volume(i,j,k) = (g%US%L_to_m**2 * g%areaCv(i,j)) * g%mask2dCv(i,j)
            stuff(i,j,k) = volume(i,j,k) * field(i,j1,k)
          enddo ; enddo
        else ! Intensive
          do j=g%jsc, g%jec ; do i=g%isc, g%iec
            j1 = j - g%jsdB + 1
            height = 0.5 * (h(i,j,k) + h(i,j+1,k))
            volume(i,j,k) = (g%US%L_to_m**2 * g%areaCv(i,j)) &
                * (gv%H_to_m * height) * g%mask2dCv(i,j)
            stuff(i,j,k) = volume(i,j,k) * field(i,j1,k)
          enddo ; enddo
        endif
      enddo
    else ! Interface
      do k=1,nz
        do j=g%jsc, g%jec ; do i=g%isc, g%iec
          j1 = j - g%jsdB + 1
          volume(i,j,k) = (g%US%L_to_m**2 * g%areaCv(i,j)) * g%mask2dCv(i,j)
          stuff(i,j,k) = volume(i,j,k) * field(i,j1,k)
        enddo ; enddo
      enddo
    endif
  elseif ((.not. staggered_in_x) .and. (.not. staggered_in_y)) then
    if (is_layer) then
      ! H-points
      do k=1,nz
        if (is_extensive) then
          do j=g%jsc, g%jec ; do i=g%isc, g%iec
            if (h(i,j,k) > 0.) then
              volume(i,j,k) = (g%US%L_to_m**2 * g%areaT(i,j)) * g%mask2dT(i,j)
              stuff(i,j,k) = volume(i,j,k) * field(i,j,k)
            else
              volume(i,j,k) = 0.
              stuff(i,j,k) = 0.
            endif
          enddo ; enddo
        else ! Intensive
          do j=g%jsc, g%jec ; do i=g%isc, g%iec
            volume(i,j,k) = (g%US%L_to_m**2 * g%areaT(i,j)) &
                * (gv%H_to_m * h(i,j,k)) * g%mask2dT(i,j)
            stuff(i,j,k) = volume(i,j,k) * field(i,j,k)
          enddo ; enddo
        endif
      enddo
    else ! Interface
      do k=1,nz
        do j=g%jsc, g%jec ; do i=g%isc, g%iec
          volume(i,j,k) = (g%US%L_to_m**2 * g%areaT(i,j)) * g%mask2dT(i,j)
          stuff(i,j,k) = volume(i,j,k) * field(i,j,k)
        enddo ; enddo
      enddo
    endif
  else
    call assert(.false., 'horizontally_average_diag_field: Q point averaging is not coded yet.')
  endif
 
  ! Packing the sums into a single array with a single call to sum across PEs saves reduces
  ! the costs of communication.
  do k=1,nz
    sums_efp(2*k-1) = reproducing_sum_efp(volume(:,:,k), only_on_pe=.true.)
    sums_efp(2*k)   = reproducing_sum_efp(stuff(:,:,k), only_on_pe=.true.)
  enddo
  call efp_sum_across_pes(sums_efp, 2*nz)
  do k=1,nz
    vol_sum(k) = efp_to_real(sums_efp(2*k-1))
    stuff_sum(k) = efp_to_real(sums_efp(2*k))
  enddo
 
  averaged_mask(:) = .true.
  do k=1,nz
    if (vol_sum(k) > 0.) then
      averaged_field(k) = stuff_sum(k) / vol_sum(k)
    else
      averaged_field(k) = 0.
      averaged_mask(k) = .false.
    endif
  enddo
 

vertically_interpolate_diag_field()

subroutine, public mom_diag_remap::vertically_interpolate_diag_field	(	type(diag_remap_ctrl), intent(in)	remap_cs,
		type(ocean_grid_type), intent(in)	G,
		real, dimension(:,:,:), intent(in)	h,
		logical, intent(in)	staggered_in_x,
		logical, intent(in)	staggered_in_y,
		real, dimension(:,:,:), pointer	mask,
		real, dimension(:,:,:), intent(in)	field,
		real, dimension(:,:,:), intent(inout)	interpolated_field
	)

Vertically interpolate diagnostic field to alternative vertical grid.

Parameters

[in]	remap_cs	Diagnostic coodinate control structure
[in]	g	Ocean grid structure
[in]	h	The current thicknesses [H ~> m or kg m-2]
[in]	staggered_in_x	True is the x-axis location is at u or q points
[in]	staggered_in_y	True is the y-axis location is at v or q points
	mask	A mask for the field [nondim]
[in]	field	The diagnostic field to be remapped [A]
[in,out]	interpolated_field	Field argument remapped to alternative coordinate [A]

Definition at line 571 of file MOM_diag_remap.F90.

  type(diag_remap_ctrl),  intent(in) :: remap_cs !< Diagnostic coodinate control structure
  type(ocean_grid_type),  intent(in) :: G   !< Ocean grid structure
  real, dimension(:,:,:), intent(in) :: h   !< The current thicknesses [H ~> m or kg m-2]
  logical,                intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points
  logical,                intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points
  real, dimension(:,:,:), pointer    :: mask !< A mask for the field [nondim]
  real, dimension(:,:,:), intent(in) :: field !<  The diagnostic field to be remapped [A]
  real, dimension(:,:,:), intent(inout) :: interpolated_field !< Field argument remapped to alternative coordinate [A]
  ! Local variables
  real, dimension(remap_cs%nz) :: h_dest ! Destination thicknesses [H ~> m or kg m-2]
  real, dimension(size(h,3)) :: h_src    ! A column of source thicknesses [H ~> m or kg m-2]
  integer :: nz_src, nz_dest
  integer :: i, j, k                !< Grid index
  integer :: i1, j1                 !< 1-based index
  integer :: i_lo, i_hi, j_lo, j_hi !< (uv->h) interpolation indices
  integer :: shift                  !< Symmetric offset for 1-based indexing
 
  call assert(remap_cs%initialized, 'vertically_interpolate_diag_field: remap_cs not initialized.')
  call assert(size(field, 3) == size(h, 3)+1, &
              'vertically_interpolate_diag_field: Remap field and thickness z-axes do not match.')
 
  interpolated_field(:,:,:) = 0.
 
  nz_src = size(h,3)
  nz_dest = remap_cs%nz
 
  ! Symmetric grid offset under 1-based indexing; see header for details.
  shift = 0; if (g%symmetric) shift = 1
 
  if (staggered_in_x .and. .not. staggered_in_y) then
    ! U-points
    do j=g%jsc, g%jec
      do i=g%iscB, g%iecB
        i1 = i - g%isdB + 1
        i_lo = i1 - shift; i_hi = i_lo + 1
        if (associated(mask)) then
          if (mask(i,j,1) == 0.) cycle
        endif
        h_src(:) = 0.5 * (h(i_lo,j,:) + h(i_hi,j,:))
        h_dest(:) = 0.5 * (remap_cs%h(i_lo,j,:) + remap_cs%h(i_hi,j,:))
        call interpolate_column(nz_src, h_src, field(i1,j,:), &
                                nz_dest, h_dest, 0., interpolated_field(i1,j,:))
      enddo
    enddo
  elseif (staggered_in_y .and. .not. staggered_in_x) then
    ! V-points
    do j=g%jscB, g%jecB
      j1 = j - g%jsdB + 1
      j_lo = j1 - shift; j_hi = j_lo + 1
      do i=g%isc, g%iec
        if (associated(mask)) then
          if (mask(i,j,1) == 0.) cycle
        endif
        h_src(:) = 0.5 * (h(i,j_lo,:) + h(i,j_hi,:))
        h_dest(:) = 0.5 * (remap_cs%h(i,j_lo,:) + remap_cs%h(i,j_hi,:))
        call interpolate_column(nz_src, h_src, field(i,j1,:), &
                                nz_dest, h_dest, 0., interpolated_field(i,j1,:))
      enddo
    enddo
  elseif ((.not. staggered_in_x) .and. (.not. staggered_in_y)) then
    ! H-points
    do j=g%jsc, g%jec
      do i=g%isc, g%iec
        if (associated(mask)) then
          if (mask(i,j,1) == 0.) cycle
        endif
        h_src(:) = h(i,j,:)
        h_dest(:) = remap_cs%h(i,j,:)
        call interpolate_column(nz_src, h_src, field(i,j,:), &
                                nz_dest, h_dest, 0., interpolated_field(i,j,:))
      enddo
    enddo
  else
    call assert(.false., 'vertically_interpolate_diag_field: Q point remapping is not coded yet.')
  endif
 

vertically_reintegrate_diag_field()

subroutine, public mom_diag_remap::vertically_reintegrate_diag_field	(	type(diag_remap_ctrl), intent(in)	remap_cs,
		type(ocean_grid_type), intent(in)	G,
		real, dimension(:,:,:), intent(in)	h,
		real, dimension(:,:,:), intent(in)	h_target,
		logical, intent(in)	staggered_in_x,
		logical, intent(in)	staggered_in_y,
		real, dimension(:,:,:), pointer	mask,
		real, dimension(:,:,:), intent(in)	field,
		real, dimension(:,:,:), intent(inout)	reintegrated_field
	)

Vertically re-grid an already vertically-integrated diagnostic field to alternative vertical grid.

Parameters

[in]	remap_cs	Diagnostic coodinate control structure
[in]	g	Ocean grid structure
[in]	h	The thicknesses of the source grid [H ~> m or kg m-2]
[in]	h_target	The thicknesses of the target grid [H ~> m or kg m-2]
[in]	staggered_in_x	True is the x-axis location is at u or q points
[in]	staggered_in_y	True is the y-axis location is at v or q points
	mask	A mask for the field [nondim]
[in]	field	The diagnostic field to be remapped [A]
[in,out]	reintegrated_field	Field argument remapped to alternative coordinate [A]

Definition at line 490 of file MOM_diag_remap.F90.

  type(diag_remap_ctrl),  intent(in) :: remap_cs !< Diagnostic coodinate control structure
  type(ocean_grid_type),  intent(in) :: G        !< Ocean grid structure
  real, dimension(:,:,:), intent(in) :: h        !< The thicknesses of the source grid [H ~> m or kg m-2]
  real, dimension(:,:,:), intent(in) :: h_target !< The thicknesses of the target grid [H ~> m or kg m-2]
  logical,                intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points
  logical,                intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points
  real, dimension(:,:,:), pointer    :: mask     !< A mask for the field [nondim]
  real, dimension(:,:,:), intent(in) :: field    !<  The diagnostic field to be remapped [A]
  real, dimension(:,:,:), intent(inout) :: reintegrated_field !< Field argument remapped to alternative coordinate [A]
  ! Local variables
  real, dimension(remap_cs%nz) :: h_dest ! Destination thicknesses [H ~> m or kg m-2]
  real, dimension(size(h,3)) :: h_src    ! A column of source thicknesses [H ~> m or kg m-2]
  integer :: nz_src, nz_dest
  integer :: i, j, k                !< Grid index
  integer :: i1, j1                 !< 1-based index
  integer :: i_lo, i_hi, j_lo, j_hi !< (uv->h) interpolation indices
  integer :: shift                  !< Symmetric offset for 1-based indexing
 
  call assert(remap_cs%initialized, 'vertically_reintegrate_diag_field: remap_cs not initialized.')
  call assert(size(field, 3) == size(h, 3), &
              'vertically_reintegrate_diag_field: Remap field and thickness z-axes do not match.')
 
  nz_src = size(field,3)
  nz_dest = remap_cs%nz
  reintegrated_field(:,:,:) = 0.
 
  ! Symmetric grid offset under 1-based indexing; see header for details.
  shift = 0; if (g%symmetric) shift = 1
 
  if (staggered_in_x .and. .not. staggered_in_y) then
    ! U-points
    do j=g%jsc, g%jec
      do i=g%iscB, g%iecB
        i1 = i - g%isdB + 1
        i_lo = i1 - shift; i_hi = i_lo + 1
        if (associated(mask)) then
          if (mask(i,j,1) == 0.) cycle
        endif
        h_src(:) = 0.5 * (h(i_lo,j,:) + h(i_hi,j,:))
        h_dest(:) = 0.5 * (h_target(i_lo,j,:) + h_target(i_hi,j,:))
        call reintegrate_column(nz_src, h_src, field(i1,j,:), &
                                nz_dest, h_dest, 0., reintegrated_field(i1,j,:))
      enddo
    enddo
  elseif (staggered_in_y .and. .not. staggered_in_x) then
    ! V-points
    do j=g%jscB, g%jecB
      j1 = j - g%jsdB + 1
      j_lo = j1 - shift; j_hi = j_lo + 1
      do i=g%isc, g%iec
        if (associated(mask)) then
          if (mask(i,j,1) == 0.) cycle
        endif
        h_src(:) = 0.5 * (h(i,j_lo,:) + h(i,j_hi,:))
        h_dest(:) = 0.5 * (h_target(i,j_lo,:) + h_target(i,j_hi,:))
        call reintegrate_column(nz_src, h_src, field(i,j1,:), &
                                nz_dest, h_dest, 0., reintegrated_field(i,j1,:))
      enddo
    enddo
  elseif ((.not. staggered_in_x) .and. (.not. staggered_in_y)) then
    ! H-points
    do j=g%jsc, g%jec
      do i=g%isc, g%iec
        if (associated(mask)) then
          if (mask(i,j,1) == 0.) cycle
        endif
        h_src(:) = h(i,j,:)
        h_dest(:) = h_target(i,j,:)
        call reintegrate_column(nz_src, h_src, field(i,j,:), &
                                nz_dest, h_dest, 0., reintegrated_field(i,j,:))
      enddo
    enddo
  else
    call assert(.false., 'vertically_reintegrate_diag_field: Q point remapping is not coded yet.')
  endif