Skip to content

Commit

Permalink
+tracer_Z_init can rescale tracers
Browse files Browse the repository at this point in the history 
  Added the new optional scale argument to tracer_Z_init to allow the units of
the tracers that it initializes from a z-space file to be rescaled, similarly to
what is already done in tracer_z_init_array.  This also required the addition of
the new missing_scale argument to the private routine read_Z_edges to specify
how the missing values are rescaled. Comments were also added or modified
describing 28 real variables or their units in this module.  By default all
answers are bitwise identical, but when this new option is exercised, answers
could change at roundoff for some non-Boussinesq cases.
  • Loading branch information
@Hallberg-NOAA @marshallward
Hallberg-NOAA authored and marshallward committed Feb 21, 2024
1 parent f14a681 commit 3f7465a
Showing 1 changed file with 54 additions and 35 deletions.
89 changes: 54 additions & 35 deletions src/tracer/MOM_tracer_Z_init.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,50 +27,58 @@ module MOM_tracer_Z_init

!> This function initializes a tracer by reading a Z-space file, returning
!! .true. if this appears to have been successful, and false otherwise.
function tracer_Z_init(tr, h, filename, tr_name, G, GV, US, missing_val, land_val)
function tracer_Z_init(tr, h, filename, tr_name, G, GV, US, missing_val, land_val, scale)
logical :: tracer_Z_init !< A return code indicating if the initialization has been successful
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
intent(out) :: tr !< The tracer to initialize
intent(out) :: tr !< The tracer to initialize [CU ~> conc]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
character(len=*), intent(in) :: filename !< The name of the file to read from
character(len=*), intent(in) :: tr_name !< The name of the tracer in the file
real, optional, intent(in) :: missing_val !< The missing value for the tracer
real, optional, intent(in) :: land_val !< A value to use to fill in land points

! This include declares and sets the variable "version".
# include "version_variable.h"
intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] or other
!! arbitrary units such as [Z ~> m]
character(len=*), intent(in) :: filename !< The name of the file to read from
character(len=*), intent(in) :: tr_name !< The name of the tracer in the file
real, optional, intent(in) :: missing_val !< The missing value for the tracer [CU ~> conc]
real, optional, intent(in) :: land_val !< A value to use to fill in land points [CU ~> conc]
real, optional, intent(in) :: scale !< A factor by which to scale the output tracers from the
!! their units in the file [CU conc-1 ~> 1]

! Local variables
real, allocatable, dimension(:,:,:) :: &
tr_in ! The z-space array of tracer concentrations that is read in.
tr_in ! The z-space array of tracer concentrations that is read in [CU ~> conc]
real, allocatable, dimension(:) :: &
z_edges, & ! The depths of the cell edges or cell centers (depending on
! the value of has_edges) in the input z* data [Z ~> m].
tr_1d, & ! A copy of the input tracer concentrations in a column.
tr_1d, & ! A copy of the input tracer concentrations in a column [CU ~> conc]
wt, & ! The fractional weight for each layer in the range between
! k_top and k_bot, nondim.
z1, & ! z1 and z2 are the depths of the top and bottom limits of the part
z2 ! of a z-cell that contributes to a layer, relative to the cell
! center and normalized by the cell thickness, nondim.
! k_top and k_bot [nondim]
z1, z2 ! z1 and z2 are the depths of the top and bottom limits of the part
! of a z-cell that contributes to a layer, relative to the cell
! center and normalized by the cell thickness [nondim].
! Note that -1/2 <= z1 <= z2 <= 1/2.
real :: e(SZK_(GV)+1) ! The z-star interface heights [Z ~> m].
real :: landval ! The tracer value to use in land points.
real :: landval ! The tracer value to use in land points [CU ~> conc]
real :: sl_tr ! The normalized slope of the tracer
! within the cell, in tracer units.
! within the cell, in tracer units [CU ~> conc]
real :: htot(SZI_(G)) ! The vertical sum of h [H ~> m or kg m-2].
real :: dilate ! The amount by which the thicknesses are dilated to
! create a z-star coordinate, nondim or in m3 kg-1.
real :: missing ! The missing value for the tracer.

