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)
 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 
)

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

Definition at line 28 of file MOM_isopycnal_slopes.F90.

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

References vert_fill_ts().

Referenced by mom_lateral_mixing_coeffs::calc_slope_functions().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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 334 of file MOM_isopycnal_slopes.F90.

334  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
335  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
336  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
337  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: T_in !< Input temperature [degC]
338  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: S_in !< Input salinity [ppt]
339  real, intent(in) :: kappa_dt !< A vertical diffusivity to use for smoothing
340  !! times a smoothing timescale [Z2 ~> m2].
341  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out) :: T_f !< Filled temperature [degC]
342  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out) :: S_f !< Filled salinity [ppt]
343  integer, optional, intent(in) :: halo_here !< Number of halo points to work on,
344  !! 0 by default
345  logical, optional, intent(in) :: larger_h_denom !< Present and true, add a large
346  !! enough minimal thickness in the denominator of
347  !! the flux calculations so that the fluxes are
348  !! never so large as eliminate the transmission
349  !! of information across groups of massless layers.
350  ! Local variables
351  real :: ent(SZI_(G),SZK_(G)+1) ! The diffusive entrainment (kappa*dt)/dz
352  ! between layers in a timestep [H ~> m or kg m-2].
353  real :: b1(SZI_(G)), d1(SZI_(G)) ! b1, c1, and d1 are variables used by the
354  real :: c1(SZI_(G),SZK_(G)) ! tridiagonal solver.
355  real :: kap_dt_x2 ! The 2*kappa_dt converted to H units [H2 ~> m2 or kg2 m-4].
356  real :: h_neglect ! A negligible thickness [H ~> m or kg m-2], to allow for zero thicknesses.
357  real :: h0 ! A negligible thickness to allow for zero thickness layers without
358  ! completely decouping groups of layers [H ~> m or kg m-2].
359  ! Often 0 < h_neglect << h0.
360  real :: h_tr ! h_tr is h at tracer points with a tiny thickness
361  ! added to ensure positive definiteness [H ~> m or kg m-2].
362  integer :: i, j, k, is, ie, js, je, nz, halo
363 
364  halo=0 ; if (present(halo_here)) halo = max(halo_here,0)
365 
366  is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo ; nz = gv%ke
367 
368  h_neglect = gv%H_subroundoff
369  kap_dt_x2 = (2.0*kappa_dt)*gv%Z_to_H**2
370  h0 = h_neglect
371  if (present(larger_h_denom)) then
372  if (larger_h_denom) h0 = 1.0e-16*sqrt(kappa_dt)*gv%Z_to_H
373  endif
374 
375  if (kap_dt_x2 <= 0.0) then
376  !$OMP parallel do default(shared)
377  do k=1,nz ; do j=js,je ; do i=is,ie
378  t_f(i,j,k) = t_in(i,j,k) ; s_f(i,j,k) = s_in(i,j,k)
379  enddo ; enddo ; enddo
380  else
381  !$OMP parallel do default(shared) private(ent,b1,d1,c1,h_tr)
382  do j=js,je
383  do i=is,ie
384  ent(i,2) = kap_dt_x2 / ((h(i,j,1)+h(i,j,2)) + h0)
385  h_tr = h(i,j,1) + h_neglect
386  b1(i) = 1.0 / (h_tr + ent(i,2))
387  d1(i) = b1(i) * h_tr
388  t_f(i,j,1) = (b1(i)*h_tr)*t_in(i,j,1)
389  s_f(i,j,1) = (b1(i)*h_tr)*s_in(i,j,1)
390  enddo
391  do k=2,nz-1 ; do i=is,ie
392  ent(i,k+1) = kap_dt_x2 / ((h(i,j,k)+h(i,j,k+1)) + h0)
393  h_tr = h(i,j,k) + h_neglect
394  c1(i,k) = ent(i,k) * b1(i)
395  b1(i) = 1.0 / ((h_tr + d1(i)*ent(i,k)) + ent(i,k+1))
396  d1(i) = b1(i) * (h_tr + d1(i)*ent(i,k))
397  t_f(i,j,k) = b1(i) * (h_tr*t_in(i,j,k) + ent(i,k)*t_f(i,j,k-1))
398  s_f(i,j,k) = b1(i) * (h_tr*s_in(i,j,k) + ent(i,k)*s_f(i,j,k-1))
399  enddo ; enddo
400  do i=is,ie
401  c1(i,nz) = ent(i,nz) * b1(i)
402  h_tr = h(i,j,nz) + h_neglect
403  b1(i) = 1.0 / (h_tr + d1(i)*ent(i,nz))
404  t_f(i,j,nz) = b1(i) * (h_tr*t_in(i,j,nz) + ent(i,nz)*t_f(i,j,nz-1))
405  s_f(i,j,nz) = b1(i) * (h_tr*s_in(i,j,nz) + ent(i,nz)*s_f(i,j,nz-1))
406  enddo
407  do k=nz-1,1,-1 ; do i=is,ie
408  t_f(i,j,k) = t_f(i,j,k) + c1(i,k+1)*t_f(i,j,k+1)
409  s_f(i,j,k) = s_f(i,j,k) + c1(i,k+1)*s_f(i,j,k+1)
410  enddo ; enddo
411  enddo
412  endif
413 

Referenced by calc_isoneutral_slopes(), mom_set_diffusivity::set_diffusivity(), mom_int_tide_input::set_int_tide_input(), and mom_thickness_diffuse::thickness_diffuse_full().

Here is the caller graph for this function: