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 a few instances 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
   to uniform values derived from `MAXIMUM_DEPTH`
 - the setting of MLD_EN_VALS in MOM_diabatic_driver.F90 which was
   previously using a work around to set defaults to 25, 2500, 250000 J/m2.
 - two vectors of 4 values in user/user_change_diffusivity.F90

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

Author: adcroft

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,68 @@ 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")

src/parameterizations/vertical/MOM_diabatic_driver.F90

@@ -3256,13 +3256,8 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di
         'Mixed layer depth (delta rho = 0.125)', 'm', conversion=US%Z_to_m)
     call get_param(param_file, mdl, "MLD_EN_VALS", CS%MLD_En_vals, &
          "The energy values used to compute MLDs.  If not set (or all set to 0.), the "//&
-         "default will overwrite to 25., 2500., 250000.", &
-         units='J/m2', default=0., scale=US%W_m2_to_RZ3_T3*US%s_to_T)
-    if ((CS%MLD_En_vals(1)==0.).and.(CS%MLD_En_vals(2)==0.).and.(CS%MLD_En_vals(3)==0.)) then
-      CS%MLD_En_vals = (/ 25.*US%W_m2_to_RZ3_T3*US%s_to_T, &
-                        2500.*US%W_m2_to_RZ3_T3*US%s_to_T, &
-                      250000.*US%W_m2_to_RZ3_T3*US%s_to_T /)
-    endif
+         "default will overwrite to 25., 2500., 250000.", units='J/m2', &
+         defaults=(/25., 2500., 250000./), scale=US%W_m2_to_RZ3_T3*US%s_to_T)
     write(EN1,'(F10.2)') CS%MLD_En_vals(1)*US%RZ3_T3_to_W_m2*US%T_to_s
     write(EN2,'(F10.2)') CS%MLD_En_vals(2)*US%RZ3_T3_to_W_m2*US%T_to_s
     write(EN3,'(F10.2)') CS%MLD_En_vals(3)*US%RZ3_T3_to_W_m2*US%T_to_s

src/user/user_change_diffusivity.F90

@@ -230,14 +230,15 @@ subroutine user_change_diff_init(Time, G, GV, US, param_file, diag, CS)
                  "applied.  The four values specify the latitudes at "//&
                  "which the extra diffusivity starts to increase from 0, "//&
                  "hits its full value, starts to decrease again, and is "//&
-                 "back to 0.", units="degrees_N", default=-1.0e9)
+                 "back to 0.", units="degrees_N", defaults=(/-1.0e9,-1.0e9,-1.0e9,-1.0e9/))
     call get_param(param_file, mdl, "USER_KD_ADD_RHO_RANGE", CS%rho_range(:), &
                  "Four successive values that define a range of potential "//&
                  "densities over which the user-given extra diffusivity "//&
                  "is applied.  The four values specify the density at "//&
                  "which the extra diffusivity starts to increase from 0, "//&
                  "hits its full value, starts to decrease again, and is "//&
-                 "back to 0.", units="kg m-3", default=-1.0e9, scale=US%kg_m3_to_R)
+                 "back to 0.", units="kg m-3", defaults=(/-1.0e9,-1.0e9,-1.0e9,-1.0e9/),&
+                 scale=US%kg_m3_to_R)
     call get_param(param_file, mdl, "USER_KD_ADD_USE_ABS_LAT", CS%use_abs_lat, &
                  "If true, use the absolute value of latitude when "//&
                  "checking whether a point fits into range of latitudes.", &