Adds a vector of default values to get_param_real_array()

The `default=` optional argument to get_param() only provides a uniform
value to initialize an array of reals. This commit adds the optional
`defaults=` argument that must have the same length as the `values`
argument.

I've also added an instance of this optional argument by adding the
`initialize_thickness_param()` procedure, selected by `THICKNESS_CONFIG
= "param"`. The procedure was based on the "uniform" method, and uses
the parameter `THICKNESS_INIT_VALUES` which defaults iniform values
derived from `MAXIMUM_DEPTH`

I modified tc2 to use `THICKNESS_INIT_VALUES="param"` instead of
`THICKNESS_INIT_VALUES="uniform"`.

There will be some doc file changes, but no answer changes.

Author: adcroft

.testing/tc2/MOM_input

@@ -193,7 +193,7 @@ MINIMUM_DEPTH = 1.0             !   [m] default = 0.0
                                 ! The minimum depth of the ocean.
 MAXIMUM_DEPTH = 4500.0          !   [m]
                                 ! The maximum depth of the ocean.
-THICKNESS_CONFIG = "uniform"    !
+THICKNESS_CONFIG = "param"    !
                                 ! A string that determines how the initial layer
                                 ! thicknesses are specified for a new run:
                                 !     file - read interface heights from the file specified

src/framework/MOM_document.F90

