(*)Refactor and document soliton_initialization

Refactored the code in soliton_initialization.F90 to more accurately reflect the nondimensionalization that was applied in developing this test case. This change includes reading in the maximum depth and beta, and using them to calculate the equatorial deformation radius and the external gravity wave speed. The references to the papers describing the test case in the module were added to the doxygen comments describing the routines in this module. This commit also includes adding comments documenting the nature and units of all the internal variables in this module. There are two new arguments each (param_file and just_read) to soliton_initialize_thickness and soliton_initialize_velocity to accommodate these changes, bringing them into line with the interfaces for other similar user initialization routines, and MOM_initialize_state was changed accordingly. This change could change answers in general, but in the specific examples that use this code, both beta and the external wave speed are deliberately set to 1 in MKS units, so this commit does not change answers for that specific case.
NOAA-GFDL · May 30, 2024 · d90ff6b · d90ff6b
1 parent 7305528
commit d90ff6b
Show file tree

Hide file tree

Showing 2 changed files with 94 additions and 19 deletions.
diff --git a/src/initialization/MOM_state_initialization.F90 b/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))

diff --git a/src/user/soliton_initialization.F90 b/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)