Skip to content

Commit

Permalink
+Revise interfaces to horiz_interp_and_extract
Browse files Browse the repository at this point in the history 
  Revised the interface to the two horiz_interp_and_extract to eliminate the
unused reentrant_x and tripolar_n arguments and to rename the conversion
argument to scale and make it the last mandatory argument in anticipation that
we might decide to make it optional, similarly to other routines like
MOM_read_data.  The renaming is because we use 'conversion' to indicate how the
internal representation is to be rescaled for output, most prominently in
register_diag_field, whereas everywhere else where we are rescaling input to the
model's internal units, we use a scale argument.  This makes the convention
self-consistent across the MOM6 code, and should avoid some confusion.  The
calls to these routines were updated in 8 places in 4 modules.  In 8 places the
get_param calls for TRIPOLAR_N or REENTRANT_X that are no longer needed were
eliminated as were the associated internal variables.  This commit also includes
some additions to comments near the changes that were directly tied to this
commit.  All answers are bitwise identical, but there are changes to a publicly
visible interface; code that tries to use the old interface will not compile.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Dec 20, 2022
1 parent 022531c commit 0e2490d
Show file tree
Hide file tree
Showing 5 changed files with 97 additions and 104 deletions.
59 changes: 30 additions & 29 deletions src/framework/MOM_horizontal_regridding.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -250,14 +250,13 @@ subroutine fill_miss_2d(aout, good, fill, prev, G, acrit, num_pass, relc, debug,
end subroutine fill_miss_2d

!> Extrapolate and interpolate from a file record
subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, recnum, G, tr_z, mask_z, &
z_in, z_edges_in, missing_value, reentrant_x, tripolar_n, &
subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, recnum, G, tr_z, mask_z, &
z_in, z_edges_in, missing_value, scale, &
homogenize, m_to_Z, answers_2018, ongrid, tr_iter_tol, answer_date)

character(len=*), intent(in) :: filename !< Path to file containing tracer to be
!! interpolated.
character(len=*), intent(in) :: varnam !< Name of tracer in file.
real, intent(in) :: conversion !< Conversion factor for tracer [CU conc-1 ~> 1]
integer, intent(in) :: recnum !< Record number of tracer to be read.
type(ocean_grid_type), intent(inout) :: G !< Grid object
real, allocatable, dimension(:,:,:), intent(out) :: tr_z
Expand All @@ -271,10 +270,10 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion,
real, allocatable, dimension(:), intent(out) :: z_edges_in
!< Cell grid edge values for input data [Z ~> m]
real, intent(out) :: missing_value !< The missing value in the returned array, scaled
!! with conversion to avoid accidentally having valid
!! values match missing values [CU ~> conc]
logical, intent(in) :: reentrant_x !< If true, this grid is reentrant in the x-direction
logical, intent(in) :: tripolar_n !< If true, this is a northern tripolar grid
!! to avoid accidentally having valid values match
!! missing values [CU ~> conc]
real, intent(in) :: scale !< Scaling factor for tracer into the internal
!! units of the model [CU conc-1 ~> 1]
logical, optional, intent(in) :: homogenize !< If present and true, horizontally homogenize data
!! to produce perfectly "flat" initial conditions
real, optional, intent(in) :: m_to_Z !< A conversion factor from meters to the units
Expand Down Expand Up @@ -311,7 +310,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion,
real, dimension(:), allocatable :: lon_in ! The longitudes in the input file [degreesE] then [radians]
real, dimension(:), allocatable :: lat_in ! The latitudes in the input file [degreesN] then [radians]
real, dimension(:), allocatable :: lat_inp ! The input file latitudes expanded to the pole [degreesN] then [radians]
real :: max_lat ! The maximum latitude on the input grid [degreeN]
real :: max_lat ! The maximum latitude on the input grid [degreesN]
real :: pole ! The sum of tracer values at the pole [conc]
real :: max_depth ! The maximum depth of the ocean [Z ~> m]
real :: npole ! The number of points contributing to the pole value [nondim]
Expand All @@ -330,7 +329,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion,
integer :: isd, ied, jsd, jed ! data domain indices
integer :: id_clock_read
logical :: debug=.false.
real :: I_scale ! The inverse of the conversion factor for diagnostic output [conc CU-1 ~> 1]
real :: I_scale ! The inverse of the scale factor for diagnostic output [conc CU-1 ~> 1]
real :: dtr_iter_stop ! The tolerance for changes in tracer concentrations between smoothing
! iterations that determines when to stop iterating [CU ~> conc]
real :: npoints ! The number of points in an average [nondim]
Expand All @@ -355,10 +354,10 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion,
is_ongrid = .false.
if (present(ongrid)) is_ongrid = ongrid

