(*)Refactor and document soliton_initialization

src/initialization/MOM_state_initialization.F90

@@ -343,7 +343,8 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
                                   just_read=just_read)
       case ("dumbbell"); call dumbbell_initialize_thickness(dz, depth_tot, G, GV, US, PF, &
                                   just_read=just_read)
-      case ("soliton"); call soliton_initialize_thickness(dz, depth_tot, G, GV, US)
+      case ("soliton"); call soliton_initialize_thickness(dz, depth_tot, G, GV, US, PF, &
+                                  just_read=just_read)
       case ("phillips"); call Phillips_initialize_thickness(dz, depth_tot, G, GV, US, PF, &
                                   just_read=just_read)
       case ("rossby_front")
@@ -508,7 +509,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
     case ("phillips"); call Phillips_initialize_velocity(u, v, G, GV, US, PF, just_read)
     case ("rossby_front"); call Rossby_front_initialize_velocity(u, v, h, &
                                      G, GV, US, PF, just_read)
-    case ("soliton"); call soliton_initialize_velocity(u, v, G, GV, US)
+    case ("soliton"); call soliton_initialize_velocity(u, v, G, GV, US, PF, just_read)
     case ("USER"); call user_initialize_velocity(u, v, G, GV, US, PF, just_read)
     case default ; call MOM_error(FATAL,  "MOM_initialize_state: "//&
           "Unrecognized velocity configuration "//trim(config))

src/user/soliton_initialization.F90

@@ -26,73 +26,147 @@ module soliton_initialization
 
 contains
 
-!> Initialization of thicknesses in Equatorial Rossby soliton test
-subroutine soliton_initialize_thickness(h, depth_tot, G, GV, US)
+!> Initialization of thicknesses in equatorial Rossby soliton test, as described in section
+!! 6.1 of Haidvogel and Beckman (1990) and in Boyd (1980, JPO) and Boyd (1985, JPO).
+subroutine soliton_initialize_thickness(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
+  real    :: max_depth  ! Maximum depth of the model bathymetry [Z ~> m]
+  real    :: cg_max     ! The external wave speed based on max_depth [L T-1 ~> m s-1]
+  real    :: beta       ! The meridional gradient of the Coriolis parameter [T-1 L-1 ~> s-1 m-1]
+  real    :: L_eq       ! The equatorial deformation radius used in nondimensionalizing this problem [L ~> m]
+  real    :: scale_pos  ! A conversion factor to nondimensionalize the axis units, usually [m-1]
+  real    :: x0    ! Initial x-position of the soliton in the same units as geoLonT, often [m].
+  real    :: y0    ! Initial y-position of the soliton in the same units as geoLatT, often [m].
+  real    :: x, y  ! Nondimensionalized positions [nondim]
+  real    :: I_nz  ! The inverse of the number of layers [nondim]
+  real    :: val1  ! A nondimensionlized zonal decay scale [nondim]
+  real    :: val2  ! An overall surface height anomaly amplitude [L T-1 ~> m s-1]
+  real    :: val3  ! A decay factor [nondim]
+  real    :: val4  ! The local velocity amplitude [L T-1 ~> m s-1]
+  ! This include declares and sets the variable "version".
+# include "version_variable.h"
   integer :: i, j, k, is, ie, js, je, nz
-  real    :: x, y, x0, y0
-  real    :: val1, val2, val3, val4
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
 
-  call MOM_mesg("soliton_initialization.F90, soliton_initialize_thickness: setting thickness")
+  if (.not.just_read) &
+    call MOM_mesg("soliton_initialization.F90, soliton_initialize_thickness: setting thickness")
+
+  if (.not.just_read) call log_version(param_file, mdl, version, "")
+  call get_param(param_file, mdl, "MAXIMUM_DEPTH", max_depth, &
+                 units="m", default=-1.e9, scale=US%m_to_Z, do_not_log=.true.)
+  call get_param(param_file, mdl, "BETA", beta, &
+                 "The northward gradient of the Coriolis parameter with the betaplane option.", &
+                 units="m-1 s-1", default=0.0, scale=US%T_to_s*US%L_to_m, do_not_log=.true.)
+
+  if (just_read) return ! All run-time parameters have been read, so return.
+
+  if (max_depth <= 0.0) call MOM_error(FATAL, &
+      "soliton_initialization, soliton_initialize_thickness: "//&
+      "This module requires a positive value of MAXIMUM_DEPTH.")
+  if (abs(beta) <= 0.0) call MOM_error(FATAL, &
+      "soliton_initialization, soliton_initialize_thickness: "//&
+      "This module requires a non-zero value of BETA.")
+
+  cg_max = sqrt(GV%g_Earth * max_depth)
+  L_eq = sqrt(cg_max / abs(beta))
+  scale_pos = US%m_to_L / L_eq
+  I_nz = 1.0 / real(nz)
 
   x0 = 2.0*G%len_lon/3.0
   y0 = 0.0
   val1 = 0.395
-  val2 = US%m_to_Z * 0.771*(val1*val1)
+  val2 = max_depth * 0.771*(val1*val1)
 
   do j = G%jsc,G%jec ; do i = G%isc,G%iec
     do k = 1, nz
-      x = G%geoLonT(i,j)-x0
-      y = G%geoLatT(i,j)-y0
+      x = (G%geoLonT(i,j)-x0) * scale_pos
+      y = (G%geoLatT(i,j)-y0) * scale_pos
       val3 = exp(-val1*x)
       val4 = val2 * ( 2.0*val3 / (1.0 + (val3*val3)) )**2
-      h(i,j,k) = (0.25*val4*(6.0*y*y + 3.0) * exp(-0.5*y*y) + depth_tot(i,j))
+      h(i,j,k) = (0.25*val4*(6.0*y*y + 3.0) * exp(-0.5*y*y) + depth_tot(i,j)) * I_nz
     enddo
   enddo ; enddo
 
 end subroutine soliton_initialize_thickness
 
 
-!> Initialization of u and v in the equatorial Rossby soliton test
-subroutine soliton_initialize_velocity(u, v, G, GV, US)
+!> Initialization of u and v in the equatorial Rossby soliton test, as described in section
+!! 6.1 of Haidvogel and Beckman (1990) and in Boyd (1980, JPO) and Boyd (1985, JPO).
+subroutine soliton_initialize_velocity(u, v, G, GV, US, param_file, just_read)
   type(ocean_grid_type),                      intent(in)  :: G  !< Grid structure
   type(verticalGrid_type),                    intent(in)  :: GV !< The ocean's vertical grid structure
   real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(out) :: u  !< i-component of velocity [L T-1 ~> m s-1]
   real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(out) :: v  !< j-component of velocity [L T-1 ~> m s-1]
   type(unit_scale_type),                      intent(in)  :: US !< A dimensional unit scaling type
+  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
-  real    :: x, x0 ! Positions in the same units as geoLonT.
-  real    :: y, y0 ! Positions in the same units as geoLatT.
-  real    :: val1  ! A zonal decay scale in the inverse of the units of geoLonT.
+  real    :: max_depth  ! Maximum depth of the model bathymetry [Z ~> m]
+  real    :: cg_max     ! The external wave speed based on max_depth [L T-1 ~> m s-1]
+  real    :: beta       ! The meridional gradient of the Coriolis parameter [T-1 L-1 ~> s-1 m-1]
+  real    :: L_eq       ! The equatorial deformation radius used in nondimensionalizing this problem [L ~> m]
+  real    :: scale_pos  ! A conversion factor to nondimensionalize the axis units, usually [m-1]
+  real    :: x0    ! Initial x-position of the soliton in the same units as geoLonT, often [m].
+  real    :: y0    ! Initial y-position of the soliton in the same units as geoLatT, often [m].
+  real    :: x, y  ! Nondimensionalized positions [nondim]
+  real    :: val1  ! A nondimensionlized zonal decay scale [nondim]
   real    :: val2  ! An overall velocity amplitude [L T-1 ~> m s-1]
   real    :: val3  ! A decay factor [nondim]
   real    :: val4  ! The local velocity amplitude [L T-1 ~> m s-1]
   integer :: i, j, k, is, ie, js, je, nz
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
 
+  if (.not.just_read) &
+    call MOM_mesg("soliton_initialization.F90, soliton_initialize_thickness: setting thickness")
+
+  call get_param(param_file, mdl, "MAXIMUM_DEPTH", max_depth, &
+                 units="m", default=-1.e9, scale=US%m_to_Z, do_not_log=.true.)
+  call get_param(param_file, mdl, "BETA", beta, &
+                 "The northward gradient of the Coriolis parameter with the betaplane option.", &
+                 units="m-1 s-1", default=0.0, scale=US%T_to_s*US%L_to_m, do_not_log=.true.)
+
+  if (just_read) return ! All run-time parameters have been read, so return.
+
+  if (max_depth <= 0.0) call MOM_error(FATAL, &
+      "soliton_initialization, soliton_initialize_velocity: "//&
+      "This module requires a positive value of MAXIMUM_DEPTH.")
+  if (abs(beta) <= 0.0) call MOM_error(FATAL, &
+      "soliton_initialization, soliton_initialize_velocity: "//&
+      "This module requires a non-zero value of BETA.")
+
+  cg_max = sqrt(GV%g_Earth * max_depth)
+  L_eq = sqrt(cg_max / abs(beta))
+  scale_pos = US%m_to_L / L_eq
+
   x0 = 2.0*G%len_lon/3.0
   y0 = 0.0
   val1 = 0.395
-  val2 = US%m_s_to_L_T * 0.771*(val1*val1)
+  val2 = cg_max * 0.771*(val1*val1)
 
   v(:,:,:) = 0.0
   u(:,:,:) = 0.0
 
   do j = G%jsc,G%jec ; do I = G%isc-1,G%iec+1
     do k = 1, nz
-      x = 0.5*(G%geoLonT(i+1,j)+G%geoLonT(i,j))-x0
-      y = 0.5*(G%geoLatT(i+1,j)+G%geoLatT(i,j))-y0
+      x = (0.5*(G%geoLonT(i+1,j)+G%geoLonT(i,j))-x0) * scale_pos
+      y = (0.5*(G%geoLatT(i+1,j)+G%geoLatT(i,j))-y0) * scale_pos
       val3 = exp(-val1*x)
       val4 = val2*((2.0*val3/(1.0+(val3*val3)))**2)
       u(I,j,k) = 0.25*val4*(6.0*y*y-9.0) * exp(-0.5*y*y)