MOM6
Data Types | Functions/Subroutines
mom_interface_heights Module Reference

Detailed Description

Functions for calculating interface heights, including free surface height.

Data Types

interface  find_eta
 Calculates the heights of sruface or all interfaces from layer thicknesses. More...
 

Functions/Subroutines

subroutine find_eta_3d (h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m)
 Calculates the heights of all interfaces between layers, using the appropriate form for consistency with the calculation of the pressure gradient forces. Additionally, these height may be dilated for consistency with the corresponding time-average quantity from the barotropic calculation. More...
 
subroutine find_eta_2d (h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m)
 Calculates the free surface height, using the appropriate form for consistency with the calculation of the pressure gradient forces. Additionally, the sea surface height may be adjusted for consistency with the corresponding time-average quantity from the barotropic calculation. More...
 

Function/Subroutine Documentation

◆ find_eta_2d()

subroutine mom_interface_heights::find_eta_2d ( real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h,
type(thermo_var_ptrs), intent(in)  tv,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(out)  eta,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in), optional  eta_bt,
integer, intent(in), optional  halo_size,
real, intent(in), optional  eta_to_m 
)
private

Calculates the free surface height, using the appropriate form for consistency with the calculation of the pressure gradient forces. Additionally, the sea surface height may be adjusted for consistency with the corresponding time-average quantity from the barotropic calculation.

Parameters
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]usA dimensional unit scaling type
[in]hLayer thicknesses [H ~> m or kg m-2]
[in]tvA structure pointing to various thermodynamic variables.
[out]etafree surface height relative to mean sea level (z=0) often [Z ~> m].
[in]eta_btoptional barotropic variable that gives the "correct" free surface height (Boussinesq) or total water column mass per unit area (non-Boussinesq) [H ~> m or kg m-2].
[in]halo_sizewidth of halo points on which to calculate eta.
[in]eta_to_mThe conversion factor from the units of eta to m; by default this is USZ_to_m.

Definition at line 148 of file MOM_interface_heights.F90.

149  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
150  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
151  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
152  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
153  type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
154  !! thermodynamic variables.
155  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: eta !< free surface height relative to
156  !! mean sea level (z=0) often [Z ~> m].
157  real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: eta_bt !< optional barotropic
158  !! variable that gives the "correct" free surface height (Boussinesq) or total
159  !! water column mass per unit area (non-Boussinesq) [H ~> m or kg m-2].
160  integer, optional, intent(in) :: halo_size !< width of halo points on
161  !! which to calculate eta.
162  real, optional, intent(in) :: eta_to_m !< The conversion factor from
163  !! the units of eta to m; by default this is US%Z_to_m.
164  ! Local variables
165  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &
166  p ! Hydrostatic pressure at each interface [R L2 T-2 ~> Pa]
167  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
168  dz_geo ! The change in geopotential height across a layer [L2 T-2 ~> m2 s-2].
169  real :: htot(SZI_(G)) ! The sum of all layers' thicknesses [H ~> m or kg m-2].
170  real :: I_gEarth ! The inverse of the gravitational acceleration times the
171  ! rescaling factor derived from eta_to_m [T2 Z L-2 ~> s2 m-1]
172  real :: Z_to_eta, H_to_eta, H_to_rho_eta ! Unit conversion factors with obvious names.
173  integer i, j, k, is, ie, js, je, nz, halo
174 
175  halo = 0 ; if (present(halo_size)) halo = max(0,halo_size)
176  is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
177  nz = g%ke
178 
179  z_to_eta = 1.0 ; if (present(eta_to_m)) z_to_eta = us%Z_to_m / eta_to_m
180  h_to_eta = gv%H_to_Z * z_to_eta
181  h_to_rho_eta = gv%H_to_RZ * z_to_eta
182  i_gearth = z_to_eta / gv%g_Earth
183 
184 !$OMP parallel default(shared) private(htot)
185 !$OMP do
186  do j=js,je ; do i=is,ie ; eta(i,j) = -z_to_eta*g%bathyT(i,j) ; enddo ; enddo
187 
188  if (gv%Boussinesq) then
189  if (present(eta_bt)) then
190 !$OMP do
191  do j=js,je ; do i=is,ie
192  eta(i,j) = h_to_eta*eta_bt(i,j)
193  enddo ; enddo
194  else
195 !$OMP do
196  do j=js,je ; do k=1,nz ; do i=is,ie
197  eta(i,j) = eta(i,j) + h(i,j,k)*h_to_eta
198  enddo ; enddo ; enddo
199  endif
200  else
201  if (associated(tv%eqn_of_state)) then
202 !$OMP do
203  do j=js,je
204  if (associated(tv%p_surf)) then
205  do i=is,ie ; p(i,j,1) = tv%p_surf(i,j) ; enddo
206  else
207  do i=is,ie ; p(i,j,1) = 0.0 ; enddo
208  endif
209 
210  do k=1,nz ; do i=is,ie
211  p(i,j,k+1) = p(i,j,k) + gv%g_Earth*gv%H_to_RZ*h(i,j,k)
212  enddo ; enddo
213  enddo
214 !$OMP do
215  do k = 1, nz
216  call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p(:,:,k), p(:,:,k+1), 0.0, &
217  g%HI, tv%eqn_of_state, us, dz_geo(:,:,k), halo_size=halo)
218  enddo
219 !$OMP do
220  do j=js,je ; do k=1,nz ; do i=is,ie
221  eta(i,j) = eta(i,j) + i_gearth * dz_geo(i,j,k)
222  enddo ; enddo ; enddo
223  else
224 !$OMP do
225  do j=js,je ; do k=1,nz ; do i=is,ie
226  eta(i,j) = eta(i,j) + h_to_rho_eta*h(i,j,k) / gv%Rlay(k)
227  enddo ; enddo ; enddo
228  endif
229  if (present(eta_bt)) then
230  ! Dilate the water column to agree with the time-averaged column
231  ! mass from the barotropic solution.
232 !$OMP do
233  do j=js,je
234  do i=is,ie ; htot(i) = gv%H_subroundoff ; enddo
235  do k=1,nz ; do i=is,ie ; htot(i) = htot(i) + h(i,j,k) ; enddo ; enddo
236  do i=is,ie
237  eta(i,j) = (eta_bt(i,j) / htot(i)) * (eta(i,j) + z_to_eta*g%bathyT(i,j)) - &
238  z_to_eta*g%bathyT(i,j)
239  enddo
240  enddo
241  endif
242  endif
243 !$OMP end parallel
244 