dtr_iter_stop = 1.0e-3*conversion
dtr_iter_stop = 1.0e-3*scale
if (present(tr_iter_tol)) dtr_iter_stop = tr_iter_tol

I_scale = 1.0 / conversion
I_scale = 1.0 / scale

PI_180 = atan(1.0)/45.

Expand All @@ -371,6 +370,9 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion,

call cpu_clock_begin(id_clock_read)

! A note by MJH copied from elsewhere suggests that this code may be using the model connectivity
! (e.g., reentrant or tripolar) but should use the dataset's connectivity instead.

call get_var_axes_info(trim(filename), trim(varnam), axes_info)

if (allocated(z_in)) deallocate(z_in)
Expand Down Expand Up @@ -418,7 +420,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion,
if (.not. found_attr) call MOM_error(FATAL, &
"error finding missing value for " // trim(varnam) // &
" in file " // trim(filename) // " in hinterp_extrap")
missing_value = conversion * missing_val_in
missing_value = scale * missing_val_in

call read_attribute(trim(filename), "scale_factor", scale_factor, &
varname=trim(varnam), found=found_attr)
Expand Down Expand Up @@ -471,7 +473,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion,
do i=is,ie
if (abs(tr_in(i,j)-missing_val_in) > abs(roundoff*missing_val_in)) then
mask_in(i,j) = 1.0
tr_in(i,j) = (tr_in(i,j)*scale_factor+add_offset) * conversion
tr_in(i,j) = (tr_in(i,j)*scale_factor+add_offset) * scale
else
tr_in(i,j) = missing_value
endif
Expand Down Expand Up @@ -511,7 +513,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion,
do j=1,jdp ; do i=1,id
if (abs(tr_inp(i,j)-missing_val_in) > abs(roundoff*missing_val_in)) then
mask_in(i,j) = 1.0
tr_inp(i,j) = (tr_inp(i,j)*scale_factor+add_offset) * conversion
tr_inp(i,j) = (tr_inp(i,j)*scale_factor+add_offset) * scale
else
tr_inp(i,j) = missing_value
endif
Expand Down Expand Up @@ -596,18 +598,17 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion,
end subroutine horiz_interp_and_extrap_tracer_record

!> Extrapolate and interpolate using a FMS time interpolation handle
subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, tr_z, mask_z, &
z_in, z_edges_in, missing_value, reentrant_x, tripolar_n, &
subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, G, tr_z, mask_z, &
z_in, z_edges_in, missing_value, scale, &
homogenize, spongeOngrid, m_to_Z, &
answers_2018, tr_iter_tol, answer_date)

integer, intent(in) :: fms_id !< A unique id used by the FMS time interpolator
type(time_type), intent(in) :: Time !< A FMS time type
real, intent(in) :: conversion !< Conversion factor for tracer.
type(ocean_grid_type), intent(inout) :: G !< Grid object
real, allocatable, dimension(:,:,:), intent(out) :: tr_z
!< Allocatable tracer array on the horizontal
!! model grid and input-file vertical levels. [CU ~> conc]
!! model grid and input-file vertical levels [CU ~> conc]
real, allocatable, dimension(:,:,:), intent(out) :: mask_z
!< Allocatable tracer mask array on the horizontal
!! model grid and input-file vertical levels [nondim]
Expand All @@ -616,10 +617,10 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t
real, allocatable, dimension(:), intent(out) :: z_edges_in
!< Cell grid edge values for input data [Z ~> m]
real, intent(out) :: missing_value !< The missing value in the returned array, scaled
!! with conversion to avoid accidentally having valid
!! values match missing values [CU ~> conc]
logical, intent(in) :: reentrant_x !< If true, this grid is reentrant in the x-direction
logical, intent(in) :: tripolar_n !< If true, this is a northern tripolar grid
!! to avoid accidentally having valid values match
!! missing values [CU ~> conc]
real, intent(in) :: scale !< Scaling factor for tracer into the internal
!! units of the model [CU conc-1 ~> 1]
logical, optional, intent(in) :: homogenize !< If present and true, horizontally homogenize data
!! to produce perfectly "flat" initial conditions
logical, optional, intent(in) :: spongeOngrid !< If present and true, the sponge data are on the model grid
Expand Down Expand Up @@ -655,7 +656,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t
real, dimension(:), allocatable :: lon_in ! The longitudes in the input file [degreesE] then [radians]
real, dimension(:), allocatable :: lat_in ! The latitudes in the input file [degreesN] then [radians]
real, dimension(:), allocatable :: lat_inp ! The input file latitudes expanded to the pole [degreesN] then [radians]
real :: max_lat ! The maximum latitude on the input grid [degreeN]
real :: max_lat ! The maximum latitude on the input grid [degreesN]
real :: pole ! The sum of tracer values at the pole [conc]
real :: max_depth ! The maximum depth of the ocean [Z ~> m]
real :: npole ! The number of points contributing to the pole value [nondim]
Expand All @@ -673,7 +674,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t
logical :: debug=.false.
logical :: is_ongrid
integer :: ans_date ! The vintage of the expressions and order of arithmetic to use
real :: I_scale ! The inverse of the conversion factor for diagnostic output [conc CU-1 ~> 1]
real :: I_scale ! The inverse of the scale factor for diagnostic output [conc CU-1 ~> 1]
real :: dtr_iter_stop ! The tolerance for changes in tracer concentrations between smoothing
! iterations that determines when to stop iterating [CU ~> conc]
real :: npoints ! The number of points in an average [nondim]
Expand All @@ -699,10 +700,10 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t