! create a z-star coordinate [Z H-1 ~> nondim or m3 kg-1]
! or other units reflecting those of h
real :: missing ! The missing value for the tracer [CU ~> conc]
real :: scale_fac ! A factor by which to scale the output tracers from the units in the
! input file [CU conc-1 ~> 1]
! This include declares and sets the variable "version".
# include "version_variable.h"
logical :: has_edges, use_missing, zero_surface
character(len=80) :: loc_msg
integer :: k_top, k_bot, k_bot_prev, k_start
integer :: i, j, k, kz, is, ie, js, je, nz, nz_in

is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke

scale_fac = 1.0 ; if (present(scale)) then ; scale_fac = scale ; endif

landval = 0.0 ; if (present(land_val)) landval = land_val

zero_surface = .false. ! Make this false for errors to be fatal.
Expand All @@ -83,15 +91,15 @@ function tracer_Z_init(tr, h, filename, tr_name, G, GV, US, missing_val, land_va
! Find out the number of input levels and read the depth of the edges,
! also modifying their sign convention to be monotonically decreasing.
call read_Z_edges(filename, tr_name, z_edges, nz_in, has_edges, use_missing, &
missing, scale=US%m_to_Z)
missing, scale=US%m_to_Z, missing_scale=scale_fac)
if (nz_in < 1) then
tracer_Z_init = .false.
return
endif

allocate(tr_in(G%isd:G%ied,G%jsd:G%jed,nz_in), source=0.0)
allocate(tr_1d(nz_in), source=0.0)
call MOM_read_data(filename, tr_name, tr_in(:,:,:), G%Domain)
call MOM_read_data(filename, tr_name, tr_in(:,:,:), G%Domain, scale=scale_fac)

! Fill missing values from above? Use a "close" test to avoid problems
! from type-conversion rounoff.
Expand Down Expand Up @@ -297,9 +305,10 @@ subroutine tracer_z_init_array(tr_in, z_edges, nk_data, e, land_fill, G, nlay, n
real :: e_1d(nlay+1) ! A 1-d column of interface heights, in the same units as e [Z ~> m] or [m]
real :: sl_tr ! The tracer concentration slope times the layer thickness, in tracer units [B]
real :: wt(nk_data) ! The fractional weight for each layer in the range between z1 and z2 [nondim]
real :: z1(nk_data) ! z1 and z2 are the fractional depths of the top and bottom
real :: z2(nk_data) ! limits of the part of a z-cell that contributes to a layer, relative
! to the cell center and normalized by the cell thickness [nondim].
real :: z1(nk_data) ! The fractional depth of the top limit of the part of a z-cell that contributes to
! a layer, relative to the cell center and normalized by the cell thickness [nondim].
real :: z2(nk_data) ! The fractional depth of the bottom limit of the part of a z-cell that contributes to
! a layer, relative to the cell center and normalized by the cell thickness [nondim].
! Note that -1/2 <= z1 <= z2 <= 1/2.
real :: scale_fac ! A factor by which to scale the output tracers from the input tracers [B A-1 ~> 1]
integer :: k_top, k_bot, k_bot_prev, kstart
Expand Down Expand Up @@ -380,18 +389,21 @@ end subroutine tracer_z_init_array
!> This subroutine reads the vertical coordinate data for a field from a NetCDF file.
!! It also might read the missing value attribute for that same field.
subroutine read_Z_edges(filename, tr_name, z_edges, nz_out, has_edges, &
use_missing, missing, scale)
use_missing, missing, scale, missing_scale)
character(len=*), intent(in) :: filename !< The name of the file to read from.
character(len=*), intent(in) :: tr_name !< The name of the tracer in the file.
real, dimension(:), allocatable, &
intent(out) :: z_edges !< The depths of the vertical edges of the tracer array
intent(out) :: z_edges !< The depths of the vertical edges of the tracer array [Z ~> m]
integer, intent(out) :: nz_out !< The number of vertical layers in the tracer array
logical, intent(out) :: has_edges !< If true the values in z_edges are the edges of the
!! tracer cells, otherwise they are the cell centers
logical, intent(inout) :: use_missing !< If false on input, see whether the tracer has a
!! missing value, and if so return true
real, intent(inout) :: missing !< The missing value, if one has been found
real, intent(in) :: scale !< A scaling factor for z_edges into new units.
real, intent(inout) :: missing !< The missing value, if one has been found [CU ~> conc]
real, intent(in) :: scale !< A scaling factor for z_edges into new units [Z m-1 ~> 1]
real, intent(in) :: missing_scale !< A scaling factor to use to convert the
!! tracers and their missing value from the units in
!! the file into their internal units [CU conc-1 ~> 1]

