MOM6
Functions/Subroutines
mom_isopycnal_slopes Module Reference

Detailed Description

Calculations of isoneutral slopes and stratification.

Functions/Subroutines

subroutine, public calc_isoneutral_slopes (G, GV, US, h, e, tv, dt_kappa_smooth, slope_x, slope_y, N2_u, N2_v, halo, OBC)
 Calculate isopycnal slopes, and optionally return N2 used in calculation. More...
 
subroutine, public vert_fill_ts (h, T_in, S_in, kappa_dt, T_f, S_f, G, GV, halo_here, larger_h_denom)
 Returns tracer arrays (nominally T and S) with massless layers filled with sensible values, by diffusing vertically with a small but constant diffusivity. More...
 

Function/Subroutine Documentation

◆ calc_isoneutral_slopes()

subroutine, public mom_isopycnal_slopes::calc_isoneutral_slopes ( 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), intent(in)  h,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke+1), intent(in)  e,
type(thermo_var_ptrs), intent(in)  tv,
real, intent(in)  dt_kappa_smooth,
real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke+1), intent(inout)  slope_x,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke+1), intent(inout)  slope_y,
real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke+1), intent(inout), optional  N2_u,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke+1), intent(inout), optional  N2_v,
integer, intent(in), optional  halo,
type(ocean_obc_type), optional, pointer  OBC 
)

Calculate isopycnal slopes, and optionally return N2 used in 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]eInterface heights [Z ~> m] or units given by 1/eta_to_m)
[in]tvA structure pointing to various thermodynamic variables
[in]dt_kappa_smoothA smoothing vertical diffusivity times a smoothing timescale [Z2 ~> m2].
[in,out]slope_xIsopycnal slope in i-direction [nondim]
[in,out]slope_yIsopycnal slope in j-direction [nondim]
[in,out]n2_uBrunt-Vaisala frequency squared at
[in,out]n2_vBrunt-Vaisala frequency squared at
[in]haloHalo width over which to compute
obcOpen boundaries control structure.

Definition at line 30 of file MOM_isopycnal_slopes.F90.

