diff --git a/src/tracer/MOM_tracer_Z_init.F90 b/src/tracer/MOM_tracer_Z_init.F90 index c404c560f3..2cf0ba1efe 100644 --- a/src/tracer/MOM_tracer_Z_init.F90 +++ b/src/tracer/MOM_tracer_Z_init.F90 @@ -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. @@ -83,7 +91,7 @@ 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 @@ -91,7 +99,7 @@ function tracer_Z_init(tr, h, filename, tr_name, G, GV, US, missing_val, land_va 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. @@ -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 @@ -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. @@ -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) @@ -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 @@ -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)) @@ -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 @@ -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