soliton_initialization.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Initial conditions for the Equatorial Rossby soliton test (Boyd).
7
8use mom_error_handler, only : mom_mesg, mom_error, fatal, is_root_pe
9use mom_file_parser, only : get_param, log_version, param_file_type
10use mom_get_input, only : directories
11use mom_grid, only : ocean_grid_type
15use regrid_consts, only : coordinatemode, default_coordinate_mode
16use regrid_consts, only : regridding_layer, regridding_zstar
17use regrid_consts, only : regridding_rho, regridding_sigma
18
19implicit none ; private
20
21#include <MOM_memory.h>
22
23! Private (module-wise) parameters
24character(len=40) :: mdl = "soliton_initialization" !< This module's name.
25
28
29contains
30
31!> Initialization of thicknesses in equatorial Rossby soliton test, as described in section
32!! 6.1 of Haidvogel and Beckman (1990) and in Boyd (1980, JPO) and Boyd (1985, JPO).
33subroutine soliton_initialize_thickness(h, depth_tot, G, GV, US, param_file, just_read)
34 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
35 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
36 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
37 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
38 intent(out) :: h !< The thickness that is being initialized [Z ~> m]
39 real, dimension(SZI_(G),SZJ_(G)), &
40 intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
41 type(param_file_type), intent(in) :: param_file !< A structure indicating the open file
42 !! to parse for model parameter values.
43 logical, intent(in) :: just_read !< If true, this call will only read
44 !! parameters without changing h.
45
46 ! Local variables
47 real :: max_depth ! Maximum depth of the model bathymetry [Z ~> m]
48 real :: cg_max ! The external wave speed based on max_depth [L T-1 ~> m s-1]
49 real :: beta ! The meridional gradient of the Coriolis parameter [T-1 L-1 ~> s-1 m-1]
50 real :: l_eq ! The equatorial deformation radius used in nondimensionalizing this problem [L ~> m]
51 real :: scale_pos ! A conversion factor to nondimensionalize the axis units, usually [m-1]
52 real :: x0 ! Initial x-position of the soliton in the same units as geoLonT, often [m].
53 real :: y0 ! Initial y-position of the soliton in the same units as geoLatT, often [m].
54 real :: x, y ! Nondimensionalized positions [nondim]
55 real :: i_nz ! The inverse of the number of layers [nondim]
56 real :: val1 ! A nondimensionlized zonal decay scale [nondim]
57 real :: val2 ! An overall surface height anomaly amplitude [L T-1 ~> m s-1]
58 real :: val3 ! A decay factor [nondim]
59 real :: val4 ! The local velocity amplitude [L T-1 ~> m s-1]
60 ! This include declares and sets the variable "version".
61# include "version_variable.h"
62 integer :: i, j, k, is, ie, js, je, nz
63
64 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
65
66 if (.not.just_read) &
67 call mom_mesg("soliton_initialization.F90, soliton_initialize_thickness: setting thickness")
68
69 if (.not.just_read) call log_version(param_file, mdl, version, "")
70 call get_param(param_file, mdl, "MAXIMUM_DEPTH", max_depth, &
71 units="m", default=-1.e9, scale=us%m_to_Z, do_not_log=.true.)
72 call get_param(param_file, mdl, "BETA", beta, &
73 "The northward gradient of the Coriolis parameter with the betaplane option.", &
74 units="m-1 s-1", default=0.0, scale=us%T_to_s*us%L_to_m, do_not_log=.true.)
75
76 if (just_read) return ! All run-time parameters have been read, so return.
77
78 if (max_depth <= 0.0) call mom_error(fatal, &
79 "soliton_initialization, soliton_initialize_thickness: "//&
80 "This module requires a positive value of MAXIMUM_DEPTH.")
81 if (abs(beta) <= 0.0) call mom_error(fatal, &
82 "soliton_initialization, soliton_initialize_thickness: "//&
83 "This module requires a non-zero value of BETA.")
84 if (g%grid_unit_to_L <= 0.) call mom_error(fatal, "soliton_initialization.F90: "//&
85 "soliton_initialize_thickness() is only set to work with Cartesian axis units.")
86
87 cg_max = sqrt(gv%g_Earth * max_depth)
88 l_eq = sqrt(cg_max / abs(beta))
89 scale_pos = g%grid_unit_to_L / l_eq
90 i_nz = 1.0 / real(nz)
91
92 x0 = 2.0*g%len_lon/3.0
93 y0 = 0.0
94 val1 = 0.395
95 val2 = max_depth * 0.771*(val1*val1)
96
97 do j = g%jsc,g%jec ; do i = g%isc,g%iec
98 do k = 1, nz
99 x = (g%geoLonT(i,j)-x0) * scale_pos
100 y = (g%geoLatT(i,j)-y0) * scale_pos
101 val3 = exp(-val1*x)
102 val4 = val2 * ( 2.0*val3 / (1.0 + (val3*val3)) )**2
103 h(i,j,k) = (0.25*val4*(6.0*y*y + 3.0) * exp(-0.5*y*y) + depth_tot(i,j)) * i_nz
104 enddo
105 enddo ; enddo
106
107end subroutine soliton_initialize_thickness
108
109
110!> Initialization of u and v in the equatorial Rossby soliton test, as described in section
111!! 6.1 of Haidvogel and Beckman (1990) and in Boyd (1980, JPO) and Boyd (1985, JPO).
112subroutine soliton_initialize_velocity(u, v, G, GV, US, param_file, just_read)
113 type(ocean_grid_type), intent(in) :: g !< Grid structure
114 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
115 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(out) :: u !< i-component of velocity [L T-1 ~> m s-1]
116 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(out) :: v !< j-component of velocity [L T-1 ~> m s-1]
117 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
118 type(param_file_type), intent(in) :: param_file !< A structure indicating the open file
119 !! to parse for model parameter values.
120 logical, intent(in) :: just_read !< If true, this call will only read
121 !! parameters without changing h.
122
123 ! Local variables
124 real :: max_depth ! Maximum depth of the model bathymetry [Z ~> m]
125 real :: cg_max ! The external wave speed based on max_depth [L T-1 ~> m s-1]
126 real :: beta ! The meridional gradient of the Coriolis parameter [T-1 L-1 ~> s-1 m-1]
127 real :: l_eq ! The equatorial deformation radius used in nondimensionalizing this problem [L ~> m]
128 real :: scale_pos ! A conversion factor to nondimensionalize the axis units, usually [m-1]
129 real :: x0 ! Initial x-position of the soliton in the same units as geoLonT, often [m].
130 real :: y0 ! Initial y-position of the soliton in the same units as geoLatT, often [m].
131 real :: x, y ! Nondimensionalized positions [nondim]
132 real :: val1 ! A nondimensionlized zonal decay scale [nondim]
133 real :: val2 ! An overall velocity amplitude [L T-1 ~> m s-1]
134 real :: val3 ! A decay factor [nondim]
135 real :: val4 ! The local velocity amplitude [L T-1 ~> m s-1]
136 integer :: i, j, k, is, ie, js, je, nz
137
138 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
139
140 if (.not.just_read) &
141 call mom_mesg("soliton_initialization.F90, soliton_initialize_thickness: setting thickness")
142
143 call get_param(param_file, mdl, "MAXIMUM_DEPTH", max_depth, &
144 units="m", default=-1.e9, scale=us%m_to_Z, do_not_log=.true.)
145 call get_param(param_file, mdl, "BETA", beta, &
146 "The northward gradient of the Coriolis parameter with the betaplane option.", &
147 units="m-1 s-1", default=0.0, scale=us%T_to_s*us%L_to_m, do_not_log=.true.)
148
149 if (just_read) return ! All run-time parameters have been read, so return.
150
151 if (max_depth <= 0.0) call mom_error(fatal, &
152 "soliton_initialization, soliton_initialize_velocity: "//&
153 "This module requires a positive value of MAXIMUM_DEPTH.")
154 if (abs(beta) <= 0.0) call mom_error(fatal, &
155 "soliton_initialization, soliton_initialize_velocity: "//&
156 "This module requires a non-zero value of BETA.")
157 if (g%grid_unit_to_L <= 0.) call mom_error(fatal, "soliton_initialization.F90: "//&
158 "soliton_initialize_velocity() is only set to work with Cartesian axis units.")
159
160 cg_max = sqrt(gv%g_Earth * max_depth)
161 l_eq = sqrt(cg_max / abs(beta))
162 scale_pos = g%grid_unit_to_L / l_eq
163
164 x0 = 2.0*g%len_lon/3.0
165 y0 = 0.0
166 val1 = 0.395
167 val2 = cg_max * 0.771*(val1*val1)
168
169 v(:,:,:) = 0.0
170 u(:,:,:) = 0.0
171
172 do j = g%jsc,g%jec ; do i = g%isc-1,g%iec+1
173 do k = 1, nz
174 x = (0.5*(g%geoLonT(i+1,j)+g%geoLonT(i,j))-x0) * scale_pos
175 y = (0.5*(g%geoLatT(i+1,j)+g%geoLatT(i,j))-y0) * scale_pos
176 val3 = exp(-val1*x)
177 val4 = val2*((2.0*val3/(1.0+(val3*val3)))**2)
178 u(i,j,k) = 0.25*val4*(6.0*y*y-9.0) * exp(-0.5*y*y)
179 enddo
180 enddo ; enddo
181 do j = g%jsc-1,g%jec+1 ; do i = g%isc,g%iec
182 do k = 1, nz
183 x = 0.5*(g%geoLonT(i,j+1)+g%geoLonT(i,j))-x0
184 y = 0.5*(g%geoLatT(i,j+1)+g%geoLatT(i,j))-y0
185 val3 = exp(-val1*x)
186 val4 = val2*((2.0*val3/(1.0+(val3*val3)))**2)
187 v(i,j,k) = 2.0*val4*y*(-2.0*val1*tanh(val1*x)) * exp(-0.5*y*y)
188 enddo
189 enddo ; enddo
190
191end subroutine soliton_initialize_velocity
192
193
194!> \namespace soliton_initialization
195!!
196!! \section section_soliton Description of the equatorial Rossby soliton initial
197!! conditions
198!!
199
200end module soliton_initialization