id_clock_read = cpu_clock_id('(Initialize tracer from Z) read', grain=CLOCK_LOOP)

dtr_iter_stop = 1.0e-3*conversion
dtr_iter_stop = 1.0e-3*scale
if (present(tr_iter_tol)) dtr_iter_stop = tr_iter_tol

I_scale = 1.0 / conversion
I_scale = 1.0 / scale

PI_180 = atan(1.0)/45.

Expand All @@ -716,7 +717,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t
call cpu_clock_begin(id_clock_read)

call get_external_field_info(fms_id, size=fld_sz, axes=axes_data, missing=missing_val_in)
missing_value = conversion*missing_val_in
missing_value = scale*missing_val_in

verbosity = MOM_get_verbosity()

Expand Down Expand Up @@ -823,7 +824,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t
do j=1,jdp ; do i=1,id
if (abs(tr_inp(i,j)-missing_val_in) > abs(roundoff*missing_val_in)) then
mask_in(i,j) = 1.0
tr_inp(i,j) = tr_inp(i,j) * conversion
tr_inp(i,j) = tr_inp(i,j) * scale
else
tr_inp(i,j) = missing_value
endif
Expand Down Expand Up @@ -909,7 +910,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t
do k=1,kd
do j=js,je
do i=is,ie
tr_z(i,j,k) = data_in(i,j,k) * conversion
tr_z(i,j,k) = data_in(i,j,k) * scale
if (ans_date >= 20190101) mask_z(i,j,k) = 1.
if (abs(tr_z(i,j,k)-missing_value) < abs(roundoff*missing_value)) mask_z(i,j,k) = 0.
enddo
Expand Down
43 changes: 23 additions & 20 deletions src/initialization/MOM_state_initialization.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2483,8 +2483,12 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just
real :: PI_180 ! for conversion from degrees to radians [radian degree-1]
real :: Hmix_default ! The default initial mixed layer depth [Z ~> m].
real :: Hmix_depth ! The mixed layer depth in the initial condition [Z ~> m].
real :: missing_value_temp ! The missing value in the input temperature field
real :: missing_value_salt ! The missing value in the input salinity field
real :: missing_value_temp ! The missing value in the input temperature field [C ~> degC]
real :: missing_value_salt ! The missing value in the input salinity field [S ~> ppt]
real :: tol_temp ! The tolerance for changes in temperature during the horizontal
! interpolation from an input dataset [C ~> degC]
real :: tol_sal ! The tolerance for changes in salinity during the horizontal
! interpolation from an input dataset [S ~> ppt]
logical :: correct_thickness
real :: h_tolerance ! A parameter that controls the tolerance when adjusting the
! thickness to fit the bathymetry [Z ~> m].
Expand All @@ -2494,19 +2498,19 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just
logical :: adjust_temperature = .true. ! fit t/s to target densities
real :: temp_land_fill ! A temperature value to use for land points [C ~> degC]
real :: salt_land_fill ! A salinity value to use for land points [C ~> degC]
logical :: reentrant_x, tripolar_n

