+tracer_Z_init can rescale tracers

src/tracer/MOM_tracer_Z_init.F90

@@ -27,50 +27,58 @@
 
 !>   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,15 +91,15 @@
   ! 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.
@@ -297,9 +305,10 @@
   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 @@
 !> 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 @@
 
   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 @@
   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 @@
     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 @@
     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 @@
 !> 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