◆ find_eta_3d()

subroutine mom_interface_heights::find_eta_3d ( real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h,
type(thermo_var_ptrs), intent(in)  tv,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke+1), intent(out)  eta,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in), optional  eta_bt,
integer, intent(in), optional  halo_size,
real, intent(in), optional  eta_to_m 
)
private

Calculates the heights of all interfaces between layers, using the appropriate form for consistency with the calculation of the pressure gradient forces. Additionally, these height may be dilated for consistency with the corresponding time-average quantity from the barotropic calculation.

Parameters
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]usA dimensional unit scaling type
[in]hLayer thicknesses [H ~> m or kg m-2]
[in]tvA structure pointing to various thermodynamic variables.
[out]etalayer interface heights [Z ~> m] or 1/eta_to_m m).
[in]eta_btoptional barotropic variable that gives the "correct" free surface height (Boussinesq) or total water column mass per unit area (non-Boussinesq). This is used to dilate the layer. thicknesses when calculating interfaceheights [H ~> m or kg m-2].
[in]halo_sizewidth of halo points on which to calculate eta.
[in]eta_to_mThe conversion factor from the units of eta to m; by default this is USZ_to_m.

Definition at line 31 of file MOM_interface_heights.F90.

32  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
33  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
34  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
35  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
36  type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
37  !! thermodynamic variables.
38  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: eta !< layer interface heights
39  !! [Z ~> m] or 1/eta_to_m m).
40  real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: eta_bt !< optional barotropic
41  !! variable that gives the "correct" free surface height (Boussinesq) or total water
42  !! column mass per unit area (non-Boussinesq). This is used to dilate the layer.
43  !! thicknesses when calculating interfaceheights [H ~> m or kg m-2].
44  integer, optional, intent(in) :: halo_size !< width of halo points on
45  !! which to calculate eta.
46  real, optional, intent(in) :: eta_to_m !< The conversion factor from
47  !! the units of eta to m; by default this is US%Z_to_m.
48 
49  ! Local variables
50  real :: p(SZI_(G),SZJ_(G),SZK_(G)+1) ! Hydrostatic pressure at each interface [R L2 T-2 ~> Pa]
51  real :: dz_geo(SZI_(G),SZJ_(G),SZK_(G)) ! The change in geopotential height
52  ! across a layer [L2 T-2 ~> m2 s-2].
53  real :: dilate(SZI_(G)) ! non-dimensional dilation factor
54  real :: htot(SZI_(G)) ! total thickness [H ~> m or kg m-2]
55  real :: I_gEarth ! The inverse of the gravitational acceleration times the
56  ! rescaling factor derived from eta_to_m [T2 Z L-2 ~> s2 m-1]
57  real :: Z_to_eta, H_to_eta, H_to_rho_eta ! Unit conversion factors with obvious names.
58  integer i, j, k, isv, iev, jsv, jev, nz, halo
59 
60  halo = 0 ; if (present(halo_size)) halo = max(0,halo_size)
61 
62  isv = g%isc-halo ; iev = g%iec+halo ; jsv = g%jsc-halo ; jev = g%jec+halo
63  nz = g%ke
64 
65  if ((isv<g%isd) .or. (iev>g%ied) .or. (jsv<g%jsd) .or. (jev>g%jed)) &
66  call mom_error(fatal,"find_eta called with an overly large halo_size.")
67 
68  z_to_eta = 1.0 ; if (present(eta_to_m)) z_to_eta = us%Z_to_m / eta_to_m
69  h_to_eta = gv%H_to_Z * z_to_eta
70  h_to_rho_eta = gv%H_to_RZ * z_to_eta
71  i_gearth = z_to_eta / gv%g_Earth
72 
73 !$OMP parallel default(shared) private(dilate,htot)
74 !$OMP do
75  do j=jsv,jev ; do i=isv,iev ; eta(i,j,nz+1) = -z_to_eta*g%bathyT(i,j) ; enddo ; enddo
76 
77  if (gv%Boussinesq) then
78 !$OMP do
79  do j=jsv,jev ; do k=nz,1,-1; do i=isv,iev
80  eta(i,j,k) = eta(i,j,k+1) + h(i,j,k)*h_to_eta
81  enddo ; enddo ; enddo
82  if (present(eta_bt)) then
83  ! Dilate the water column to agree with the free surface height
84  ! that is used for the dynamics.
85 !$OMP do
86  do j=jsv,jev
87  do i=isv,iev
88  dilate(i) = (eta_bt(i,j)*h_to_eta + z_to_eta*g%bathyT(i,j)) / &
89  (eta(i,j,1) + z_to_eta*g%bathyT(i,j))
90  enddo
91  do k=1,nz ; do i=isv,iev
92  eta(i,j,k) = dilate(i) * (eta(i,j,k) + z_to_eta*g%bathyT(i,j)) - z_to_eta*g%bathyT(i,j)
93  enddo ; enddo
94  enddo
95  endif
96  else
97  if (associated(tv%eqn_of_state)) then
98 !$OMP do
99  do j=jsv,jev
100  if (associated(tv%p_surf)) then
101  do i=isv,iev ; p(i,j,1) = tv%p_surf(i,j) ; enddo
102  else
103  do i=isv,iev ; p(i,j,1) = 0.0 ; enddo
104  endif
105  do k=1,nz ; do i=isv,iev
106  p(i,j,k+1) = p(i,j,k) + gv%g_Earth*gv%H_to_RZ*h(i,j,k)
107  enddo ; enddo
108  enddo
109 !$OMP do
110  do k=1,nz
111  call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p(:,:,k), p(:,:,k+1), &
112  0.0, g%HI, tv%eqn_of_state, us, dz_geo(:,:,k), halo_size=halo)
113  enddo
114 !$OMP do
115  do j=jsv,jev
116  do k=nz,1,-1 ; do i=isv,iev
117  eta(i,j,k) = eta(i,j,k+1) + i_gearth * dz_geo(i,j,k)
118  enddo ; enddo
119  enddo
120  else
121 !$OMP do
122  do j=jsv,jev ; do k=nz,1,-1; do i=isv,iev
123  eta(i,j,k) = eta(i,j,k+1) + h_to_rho_eta*h(i,j,k) / gv%Rlay(k)
124  enddo ; enddo ; enddo
125  endif
126  if (present(eta_bt)) then
127  ! Dilate the water column to agree with the free surface height
128  ! from the time-averaged barotropic solution.
129 !$OMP do
130  do j=jsv,jev
131  do i=isv,iev ; htot(i) = gv%H_subroundoff ; enddo
132  do k=1,nz ; do i=isv,iev ; htot(i) = htot(i) + h(i,j,k) ; enddo ; enddo
133  do i=isv,iev ; dilate(i) = eta_bt(i,j) / htot(i) ; enddo
134  do k=1,nz ; do i=isv,iev
135  eta(i,j,k) = dilate(i) * (eta(i,j,k) + z_to_eta*g%bathyT(i,j)) - z_to_eta*g%bathyT(i,j)
136  enddo ; enddo
137  enddo
138  endif
139  endif
140 !$OMP end parallel
141