! data arrays
real, dimension(:), allocatable :: z_edges_in, z_in ! Interface heights [Z ~> m]
real, dimension(:), allocatable :: Rb ! Interface densities [R ~> kg m-3]
real, dimension(:), allocatable :: z_edges_in ! Input data interface heights or depths [Z ~> m]
real, dimension(:), allocatable :: z_in ! Input data cell heights or depths [Z ~> m]
real, dimension(:), allocatable :: Rb ! Interface densities [R ~> kg m-3]
real, dimension(:,:,:), allocatable, target :: temp_z ! Input temperatures [C ~> degC]
real, dimension(:,:,:), allocatable, target :: salt_z ! Input salinities [S ~> ppt]
real, dimension(:,:,:), allocatable, target :: mask_z ! 1 for valid data points [nondim]
real, dimension(:,:,:), allocatable :: rho_z ! Densities in Z-space [R ~> kg m-3]
real, dimension(:,:,:), allocatable :: rho_z ! Densities in Z-space [R ~> kg m-3]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: zi ! Interface heights [Z ~> m].
real, dimension(SZI_(G),SZJ_(G)) :: Z_bottom ! The (usually negative) height of the seafloor
! relative to the surface [Z ~> m].
integer, dimension(SZI_(G),SZJ_(G)) :: nlevs ! The number of levels in each column with valid data
real, dimension(SZI_(G),SZJ_(G)) :: Z_bottom ! The (usually negative) height of the seafloor
! relative to the surface [Z ~> m].
integer, dimension(SZI_(G),SZJ_(G)) :: nlevs ! The number of levels in each column with valid data
real, dimension(SZI_(G)) :: press ! Pressures [R L2 T-2 ~> Pa].

! Local variables for ALE remapping
Expand Down Expand Up @@ -2569,9 +2573,6 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just

eos => tv%eqn_of_state

reentrant_x = .false. ; call get_param(PF, mdl, "REENTRANT_X", reentrant_x, default=.true.)
tripolar_n = .false. ; call get_param(PF, mdl, "TRIPOLAR_N", tripolar_n, default=.false.)

call get_param(PF, mdl, "TEMP_SALT_Z_INIT_FILE", filename, &
"The name of the z-space input file used to initialize "//&
"temperatures (T) and salinities (S). If T and S are not "//&
Expand Down Expand Up @@ -2701,6 +2702,8 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just
!### These hard-coded constants should be made into runtime parameters
temp_land_fill = 0.0*US%degC_to_C
salt_land_fill = 35.0*US%ppt_to_S
tol_temp = 1.0e-3*US%degC_to_C
tol_sal = 1.0e-3*US%ppt_to_S

eps_z = GV%Angstrom_Z
eps_rho = 1.0e-10*US%kg_m3_to_R
Expand All @@ -2720,15 +2723,15 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just
! to the North/South Pole past the limits of the input data, they are extrapolated using the average
! value at the northernmost/southernmost latitude.

call horiz_interp_and_extrap_tracer(tfilename, potemp_var, US%degC_to_C, 1, &
G, temp_z, mask_z, z_in, z_edges_in, missing_value_temp, reentrant_x, &
tripolar_n, homogenize, m_to_Z=US%m_to_Z, answer_date=hor_regrid_answer_date, &
ongrid=pre_gridded, tr_iter_tol=1.0e-3*US%degC_to_C)
call horiz_interp_and_extrap_tracer(tfilename, potemp_var, 1, &
G, temp_z, mask_z, z_in, z_edges_in, missing_value_temp, &
scale=US%degC_to_C, homogenize=homogenize, m_to_Z=US%m_to_Z, &
answer_date=hor_regrid_answer_date, ongrid=pre_gridded, tr_iter_tol=tol_temp)

call horiz_interp_and_extrap_tracer(sfilename, salin_var, US%ppt_to_S, 1, &
G, salt_z, mask_z, z_in, z_edges_in, missing_value_salt, reentrant_x, &
tripolar_n, homogenize, m_to_Z=US%m_to_Z, answer_date=hor_regrid_answer_date, &
ongrid=pre_gridded, tr_iter_tol=1.0e-3*US%ppt_to_S)
call horiz_interp_and_extrap_tracer(sfilename, salin_var, 1, &
G, salt_z, mask_z, z_in, z_edges_in, missing_value_salt, &
scale=US%ppt_to_S, homogenize=homogenize, m_to_Z=US%m_to_Z, &
answer_date=hor_regrid_answer_date, ongrid=pre_gridded, tr_iter_tol=tol_sal)

kd = size(z_in,1)

Expand Down
Loading

0 comments on commit 0e2490d

Please sign in to comment.