@@ -303,14 +303,16 @@ subroutine doc_param_real(doc, varname, desc, units, val, default, debuggingPara
 end subroutine doc_param_real
 
 !> This subroutine handles parameter documentation for arrays of reals.
-subroutine doc_param_real_array(doc, varname, desc, units, vals, default, debuggingParam, like_default)
+subroutine doc_param_real_array(doc, varname, desc, units, vals, default, defaults, &
+                                debuggingParam, like_default)
   type(doc_type),    pointer    :: doc     !< A pointer to a structure that controls where the
                                            !! documentation occurs and its formatting
   character(len=*),  intent(in) :: varname !< The name of the parameter being documented
   character(len=*),  intent(in) :: desc    !< A description of the parameter being documented
   character(len=*),  intent(in) :: units   !< The units of the parameter being documented
   real,              intent(in) :: vals(:) !< The array of values to record
-  real,    optional, intent(in) :: default !< The default value of this parameter
+  real,    optional, intent(in) :: default !< A uniform default value of this parameter
+  real,    optional, intent(in) :: defaults(:) !< The element-wise default values of this parameter
   logical, optional, intent(in) :: debuggingParam !< If present and true, this is a debugging parameter.
   logical, optional, intent(in) :: like_default !< If present and true, log this parameter as though
                                            !! it has the default value, even if there is no default.
@@ -334,6 +336,11 @@ subroutine doc_param_real_array(doc, varname, desc, units, vals, default, debugg
       do i=1,size(vals) ; if (vals(i) /= default) equalsDefault = .false. ; enddo
       mesg = trim(mesg)//" default = "//trim(real_string(default))
     endif
+    if (present(defaults)) then
+      equalsDefault = .true.
+      do i=1,size(vals) ; if (vals(i) /= defaults(i)) equalsDefault = .false. ; enddo
+      mesg = trim(mesg)//" default = "//trim(real_array_string(defaults))
+    endif
     if (present(like_default)) then ; if (like_default) equalsDefault = .true. ; endif
 
     if (mesgHasBeenDocumented(doc, varName, mesg)) return ! Avoid duplicates

src/framework/MOM_file_parser.F90

@@ -1464,7 +1464,7 @@ end subroutine log_param_real
 
 !> Log the name and values of an array of real model parameter in documentation files.
 subroutine log_param_real_array(CS, modulename, varname, value, desc, &
-                                units, default, debuggingParam, like_default, unscale)
+                                units, default, defaults, debuggingParam, like_default, unscale)
   type(param_file_type),      intent(in) :: CS      !< The control structure for the file_parser module,
                                          !! it is also a structure to parse for run-time parameters
   character(len=*),           intent(in) :: modulename !< The name of the calling module
@@ -1473,7 +1473,8 @@ subroutine log_param_real_array(CS, modulename, varname, value, desc, &
   character(len=*), optional, intent(in) :: desc    !< A description of this variable; if not
                                              !! present, this parameter is not written to a doc file
   character(len=*),           intent(in) :: units   !< The units of this parameter
-  real,             optional, intent(in) :: default !< The default value of the parameter
+  real,             optional, intent(in) :: default !< A uniform default value of the parameter
+  real,             optional, intent(in) :: defaults(:) !< The element-wise defaults of the parameter
   logical,          optional, intent(in) :: debuggingParam !< If present and true, this parameter is
                                          !! logged in the debugging parameter file
   logical,          optional, intent(in) :: like_default !< If present and true, log this parameter as
@@ -1498,7 +1499,7 @@ subroutine log_param_real_array(CS, modulename, varname, value, desc, &
 
   write(myunits(1:240),'(A)') trim(units)
   if (present(desc)) &
-    call doc_param(CS%doc, varname, desc, myunits, log_val, default, &
+    call doc_param(CS%doc, varname, desc, myunits, log_val, default, defaults, &
                    debuggingParam=debuggingParam, like_default=like_default)
 
 end subroutine log_param_real_array
@@ -1835,7 +1836,7 @@ end subroutine get_param_real
 !> This subroutine reads the values of an array of real model parameters from a parameter file
 !! and logs them in documentation files.
 subroutine get_param_real_array(CS, modulename, varname, value, desc, units, &
-               default, fail_if_missing, do_not_read, do_not_log, debuggingParam, &
+               default, defaults, fail_if_missing, do_not_read, do_not_log, debuggingParam, &
                scale, unscaled)
   type(param_file_type),      intent(in)    :: CS      !< The control structure for the file_parser module,
                                          !! it is also a structure to parse for run-time parameters
@@ -1846,7 +1847,8 @@ subroutine get_param_real_array(CS, modulename, varname, value, desc, units, &
   character(len=*), optional, intent(in)    :: desc    !< A description of this variable; if not
                                          !! present, this parameter is not written to a doc file
   character(len=*),           intent(in)    :: units   !< The units of this parameter
-  real,             optional, intent(in)    :: default !< The default value of the parameter
+  real,             optional, intent(in)    :: default !< A uniform default value of the parameter
+  real,             optional, intent(in)    :: defaults(:) !< The element-wise defaults of the parameter
   logical,          optional, intent(in)    :: fail_if_missing !< If present and true, a fatal error occurs
                                          !! if this variable is not found in the parameter file
   logical,          optional, intent(in)    :: do_not_read  !< If present and true, do not read a
@@ -1865,14 +1867,22 @@ subroutine get_param_real_array(CS, modulename, varname, value, desc, units, &
   do_read = .true. ; if (present(do_not_read)) do_read = .not.do_not_read
   do_log  = .true. ; if (present(do_not_log))  do_log  = .not.do_not_log
 
+  if (present(defaults)) then
+    if (present(default)) call MOM_error(FATAL, &
+          "get_param_real_array: Only one of default and defaults can be specified at a time.")
+    if (size(defaults) /= size(value)) call MOM_error(FATAL, &
+          "get_param_real_array: The size of defaults nad value are not the same.")
+  endif
+
   if (do_read) then
     if (present(default)) value(:) = default
+    if (present(defaults)) value(:) = defaults(:)
     call read_param_real_array(CS, varname, value, fail_if_missing)
   endif
 
   if (do_log) then
     call log_param_real_array(CS, modulename, varname, value, desc, &
-                              units, default, debuggingParam)
+                              units, default, defaults, debuggingParam)
   endif
 
   if (present(unscaled)) unscaled(:) = value(:)

src/initialization/MOM_state_initialization.F90

@@ -278,6 +278,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
              " \t uniform - uniform thickness layers evenly distributed \n"//&
              " \t\t between the surface and MAXIMUM_DEPTH. \n"//&
              " \t list - read a list of positive interface depths. \n"//&
+             " \t param - use thicknesses from parameter THICKNESS_INIT_VALUES. \n"//&
              " \t DOME - use a slope and channel configuration for the \n"//&
              " \t\t DOME sill-overflow test case. \n"//&
              " \t ISOMIP - use a configuration for the \n"//&
@@ -318,6 +319,8 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
                                  just_read=just_read)
       case ("list"); call initialize_thickness_list(dz, depth_tot, G, GV, US, PF, &
                                  just_read=just_read)
+      case ("param"); call initialize_thickness_param(dz, depth_tot, G, GV, US, PF, &
+                                 just_read=just_read)
       case ("DOME"); call DOME_initialize_thickness(dz, depth_tot, G, GV, PF, &
                               just_read=just_read)
       case ("ISOMIP"); call ISOMIP_initialize_thickness(dz, depth_tot, G, GV, US, PF, tv, &
@@ -1011,6 +1014,69 @@ subroutine initialize_thickness_list(h, depth_tot, G, GV, US, param_file, just_r
   call callTree_leave(trim(mdl)//'()')
 end subroutine initialize_thickness_list
 
+!> Initializes thickness based on a run-time parameter with nominal thickness
+!! for each layer
+subroutine initialize_thickness_param(h, depth_tot, G, GV, US, param_file, just_read)
+  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) :: h           !< The thickness that is being initialized [Z ~> m]
+  real, dimension(SZI_(G),SZJ_(G)), &
+                           intent(in)  :: depth_tot   !< The nominal total depth of the ocean [Z ~> m]
+  type(param_file_type),   intent(in)  :: param_file  !< A structure indicating the open file
+                                                      !! to parse for model parameter values.
+  logical,                 intent(in)  :: just_read   !< If true, this call will only read
+                                                      !! parameters without changing h.
+  ! Local variables
+  character(len=40)  :: mdl = "initialize_thickness_param" ! This subroutine's name.
+  real :: e0(SZK_(GV)+1)  ! The resting interface heights [Z ~> m], usually
+                          ! negative because it is positive upward.
+  real :: eta1D(SZK_(GV)+1)! Interface height relative to the sea surface,
+                          ! positive upward [Z ~> m].
+  real :: dz(SZK_(GV))    ! The nominal initial layer thickness [Z ~> m], usually
+  real :: h0_def(SZK_(GV)) ! Uniform default values for dz [Z ~> m], usually
+  integer :: i, j, k, is, ie, js, je, nz
+
+  call callTree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
+  if (G%max_depth<=0.) call MOM_error(FATAL, "initialize_thickness_param: "// &
+      "MAXIMUM_DEPTH has a nonsensical value! Was it set?")
+
+  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
+
+
+  h0_def = ( G%max_depth / real(nz) ) * US%Z_to_m
+  call get_param(param_file, mdl, "THICKNESS_INIT_VALUES", dz, &
+                 "A list of nominal thickness for each layer to initialize with", &
+                 units="m", scale=US%m_to_Z, defaults=h0_def, do_not_log=just_read)
+  if (just_read) return ! This subroutine has no run-time parameters.
+
+  e0(nz+1) = -G%max_depth
+  do k=nz, 1, -1
+    e0(K) = e0(K+1) + dz(k)
+  enddo
+
+  do j=js,je ; do i=is,ie
+    ! This sets the initial thickness (in m) of the layers.  The
+    ! thicknesses are set to insure that: 1.  each layer is at least an
+    ! Angstrom thick, and 2.  the interfaces are where they should be
+    ! based on the resting depths and interface height perturbations,
+    ! as long at this doesn't interfere with 1.
+    eta1D(nz+1) = -depth_tot(i,j)
+    do k=nz,1,-1
+      eta1D(K) = e0(K)
+      if (eta1D(K) < (eta1D(K+1) + GV%Angstrom_Z)) then
+        eta1D(K) = eta1D(K+1) + GV%Angstrom_Z
+        h(i,j,k) = GV%Angstrom_Z
+      else
+        h(i,j,k) = eta1D(K) - eta1D(K+1)
+      endif
+    enddo
+  enddo ; enddo
+
+  call callTree_leave(trim(mdl)//'()')
+end subroutine initialize_thickness_param
+
 !> Search density space for location of layers (not implemented!)
 subroutine initialize_thickness_search
   call MOM_error(FATAL,"  MOM_state_initialization.F90, initialize_thickness_search: NOT IMPLEMENTED")