30  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
31  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
32  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
33  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
34  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in) :: e !< Interface heights [Z ~> m] or units
35  !! given by 1/eta_to_m)
36  type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
37  !! thermodynamic variables
38  real, intent(in) :: dt_kappa_smooth !< A smoothing vertical diffusivity
39  !! times a smoothing timescale [Z2 ~> m2].
40  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: slope_x !< Isopycnal slope in i-direction [nondim]
41  real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), intent(inout) :: slope_y !< Isopycnal slope in j-direction [nondim]
42  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), &
43  optional, intent(inout) :: N2_u !< Brunt-Vaisala frequency squared at
44  !! interfaces between u-points [T-2 ~> s-2]
45  real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), &
46  optional, intent(inout) :: N2_v !< Brunt-Vaisala frequency squared at
47  !! interfaces between u-points [T-2 ~> s-2]
48  integer, optional, intent(in) :: halo !< Halo width over which to compute
49  type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure.
50 
51  ! real, optional, intent(in) :: eta_to_m !< The conversion factor from the units
52  ! (This argument has been tested but for now serves no purpose.) !! of eta to m; US%Z_to_m by default.
53  ! Local variables
54  real, dimension(SZI_(G), SZJ_(G), SZK_(G)) :: &
55  T, & ! The temperature [degC], with the values in
56  ! in massless layers filled vertically by diffusion.
57  s !, & ! The filled salinity [ppt], with the values in
58  ! in massless layers filled vertically by diffusion.
59 ! Rho ! Density itself, when a nonlinear equation of state is not in use [R ~> kg m-3].
60  real, dimension(SZI_(G), SZJ_(G), SZK_(G)+1) :: &
61  pres ! The pressure at an interface [R L2 T-2 ~> Pa].
62  real, dimension(SZIB_(G)) :: &
63  drho_dT_u, & ! The derivative of density with temperature at u points [R degC-1 ~> kg m-3 degC-1].
64  drho_dS_u ! The derivative of density with salinity at u points [R ppt-1 ~> kg m-3 ppt-1].
65  real, dimension(SZI_(G)) :: &
66  drho_dT_v, & ! The derivative of density with temperature at v points [R degC-1 ~> kg m-3 degC-1].
67  drho_dS_v ! The derivative of density with salinity at v points [R ppt-1 ~> kg m-3 ppt-1].
68  real, dimension(SZIB_(G)) :: &
69  T_u, & ! Temperature on the interface at the u-point [degC].
70  S_u, & ! Salinity on the interface at the u-point [ppt].
71  pres_u ! Pressure on the interface at the u-point [R L2 T-2 ~> Pa].
72  real, dimension(SZI_(G)) :: &
73  T_v, & ! Temperature on the interface at the v-point [degC].
74  S_v, & ! Salinity on the interface at the v-point [ppt].
75  pres_v ! Pressure on the interface at the v-point [R L2 T-2 ~> Pa].
76  real :: drdiA, drdiB ! Along layer zonal- and meridional- potential density
77  real :: drdjA, drdjB ! gradients in the layers above (A) and below (B) the
78  ! interface times the grid spacing [R ~> kg m-3].
79  real :: drdkL, drdkR ! Vertical density differences across an interface [R ~> kg m-3].
80  real :: hg2A, hg2B ! Squares of geometric mean thicknesses [H2 ~> m2 or kg2 m-4].
81  real :: hg2L, hg2R ! Squares of geometric mean thicknesses [H2 ~> m2 or kg2 m-4].
82  real :: haA, haB, haL, haR ! Arithmetic mean thicknesses [H ~> m or kg m-2].
83  real :: dzaL, dzaR ! Temporary thicknesses in eta units [Z ~> m].
84  real :: wtA, wtB, wtL, wtR ! Unscaled weights, with various units.
85  real :: drdx, drdy ! Zonal and meridional density gradients [R L-1 ~> kg m-4].
86  real :: drdz ! Vertical density gradient [R Z-1 ~> kg m-4].
87  real :: Slope ! The slope of density surfaces, calculated in a way
88  ! that is always between -1 and 1.
89  real :: mag_grad2 ! The squared magnitude of the 3-d density gradient [R2 L-2 ~> kg2 m-8].
90  real :: slope2_Ratio ! The ratio of the slope squared to slope_max squared.
91  real :: h_neglect ! A thickness that is so small it is usually lost
92  ! in roundoff and can be neglected [H ~> m or kg m-2].
93  real :: h_neglect2 ! h_neglect^2 [H2 ~> m2 or kg2 m-4].
94  real :: dz_neglect ! A change in interface heighs that is so small it is usually lost
95  ! in roundoff and can be neglected [Z ~> m].
96  logical :: use_EOS ! If true, density is calculated from T & S using an equation of state.
97  real :: G_Rho0 ! The gravitational acceleration divided by density [Z2 T-2 R-1 ~> m5 kg-2 s-2]
98  real :: Z_to_L ! A conversion factor between from units for e to the
99  ! units for lateral distances.
100  real :: L_to_Z ! A conversion factor between from units for lateral distances
101  ! to the units for e.
102  real :: H_to_Z ! A conversion factor from thickness units to the units of e.
103 
104  logical :: present_N2_u, present_N2_v
105  integer, dimension(2) :: EOSdom_u, EOSdom_v ! Domains for the equation of state calculations at u and v points
106  integer :: is, ie, js, je, nz, IsdB
107  integer :: i, j, k
108  integer :: l_seg
109  logical :: local_open_u_BC, local_open_v_BC
110 
111  if (present(halo)) then
112  is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
113  else
114  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
115  endif
116  nz = g%ke ; isdb = g%IsdB
117 
118  h_neglect = gv%H_subroundoff ; h_neglect2 = h_neglect**2
119  z_to_l = us%Z_to_L ; h_to_z = gv%H_to_Z
120  ! if (present(eta_to_m)) then
121  ! Z_to_L = eta_to_m*US%m_to_L ; H_to_Z = GV%H_to_m / eta_to_m
122  ! endif
123  l_to_z = 1.0 / z_to_l
124  dz_neglect = gv%H_subroundoff * h_to_z
125 
126  local_open_u_bc = .false.
127  local_open_v_bc = .false.
128  if (present(obc)) then ; if (associated(obc)) then
129  local_open_u_bc = obc%open_u_BCs_exist_globally
130  local_open_v_bc = obc%open_v_BCs_exist_globally
131  endif ; endif
132 
133  use_eos = associated(tv%eqn_of_state)
134 
135  present_n2_u = PRESENT(n2_u)
136  present_n2_v = PRESENT(n2_v)
137  g_rho0 = (us%L_to_Z*l_to_z*gv%g_Earth) / gv%Rho0
138  if (present_n2_u) then
139  do j=js,je ; do i=is-1,ie
140  n2_u(i,j,1) = 0.
141  n2_u(i,j,nz+1) = 0.
142  enddo ; enddo
143  endif
144  if (present_n2_v) then
145  do j=js-1,je ; do i=is,ie
146  n2_v(i,j,1) = 0.
147  n2_v(i,j,nz+1) = 0.
148  enddo ; enddo
149  endif
150 
151  if (use_eos) then
152  if (present(halo)) then
153  call vert_fill_ts(h, tv%T, tv%S, dt_kappa_smooth, t, s, g, gv, halo+1)
154  else
155  call vert_fill_ts(h, tv%T, tv%S, dt_kappa_smooth, t, s, g, gv, 1)
156  endif
157  endif
158 
159  ! Find the maximum and minimum permitted streamfunction.
160  if (associated(tv%p_surf)) then
161  !$OMP parallel do default(shared)
162  do j=js-1,je+1 ; do i=is-1,ie+1
163  pres(i,j,1) = tv%p_surf(i,j)
164  enddo ; enddo
165  else
166  !$OMP parallel do default(shared)
167  do j=js-1,je+1 ; do i=is-1,ie+1
168  pres(i,j,1) = 0.0
169  enddo ; enddo
170  endif
171  !$OMP parallel do default(shared)
172  do j=js-1,je+1
173  do k=1,nz ; do i=is-1,ie+1
174  pres(i,j,k+1) = pres(i,j,k) + gv%g_Earth * gv%H_to_RZ * h(i,j,k)
175  enddo ; enddo
176  enddo
177 
178  eosdom_u(1) = is-1 - (g%IsdB-1) ; eosdom_u(2) = ie - (g%IsdB-1)
179 
180  !$OMP parallel do default(none) shared(nz,is,ie,js,je,IsdB,use_EOS,G,GV,US,pres,T,S,tv,h,e, &
181  !$OMP h_neglect,dz_neglect,Z_to_L,L_to_Z,H_to_Z,h_neglect2, &
182  !$OMP present_N2_u,G_Rho0,N2_u,slope_x,EOSdom_u,local_open_u_BC, &
183  !$OMP OBC) &
184  !$OMP private(drdiA,drdiB,drdkL,drdkR,pres_u,T_u,S_u, &
185  !$OMP drho_dT_u,drho_dS_u,hg2A,hg2B,hg2L,hg2R,haA, &
186  !$OMP haB,haL,haR,dzaL,dzaR,wtA,wtB,wtL,wtR,drdz, &
187  !$OMP drdx,mag_grad2,Slope,slope2_Ratio,l_seg)
188  do j=js,je ; do k=nz,2,-1
189  if (.not.(use_eos)) then
190  drdia = 0.0 ; drdib = 0.0
191  drdkl = gv%Rlay(k)-gv%Rlay(k-1) ; drdkr = gv%Rlay(k)-gv%Rlay(k-1)
192  endif
193 
194  ! Calculate the zonal isopycnal slope.
195  if (use_eos) then
196  do i=is-1,ie
197  pres_u(i) = 0.5*(pres(i,j,k) + pres(i+1,j,k))
198  t_u(i) = 0.25*((t(i,j,k) + t(i+1,j,k)) + (t(i,j,k-1) + t(i+1,j,k-1)))
199  s_u(i) = 0.25*((s(i,j,k) + s(i+1,j,k)) + (s(i,j,k-1) + s(i+1,j,k-1)))
200  enddo
201  call calculate_density_derivs(t_u, s_u, pres_u, drho_dt_u, drho_ds_u, &
202  tv%eqn_of_state, eosdom_u)
203  endif
204 
205  do i=is-1,ie
206  if (use_eos) then
207  ! Estimate the horizontal density gradients along layers.
208  drdia = drho_dt_u(i) * (t(i+1,j,k-1)-t(i,j,k-1)) + &
209  drho_ds_u(i) * (s(i+1,j,k-1)-s(i,j,k-1))
210  drdib = drho_dt_u(i) * (t(i+1,j,k)-t(i,j,k)) + &
211  drho_ds_u(i) * (s(i+1,j,k)-s(i,j,k))
212 
213  ! Estimate the vertical density gradients times the grid spacing.
214  drdkl = (drho_dt_u(i) * (t(i,j,k)-t(i,j,k-1)) + &
215  drho_ds_u(i) * (s(i,j,k)-s(i,j,k-1)))
216  drdkr = (drho_dt_u(i) * (t(i+1,j,k)-t(i+1,j,k-1)) + &
217  drho_ds_u(i) * (s(i+1,j,k)-s(i+1,j,k-1)))
218  endif
219 
220 
221  if (use_eos) then
222  hg2a = h(i,j,k-1)*h(i+1,j,k-1) + h_neglect2
223  hg2b = h(i,j,k)*h(i+1,j,k) + h_neglect2
224  hg2l = h(i,j,k-1)*h(i,j,k) + h_neglect2
225  hg2r = h(i+1,j,k-1)*h(i+1,j,k) + h_neglect2
226  haa = 0.5*(h(i,j,k-1) + h(i+1,j,k-1))
227  hab = 0.5*(h(i,j,k) + h(i+1,j,k)) + h_neglect
228  hal = 0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect
229  har = 0.5*(h(i+1,j,k-1) + h(i+1,j,k)) + h_neglect
230  if (gv%Boussinesq) then
231  dzal = hal * h_to_z ; dzar = har * h_to_z
232  else
233  dzal = 0.5*(e(i,j,k-1) - e(i,j,k+1)) + dz_neglect
234  dzar = 0.5*(e(i+1,j,k-1) - e(i+1,j,k+1)) + dz_neglect
235  endif
236  ! Use the harmonic mean thicknesses to weight the horizontal gradients.
237  ! These unnormalized weights have been rearranged to minimize divisions.
238  wta = hg2a*hab ; wtb = hg2b*haa
239  wtl = hg2l*(har*dzar) ; wtr = hg2r*(hal*dzal)
240 
241  drdz = (wtl * drdkl + wtr * drdkr) / (dzal*wtl + dzar*wtr)
242  ! The expression for drdz above is mathematically equivalent to:
243  ! drdz = ((hg2L/haL) * drdkL/dzaL + (hg2R/haR) * drdkR/dzaR) / &
244  ! ((hg2L/haL) + (hg2R/haR))
245  ! This is the gradient of density along geopotentials.
246  drdx = ((wta * drdia + wtb * drdib) / (wta + wtb) - &
247  drdz * (e(i,j,k)-e(i+1,j,k))) * g%IdxCu(i,j)
248 
249  ! This estimate of slope is accurate for small slopes, but bounded
250  ! to be between -1 and 1.
251  mag_grad2 = drdx**2 + (l_to_z*drdz)**2
252  if (mag_grad2 > 0.0) then
253  slope_x(i,j,k) = drdx / sqrt(mag_grad2)
254  else ! Just in case mag_grad2 = 0 ever.
255  slope_x(i,j,k) = 0.0
256  endif
257 
258  if (present_n2_u) n2_u(i,j,k) = g_rho0 * drdz * g%mask2dCu(i,j) ! Square of buoyancy frequency [T-2 ~> s-2]
259 
260  else ! With .not.use_EOS, the layers are constant density.
261  slope_x(i,j,k) = (z_to_l*(e(i,j,k)-e(i+1,j,k))) * g%IdxCu(i,j)
262  endif
263  if (local_open_u_bc) then
264  l_seg = obc%segnum_u(i,j)
265  if (l_seg /= obc_none) then
266  if (obc%segment(l_seg)%open) then
267  slope_x(i,j,k) = 0.
268  ! This and/or the masking code below is to make slopes match inside
269  ! land mask. Might not be necessary except for DEBUG output.
270 ! if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_E) then
271 ! slope_x(I+1,j,K) = 0.
272 ! else
273 ! slope_x(I-1,j,K) = 0.
274 ! endif
275  endif
276  endif
277  slope_x(i,j,k) = slope_x(i,j,k) * max(g%mask2dT(i,j),g%mask2dT(i+1,j))
278  endif
279 
280  enddo ! I
281  enddo ; enddo ! end of j-loop
282 
283  eosdom_v(1) = is - (g%isd-1) ; eosdom_v(2) = ie - (g%isd-1)
284 
285  ! Calculate the meridional isopycnal slope.
286  !$OMP parallel do default(none) shared(nz,is,ie,js,je,IsdB,use_EOS,G,GV,US,pres,T,S,tv, &
287  !$OMP h,h_neglect,e,dz_neglect,Z_to_L,L_to_Z,H_to_Z, &
288  !$OMP h_neglect2,present_N2_v,G_Rho0,N2_v,slope_y,EOSdom_v, &
289  !$OMP local_open_v_BC,OBC) &
290  !$OMP private(drdjA,drdjB,drdkL,drdkR,pres_v,T_v,S_v, &
291  !$OMP drho_dT_v,drho_dS_v,hg2A,hg2B,hg2L,hg2R,haA, &
292  !$OMP haB,haL,haR,dzaL,dzaR,wtA,wtB,wtL,wtR,drdz, &
293  !$OMP drdy,mag_grad2,Slope,slope2_Ratio,l_seg)
294  do j=js-1,je ; do k=nz,2,-1
295  if (.not.(use_eos)) then
296  drdja = 0.0 ; drdjb = 0.0
297  drdkl = gv%Rlay(k)-gv%Rlay(k-1) ; drdkr = gv%Rlay(k)-gv%Rlay(k-1)
298  endif
299 
300  if (use_eos) then
301  do i=is,ie
302  pres_v(i) = 0.5*(pres(i,j,k) + pres(i,j+1,k))
303  t_v(i) = 0.25*((t(i,j,k) + t(i,j+1,k)) + (t(i,j,k-1) + t(i,j+1,k-1)))
304  s_v(i) = 0.25*((s(i,j,k) + s(i,j+1,k)) + (s(i,j,k-1) + s(i,j+1,k-1)))
305  enddo
306  call calculate_density_derivs(t_v, s_v, pres_v, drho_dt_v, drho_ds_v, tv%eqn_of_state, &
307  eosdom_v)
308  endif
309  do i=is,ie
310  if (use_eos) then
311  ! Estimate the horizontal density gradients along layers.
312  drdja = drho_dt_v(i) * (t(i,j+1,k-1)-t(i,j,k-1)) + &
313  drho_ds_v(i) * (s(i,j+1,k-1)-s(i,j,k-1))
314  drdjb = drho_dt_v(i) * (t(i,j+1,k)-t(i,j,k)) + &
315  drho_ds_v(i) * (s(i,j+1,k)-s(i,j,k))
316 
317  ! Estimate the vertical density gradients times the grid spacing.
318  drdkl = (drho_dt_v(i) * (t(i,j,k)-t(i,j,k-1)) + &
319  drho_ds_v(i) * (s(i,j,k)-s(i,j,k-1)))
320  drdkr = (drho_dt_v(i) * (t(i,j+1,k)-t(i,j+1,k-1)) + &
321  drho_ds_v(i) * (s(i,j+1,k)-s(i,j+1,k-1)))
322  endif
323 
324  if (use_eos) then
325  hg2a = h(i,j,k-1)*h(i,j+1,k-1) + h_neglect2
326  hg2b = h(i,j,k)*h(i,j+1,k) + h_neglect2
327  hg2l = h(i,j,k-1)*h(i,j,k) + h_neglect2
328  hg2r = h(i,j+1,k-1)*h(i,j+1,k) + h_neglect2
329  haa = 0.5*(h(i,j,k-1) + h(i,j+1,k-1)) + h_neglect
330  hab = 0.5*(h(i,j,k) + h(i,j+1,k)) + h_neglect
331  hal = 0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect
332  har = 0.5*(h(i,j+1,k-1) + h(i,j+1,k)) + h_neglect
333  if (gv%Boussinesq) then
334  dzal = hal * h_to_z ; dzar = har * h_to_z
335  else
336  dzal = 0.5*(e(i,j,k-1) - e(i,j,k+1)) + dz_neglect
337  dzar = 0.5*(e(i,j+1,k-1) - e(i,j+1,k+1)) + dz_neglect
338  endif
339  ! Use the harmonic mean thicknesses to weight the horizontal gradients.
340  ! These unnormalized weights have been rearranged to minimize divisions.
341  wta = hg2a*hab ; wtb = hg2b*haa
342  wtl = hg2l*(har*dzar) ; wtr = hg2r*(hal*dzal)
343 
344  drdz = (wtl * drdkl + wtr * drdkr) / (dzal*wtl + dzar*wtr)
345  ! The expression for drdz above is mathematically equivalent to:
346  ! drdz = ((hg2L/haL) * drdkL/dzaL + (hg2R/haR) * drdkR/dzaR) / &
347  ! ((hg2L/haL) + (hg2R/haR))
348  ! This is the gradient of density along geopotentials.
349  drdy = ((wta * drdja + wtb * drdjb) / (wta + wtb) - &
350  drdz * (e(i,j,k)-e(i,j+1,k))) * g%IdyCv(i,j)
351 
352  ! This estimate of slope is accurate for small slopes, but bounded
353  ! to be between -1 and 1.
354  mag_grad2 = drdy**2 + (l_to_z*drdz)**2
355  if (mag_grad2 > 0.0) then
356  slope_y(i,j,k) = drdy / sqrt(mag_grad2)
357  else ! Just in case mag_grad2 = 0 ever.
358  slope_y(i,j,k) = 0.0
359  endif
360 
361  if (present_n2_v) n2_v(i,j,k) = g_rho0 * drdz * g%mask2dCv(i,j) ! Square of buoyancy frequency [T-2 ~> s-2]
362 
363  else ! With .not.use_EOS, the layers are constant density.
364  slope_y(i,j,k) = (z_to_l*(e(i,j,k)-e(i,j+1,k))) * g%IdyCv(i,j)
365  endif
366  if (local_open_v_bc) then
367  l_seg = obc%segnum_v(i,j)
368  if (l_seg /= obc_none) then
369  if (obc%segment(l_seg)%open) then
370  slope_y(i,j,k) = 0.
371  ! This and/or the masking code below is to make slopes match inside
372  ! land mask. Might not be necessary except for DEBUG output.
373 ! if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_N) then
374 ! slope_y(i,J+1,K) = 0.
375 ! else
376 ! slope_y(i,J-1,K) = 0.
377 ! endif
378  endif
379  endif
380  slope_y(i,j,k) = slope_y(i,j,k) * max(g%mask2dT(i,j),g%mask2dT(i,j+1))
381  endif
382 
383  enddo ! i
384  enddo ; enddo ! end of j-loop
385 