! This subroutine reads the vertical coordinate data for a field from a
! NetCDF file. It also might read the missing value attribute for that same field.
Expand Down Expand Up @@ -419,6 +431,7 @@ subroutine read_Z_edges(filename, tr_name, z_edges, nz_out, has_edges, &

if (.not.use_missing) then ! Try to find the missing value from the dataset.
call read_attribute(filename, "missing_value", missing, varname=tr_name, found=use_missing, ncid_in=ncid)
if (use_missing) missing = missing * missing_scale
endif
! Find out if the Z-axis has an edges attribute
call read_attribute(filename, "edges", edge_name, varname=dim_names(3), found=has_edges, ncid_in=ncid)
Expand Down Expand Up @@ -475,8 +488,12 @@ subroutine find_overlap(e, Z_top, Z_bot, k_max, k_start, k_top, k_bot, wt, z1, z
real, dimension(:), intent(out) :: z2 !< Depths of the bottom limit of the part of
!! a layer that contributes to a depth level, relative to the cell center and normalized
!! by the cell thickness [nondim]. Note that -1/2 <= z1 < z2 <= 1/2.

! Local variables
real :: Ih, e_c, tot_wt, I_totwt
real :: Ih ! The inverse of the vertical distance across a layer, in the inverse of the units of e [Z-1 ~> m-1]
real :: e_c ! The height of the layer center, in the units of e [Z ~> m]
real :: tot_wt ! The sum of the thicknesses contributing to a layer [Z ~> m]
real :: I_totwt ! The Adcroft reciprocal of tot_wt [Z-1 ~> m-1]
integer :: k

wt(:) = 0.0 ; z1(:) = 0.0 ; z2(:) = 0.0 ; k_bot = k_max
Expand All @@ -494,6 +511,7 @@ subroutine find_overlap(e, Z_top, Z_bot, k_max, k_start, k_top, k_bot, wt, z1, z
z1(k) = (e_c - MIN(e(K), Z_top)) * Ih
z2(k) = (e_c - Z_bot) * Ih
else
! Note that in theis branch, wt temporarily has units of [Z ~> m]
wt(k) = MIN(e(K),Z_top) - e(K+1) ; tot_wt = wt(k) ! These are always > 0.
if (e(K) /= e(K+1)) then
z1(k) = (0.5*(e(K)+e(K+1)) - MIN(e(K), Z_top)) / (e(K)-e(K+1))
Expand All @@ -515,6 +533,7 @@ subroutine find_overlap(e, Z_top, Z_bot, k_max, k_start, k_top, k_bot, wt, z1, z
enddo

I_totwt = 0.0 ; if (tot_wt > 0.0) I_totwt = 1.0 / tot_wt
! This loop changes the units of wt from [Z ~> m] to [nondim].
do k=k_top,k_bot ; wt(k) = I_totwt*wt(k) ; enddo
endif

Expand All @@ -523,13 +542,13 @@ end subroutine find_overlap
!> This subroutine determines a limited slope for val to be advected with
!! a piecewise limited scheme.
function find_limited_slope(val, e, k) result(slope)
real, dimension(:), intent(in) :: val !< An column the values that are being interpolated.
real, dimension(:), intent(in) :: val !< A column of the values that are being interpolated, in arbitrary units [A]
real, dimension(:), intent(in) :: e !< A column's interface heights [Z ~> m] or other units.
integer, intent(in) :: k !< The layer whose slope is being determined.
real :: slope !< The normalized slope in the intracell distribution of val.
real :: slope !< The normalized slope in the intracell distribution of val [A]
! Local variables
real :: amn, cmn
real :: d1, d2
real :: amn, cmn ! Limited differences and curvatures in the values [A]
real :: d1, d2 ! Layer thicknesses, in the units of e [Z ~> m]

if ((val(k)-val(k-1)) * (val(k)-val(k+1)) >= 0.0) then
slope = 0.0 ! ; curvature = 0.0
Expand Down

0 comments on commit 3f7465a

Please sign in to comment.