◆ vert_fill_ts()

subroutine, public mom_isopycnal_slopes::vert_fill_ts ( real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  T_in,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  S_in,
real, intent(in)  kappa_dt,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(out)  T_f,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(out)  S_f,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
integer, intent(in), optional  halo_here,
logical, intent(in), optional  larger_h_denom 
)

Returns tracer arrays (nominally T and S) with massless layers filled with sensible values, by diffusing vertically with a small but constant diffusivity.

Parameters
[in]gThe ocean's grid structure
[in]gvThe ocean's vertical grid structure
[in]hLayer thicknesses [H ~> m or kg m-2]
[in]t_inInput temperature [degC]
[in]s_inInput salinity [ppt]
[in]kappa_dtA vertical diffusivity to use for smoothing times a smoothing timescale [Z2 ~> m2].
[out]t_fFilled temperature [degC]
[out]s_fFilled salinity [ppt]
[in]halo_hereNumber of halo points to work on, 0 by default
[in]larger_h_denomPresent and true, add a large enough minimal thickness in the denominator of the flux calculations so that the fluxes are never so large as eliminate the transmission of information across groups of massless layers.

Definition at line 391 of file MOM_isopycnal_slopes.F90.

391  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
392  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
393  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
394  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: T_in !< Input temperature [degC]
395  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: S_in !< Input salinity [ppt]
396  real, intent(in) :: kappa_dt !< A vertical diffusivity to use for smoothing
397  !! times a smoothing timescale [Z2 ~> m2].
398  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out) :: T_f !< Filled temperature [degC]
399  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out) :: S_f !< Filled salinity [ppt]
400  integer, optional, intent(in) :: halo_here !< Number of halo points to work on,
401  !! 0 by default
402  logical, optional, intent(in) :: larger_h_denom !< Present and true, add a large
403  !! enough minimal thickness in the denominator of
404  !! the flux calculations so that the fluxes are
405  !! never so large as eliminate the transmission
406  !! of information across groups of massless layers.
407  ! Local variables
408  real :: ent(SZI_(G),SZK_(G)+1) ! The diffusive entrainment (kappa*dt)/dz
409  ! between layers in a timestep [H ~> m or kg m-2].
410  real :: b1(SZI_(G)), d1(SZI_(G)) ! b1, c1, and d1 are variables used by the
411  real :: c1(SZI_(G),SZK_(G)) ! tridiagonal solver.
412  real :: kap_dt_x2 ! The 2*kappa_dt converted to H units [H2 ~> m2 or kg2 m-4].
413  real :: h_neglect ! A negligible thickness [H ~> m or kg m-2], to allow for zero thicknesses.
414  real :: h0 ! A negligible thickness to allow for zero thickness layers without
415  ! completely decouping groups of layers [H ~> m or kg m-2].
416  ! Often 0 < h_neglect << h0.
417  real :: h_tr ! h_tr is h at tracer points with a tiny thickness
418  ! added to ensure positive definiteness [H ~> m or kg m-2].
419  integer :: i, j, k, is, ie, js, je, nz, halo
420 
421  halo=0 ; if (present(halo_here)) halo = max(halo_here,0)
422 
423  is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo ; nz = gv%ke
424 
425  h_neglect = gv%H_subroundoff
426  kap_dt_x2 = (2.0*kappa_dt)*gv%Z_to_H**2
427  h0 = h_neglect
428  if (present(larger_h_denom)) then
429  if (larger_h_denom) h0 = 1.0e-16*sqrt(kappa_dt)*gv%Z_to_H
430  endif
431 
432  if (kap_dt_x2 <= 0.0) then
433  !$OMP parallel do default(shared)
434  do k=1,nz ; do j=js,je ; do i=is,ie
435  t_f(i,j,k) = t_in(i,j,k) ; s_f(i,j,k) = s_in(i,j,k)
436  enddo ; enddo ; enddo
437  else
438  !$OMP parallel do default(shared) private(ent,b1,d1,c1,h_tr)
439  do j=js,je
440  do i=is,ie
441  ent(i,2) = kap_dt_x2 / ((h(i,j,1)+h(i,j,2)) + h0)
442  h_tr = h(i,j,1) + h_neglect
443  b1(i) = 1.0 / (h_tr + ent(i,2))
444  d1(i) = b1(i) * h_tr
445  t_f(i,j,1) = (b1(i)*h_tr)*t_in(i,j,1)
446  s_f(i,j,1) = (b1(i)*h_tr)*s_in(i,j,1)
447  enddo
448  do k=2,nz-1 ; do i=is,ie
449  ent(i,k+1) = kap_dt_x2 / ((h(i,j,k)+h(i,j,k+1)) + h0)
450  h_tr = h(i,j,k) + h_neglect
451  c1(i,k) = ent(i,k) * b1(i)
452  b1(i) = 1.0 / ((h_tr + d1(i)*ent(i,k)) + ent(i,k+1))
453  d1(i) = b1(i) * (h_tr + d1(i)*ent(i,k))
454  t_f(i,j,k) = b1(i) * (h_tr*t_in(i,j,k) + ent(i,k)*t_f(i,j,k-1))
455  s_f(i,j,k) = b1(i) * (h_tr*s_in(i,j,k) + ent(i,k)*s_f(i,j,k-1))
456  enddo ; enddo
457  do i=is,ie
458  c1(i,nz) = ent(i,nz) * b1(i)
459  h_tr = h(i,j,nz) + h_neglect
460  b1(i) = 1.0 / (h_tr + d1(i)*ent(i,nz))
461  t_f(i,j,nz) = b1(i) * (h_tr*t_in(i,j,nz) + ent(i,nz)*t_f(i,j,nz-1))
462  s_f(i,j,nz) = b1(i) * (h_tr*s_in(i,j,nz) + ent(i,nz)*s_f(i,j,nz-1))
463  enddo
464  do k=nz-1,1,-1 ; do i=is,ie
465  t_f(i,j,k) = t_f(i,j,k) + c1(i,k+1)*t_f(i,j,k+1)
466  s_f(i,j,k) = s_f(i,j,k) + c1(i,k+1)*s_f(i,j,k+1)
467  enddo ; enddo
468  enddo
469  endif
470