MOM6
Functions/Subroutines
mom_full_convection Module Reference

Detailed Description

Does full convective adjustment of unstable regions via a strong diffusivity.

Functions/Subroutines

subroutine, public full_convection (G, GV, US, h, tv, T_adj, S_adj, p_surf, Kddt_smooth, Kddt_convect, halo)
 Calculate new temperatures and salinities that have been subject to full convective mixing. More...
 
logical function is_unstable (dRho_dT, dRho_dS, h_a, h_b, mix_A, mix_B, T_a, T_b, S_a, S_b, Te_aa, Te_bb, Se_aa, Se_bb, d_A, d_B)
 This function returns True if the profiles around the given interface will be statically unstable after mixing above below. The arguments are the ocean state above and below, including partial calculations from a tridiagonal solver. More...
 
subroutine smoothed_drdt_drds (h, tv, Kddt, dR_dT, dR_dS, G, GV, US, j, p_surf, halo)
 Returns the partial derivatives of locally referenced potential density with temperature and salinity after the properties have been smoothed with a small constant diffusivity. More...
 

Function/Subroutine Documentation

◆ full_convection()

subroutine, public mom_full_convection::full_convection ( 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,
type(thermo_var_ptrs), intent(in)  tv,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(out)  T_adj,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(out)  S_adj,
real, dimension(:,:), pointer  p_surf,
real, intent(in)  Kddt_smooth,
real, intent(in), optional  Kddt_convect,
integer, intent(in), optional  halo 
)

Calculate new temperatures and salinities that have been subject to full convective mixing.

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]t_adjAdjusted potential temperature [degC].
[out]s_adjAdjusted salinity [ppt].
p_surfThe pressure at the ocean surface [R L2 T-2 ~> Pa] (or NULL).
[in]kddt_smoothA smoothing vertical diffusivity times a timestep [H2 ~> m2 or kg2 m-4].
[in]kddt_convectA large convecting vertical diffusivity times a timestep [H2 ~> m2 or kg2 m-4].
[in]haloHalo width over which to compute

Definition at line 23 of file MOM_full_convection.F90.

23  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
24  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
25  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
26  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
27  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
28  type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
29  !! thermodynamic variables
30  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
31  intent(out) :: T_adj !< Adjusted potential temperature [degC].
32  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
33  intent(out) :: S_adj !< Adjusted salinity [ppt].
34  real, dimension(:,:), pointer :: p_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa] (or NULL).
35  real, intent(in) :: Kddt_smooth !< A smoothing vertical
36  !! diffusivity times a timestep [H2 ~> m2 or kg2 m-4].
37  real, optional, intent(in) :: Kddt_convect !< A large convecting vertical
38  !! diffusivity times a timestep [H2 ~> m2 or kg2 m-4].
39  integer, optional, intent(in) :: halo !< Halo width over which to compute
40 
41  ! Local variables
42  real, dimension(SZI_(G),SZK_(G)+1) :: &
43  dRho_dT, & ! The derivative of density with temperature [R degC-1 ~> kg m-3 degC-1]
44  dRho_dS ! The derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1].
45  real :: h_neglect, h0 ! A thickness that is so small it is usually lost
46  ! in roundoff and can be neglected [H ~> m or kg m-2].
47 ! logical :: use_EOS ! If true, density is calculated from T & S using an equation of state.
48  real, dimension(SZI_(G),SZK0_(G)) :: &
49  Te_a, & ! A partially updated temperature estimate including the influnce from
50  ! mixing with layers above rescaled by a factor of d_a [degC].
51  ! This array is discreted on tracer cells, but contains an extra
52  ! layer at the top for algorithmic convenience.
53  se_a ! A partially updated salinity estimate including the influnce from
54  ! mixing with layers above rescaled by a factor of d_a [ppt].
55  ! This array is discreted on tracer cells, but contains an extra
56  ! layer at the top for algorithmic convenience.
57  real, dimension(SZI_(G),SZK_(G)+1) :: &
58  Te_b, & ! A partially updated temperature estimate including the influnce from
59  ! mixing with layers below rescaled by a factor of d_b [degC].
60  ! This array is discreted on tracer cells, but contains an extra
61  ! layer at the bottom for algorithmic convenience.
62  se_b ! A partially updated salinity estimate including the influnce from
63  ! mixing with layers below rescaled by a factor of d_b [ppt].
64  ! This array is discreted on tracer cells, but contains an extra
65  ! layer at the bottom for algorithmic convenience.
66  real, dimension(SZI_(G),SZK_(G)+1) :: &
67  c_a, & ! The fractional influence of the properties of the layer below
68  ! in the final properies with a downward-first solver, nondim.
69  d_a, & ! The fractional influence of the properties of the layer in question
70  ! and layers above in the final properies with a downward-first solver, nondim.
71  ! d_a = 1.0 - c_a
72  c_b, & ! The fractional influence of the properties of the layer above
73  ! in the final properies with a upward-first solver, nondim.
74  d_b ! The fractional influence of the properties of the layer in question
75  ! and layers below in the final properies with a upward-first solver, nondim.
76  ! d_b = 1.0 - c_b
77  real, dimension(SZI_(G),SZK_(G)+1) :: &
78  mix !< The amount of mixing across the interface between layers [H ~> m or kg m-2].
79  real :: mix_len ! The length-scale of mixing, when it is active [H ~> m or kg m-2]
80  real :: h_b, h_a ! The thicknessses of the layers above and below an interface [H ~> m or kg m-2]
81  real :: b_b, b_a ! Inverse pivots used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1].
82 
83  real :: kap_dt_x2 ! The product of 2*kappa*dt [H2 ~> m2 or kg2 m-4].
84 
85  logical, dimension(SZI_(G)) :: do_i ! Do more work on this column.
86  logical, dimension(SZI_(G)) :: last_down ! The last setup pass was downward.
87  integer, dimension(SZI_(G)) :: change_ct ! The number of interfaces where the
88  ! mixing has changed this iteration.
89  integer :: changed_col ! The number of colums whose mixing changed.
90  integer :: i, j, k, is, ie, js, je, nz, itt
91 
92  if (present(halo)) then
93  is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
94  else
95  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
96  endif
97  nz = g%ke
98 
99  if (.not.associated(tv%eqn_of_state)) return
100 
101  h_neglect = gv%H_subroundoff
102  kap_dt_x2 = 0.0
103  if (present(kddt_convect)) kap_dt_x2 = 2.0*kddt_convect
104  mix_len = (1.0e20 * nz) * (g%max_depth * gv%Z_to_H)
105  h0 = 1.0e-16*sqrt(kddt_smooth) + h_neglect
106 
107  do j=js,je
108  mix(:,:) = 0.0 ; d_b(:,:) = 1.0
109  ! These would be Te_b(:,:) = tv%T(:,j,:), etc., but the values are not used
110  te_b(:,:) = 0.0 ; se_b(:,:) = 0.0
111 
112  call smoothed_drdt_drds(h, tv, kddt_smooth, drho_dt, drho_ds, g, gv, us, j, p_surf, halo)
113 
114  do i=is,ie
115  do_i(i) = (g%mask2dT(i,j) > 0.0)
116 
117  d_a(i,1) = 1.0
118  last_down(i) = .true. ! This is set for debuggers.
119  ! These are extra values are used for convenience in the stability test
120  te_a(i,0) = 0.0 ; se_a(i,0) = 0.0
121  enddo
122 
123  do itt=1,nz ! At least 2 interfaces will change with each full pass, or the
124  ! iterations stop, so the maximum count of nz is very conservative.
125 
126  do i=is,ie ; change_ct(i) = 0 ; enddo
127  ! Move down the water column, finding unstable interfaces, and building up the
128  ! temporary arrays for the tridiagonal solver.
129  do k=2,nz ; do i=is,ie ; if (do_i(i)) then
130 
131  h_a = h(i,j,k-1) + h_neglect ; h_b = h(i,j,k) + h_neglect
132  if (mix(i,k) <= 0.0) then
133  if (is_unstable(drho_dt(i,k), drho_ds(i,k), h_a, h_b, mix(i,k-1), mix(i,k+1), &
134  tv%T(i,j,k-1), tv%T(i,j,k), tv%S(i,j,k-1), tv%S(i,j,k), &
135  te_a(i,k-2), te_b(i,k+1), se_a(i,k-2), se_b(i,k+1), &
136  d_a(i,k-1), d_b(i,k+1))) then
137  mix(i,k) = mix_len
138  if (kap_dt_x2 > 0.0) mix(i,k) = kap_dt_x2 / ((h(i,j,k-1)+h(i,j,k)) + h0)
139  change_ct(i) = change_ct(i) + 1
140  endif
141  endif
142 
143  b_a = 1.0 / ((h_a + d_a(i,k-1)*mix(i,k-1)) + mix(i,k))
144  if (mix(i,k) <= 0.0) then
145  c_a(i,k) = 0.0 ; d_a(i,k) = 1.0
146  else
147  d_a(i,k) = b_a * (h_a + d_a(i,k-1)*mix(i,k-1)) ! = 1.0-c_a(i,K)
148  c_a(i,k) = 1.0 ; if (d_a(i,k) > epsilon(b_a)) c_a(i,k) = b_a * mix(i,k)
149  endif
150 
151  if (k>2) then
152  te_a(i,k-1) = b_a * (h_a*tv%T(i,j,k-1) + mix(i,k-1)*te_a(i,k-2))
153  se_a(i,k-1) = b_a * (h_a*tv%S(i,j,k-1) + mix(i,k-1)*se_a(i,k-2))
154  else
155  te_a(i,k-1) = b_a * (h_a*tv%T(i,j,k-1))
156  se_a(i,k-1) = b_a * (h_a*tv%S(i,j,k-1))
157  endif
158  endif ; enddo ; enddo
159 
160  ! Determine which columns might have further instabilities.
161  changed_col = 0
162  do i=is,ie ; if (do_i(i)) then
163  if (change_ct(i) == 0) then
164  last_down(i) = .true. ; do_i(i) = .false.
165  else
166  changed_col = changed_col + 1 ; change_ct(i) = 0
167  endif
168  endif ; enddo
169  if (changed_col == 0) exit ! No more columns are unstable.
170 
171  ! This is the same as above, but with the direction reversed (bottom to top)
172  do k=nz,2,-1 ; do i=is,ie ; if (do_i(i)) then
173 
174  h_a = h(i,j,k-1) + h_neglect ; h_b = h(i,j,k) + h_neglect
175  if (mix(i,k) <= 0.0) then
176  if (is_unstable(drho_dt(i,k), drho_ds(i,k), h_a, h_b, mix(i,k-1), mix(i,k+1), &
177  tv%T(i,j,k-1), tv%T(i,j,k), tv%S(i,j,k-1), tv%S(i,j,k), &
178  te_a(i,k-2), te_b(i,k+1), se_a(i,k-2), se_b(i,k+1), &
179  d_a(i,k-1), d_b(i,k+1))) then
180  mix(i,k) = mix_len
181  if (kap_dt_x2 > 0.0) mix(i,k) = kap_dt_x2 / ((h(i,j,k-1)+h(i,j,k)) + h0)
182  change_ct(i) = change_ct(i) + 1
183  endif
184  endif
185 
186  b_b = 1.0 / ((h_b + d_b(i,k+1)*mix(i,k+1)) + mix(i,k))
187  if (mix(i,k) <= 0.0) then
188  c_b(i,k) = 0.0 ; d_b(i,k) = 1.0
189  else
190  d_b(i,k) = b_b * (h_b + d_b(i,k+1)*mix(i,k+1)) ! = 1.0-c_b(i,K)
191  c_b(i,k) = 1.0 ; if (d_b(i,k) > epsilon(b_b)) c_b(i,k) = b_b * mix(i,k)
192  endif
193 
194  if (k<nz) then
195  te_b(i,k) = b_b * (h_b*tv%T(i,j,k) + mix(i,k+1)*te_b(i,k+1))
196  se_b(i,k) = b_b * (h_b*tv%S(i,j,k) + mix(i,k+1)*se_b(i,k+1))
197  else
198  te_b(i,k) = b_b * (h_b*tv%T(i,j,k))
199  se_b(i,k) = b_b * (h_b*tv%S(i,j,k))
200  endif
201  endif ; enddo ; enddo
202 
203  ! Determine which columns might have further instabilities.
204  changed_col = 0
205  do i=is,ie ; if (do_i(i)) then
206  if (change_ct(i) == 0) then
207  last_down(i) = .false. ; do_i(i) = .false.
208  else
209  changed_col = changed_col + 1 ; change_ct(i) = 0
210  endif
211  endif ; enddo
212  if (changed_col == 0) exit ! No more columns are unstable.
213 
214  enddo ! End of iterations, all columns are now stable.
215 
216  ! Do the final return pass on the columns where the penultimate pass was downward.
217  do i=is,ie ; do_i(i) = ((g%mask2dT(i,j) > 0.0) .and. last_down(i)) ; enddo
218  do i=is,ie ; if (do_i(i)) then
219  h_a = h(i,j,nz) + h_neglect
220  b_a = 1.0 / (h_a + d_a(i,nz)*mix(i,nz))
221  t_adj(i,j,nz) = b_a * (h_a*tv%T(i,j,nz) + mix(i,nz)*te_a(i,nz-1))
222  s_adj(i,j,nz) = b_a * (h_a*tv%S(i,j,nz) + mix(i,nz)*se_a(i,nz-1))
223  endif ; enddo
224  do k=nz-1,1,-1 ; do i=is,ie ; if (do_i(i)) then
225  t_adj(i,j,k) = te_a(i,k) + c_a(i,k+1)*t_adj(i,j,k+1)
226  s_adj(i,j,k) = se_a(i,k) + c_a(i,k+1)*s_adj(i,j,k+1)
227  endif ; enddo ; enddo
228 
229  do i=is,ie ; if (do_i(i)) then
230  k = 1 ! A hook for debugging.
231  endif ; enddo
232 
233  ! Do the final return pass on the columns where the penultimate pass was upward.
234  ! Also do a simple copy of T & S values on land points.
235  do i=is,ie
236  do_i(i) = ((g%mask2dT(i,j) > 0.0) .and. .not.last_down(i))
237  if (do_i(i)) then
238  h_b = h(i,j,1) + h_neglect
239  b_b = 1.0 / (h_b + d_b(i,2)*mix(i,2))
240  t_adj(i,j,1) = b_b * (h_b*tv%T(i,j,1) + mix(i,2)*te_b(i,2))
241  s_adj(i,j,1) = b_b * (h_b*tv%S(i,j,1) + mix(i,2)*se_b(i,2))
242  elseif (g%mask2dT(i,j) <= 0.0) then
243  t_adj(i,j,1) = tv%T(i,j,1) ; s_adj(i,j,1) = tv%S(i,j,1)
244  endif
245  enddo
246  do k=2,nz ; do i=is,ie
247  if (do_i(i)) then
248  t_adj(i,j,k) = te_b(i,k) + c_b(i,k)*t_adj(i,j,k-1)
249  s_adj(i,j,k) = se_b(i,k) + c_b(i,k)*s_adj(i,j,k-1)
250  elseif (g%mask2dT(i,j) <= 0.0) then
251  t_adj(i,j,k) = tv%T(i,j,k) ; s_adj(i,j,k) = tv%S(i,j,k)
252  endif
253  enddo ; enddo
254 
255  do i=is,ie ; if (do_i(i)) then
256  k = 1 ! A hook for debugging.
257  endif ; enddo
258 
259  enddo ! j-loop
260 
261  k = 1 ! A hook for debugging.
262 
263  ! The following set of expressions for the final values are derived from the the partial
264  ! updates for the estimated temperatures and salinities around an interface, then directly
265  ! solving for the final temperatures and salinities. They are here for later reference
266  ! and to document an intermediate step in the stability calculation.
267  ! hp_a = (h_a + d_a(i,K-1)*mix(i,K-1))
268  ! hp_b = (h_b + d_b(i,K+1)*mix(i,K+1))
269  ! b2_c = 1.0 / (hp_a*hp_b + (hp_a + hp_b) * mix(i,K))
270  ! Th_a = h_a*tv%T(i,j,k-1) + mix(i,K-1)*Te_a(i,k-2)
271  ! Th_b = h_b*tv%T(i,j,k) + mix(i,K+1)*Te_b(i,k+1)
272  ! T_fin(i,k) = ( (hp_a + mix(i,K)) * Th_b + Th_a * mix(i,K) ) * b2_c
273  ! T_fin(i,k-1) = ( (hp_b + mix(i,K)) * Th_a + Th_b * mix(i,K) ) * b2_c
274  ! Sh_a = h_a*tv%S(i,j,k-1) + mix(i,K-1)*Se_a(i,k-2)
275  ! Sh_b = h_b*tv%S(i,j,k) + mix(i,K+1)*Se_b(i,k+1)
276  ! S_fin(i,k) = ( (hp_a + mix(i,K)) * Sh_b + Sh_a * mix(i,K) ) * b2_c
277  ! S_fin(i,k-1) = ( (hp_b + mix(i,K)) * Sh_a + Sh_b * mix(i,K) ) * b2_c
278 

◆ is_unstable()

logical function mom_full_convection::is_unstable ( real, intent(in)  dRho_dT,
real, intent(in)  dRho_dS,
real, intent(in)  h_a,
real, intent(in)  h_b,
real, intent(in)  mix_A,
real, intent(in)  mix_B,
real, intent(in)  T_a,
real, intent(in)  T_b,
real, intent(in)  S_a,
real, intent(in)  S_b,
real, intent(in)  Te_aa,
real, intent(in)  Te_bb,
real, intent(in)  Se_aa,
real, intent(in)  Se_bb,
real, intent(in)  d_A,
real, intent(in)  d_B 
)
private

This function returns True if the profiles around the given interface will be statically unstable after mixing above below. The arguments are the ocean state above and below, including partial calculations from a tridiagonal solver.

Parameters
[in]drho_dtThe derivative of in situ density with temperature [R degC-1 ~> kg m-3 degC-1]
[in]drho_dsThe derivative of in situ density with salinity [R ppt-1 ~> kg m-3 ppt-1]
[in]h_aThe thickness of the layer above [H ~> m or kg m-2]
[in]h_bThe thickness of the layer below [H ~> m or kg m-2]
[in]mix_aThe time integrated mixing rate of the interface above [H ~> m or kg m-2]
[in]mix_bThe time integrated mixing rate of the interface below [H ~> m or kg m-2]
[in]t_aThe initial temperature of the layer above [degC]
[in]t_bThe initial temperature of the layer below [degC]
[in]s_aThe initial salinity of the layer below [ppt]
[in]s_bThe initial salinity of the layer below [ppt]
[in]te_aaThe estimated temperature two layers above rescaled by d_A [degC]
[in]te_bbThe estimated temperature two layers below rescaled by d_B [degC]
[in]se_aaThe estimated salinity two layers above rescaled by d_A [ppt]
[in]se_bbThe estimated salinity two layers below rescaled by d_B [ppt]
[in]d_aThe rescaling dependency across the interface above, nondim.
[in]d_bThe rescaling dependency across the interface below, nondim.
Returns
The return value, true if the profile is statically unstable around the interface in question.

Definition at line 286 of file MOM_full_convection.F90.

286  real, intent(in) :: dRho_dT !< The derivative of in situ density with temperature [R degC-1 ~> kg m-3 degC-1]
287  real, intent(in) :: dRho_dS !< The derivative of in situ density with salinity [R ppt-1 ~> kg m-3 ppt-1]
288  real, intent(in) :: h_a !< The thickness of the layer above [H ~> m or kg m-2]
289  real, intent(in) :: h_b !< The thickness of the layer below [H ~> m or kg m-2]
290  real, intent(in) :: mix_A !< The time integrated mixing rate of the interface above [H ~> m or kg m-2]
291  real, intent(in) :: mix_B !< The time integrated mixing rate of the interface below [H ~> m or kg m-2]
292  real, intent(in) :: T_a !< The initial temperature of the layer above [degC]
293  real, intent(in) :: T_b !< The initial temperature of the layer below [degC]
294  real, intent(in) :: S_a !< The initial salinity of the layer below [ppt]
295  real, intent(in) :: S_b !< The initial salinity of the layer below [ppt]
296  real, intent(in) :: Te_aa !< The estimated temperature two layers above rescaled by d_A [degC]
297  real, intent(in) :: Te_bb !< The estimated temperature two layers below rescaled by d_B [degC]
298  real, intent(in) :: Se_aa !< The estimated salinity two layers above rescaled by d_A [ppt]
299  real, intent(in) :: Se_bb !< The estimated salinity two layers below rescaled by d_B [ppt]
300  real, intent(in) :: d_A !< The rescaling dependency across the interface above, nondim.
301  real, intent(in) :: d_B !< The rescaling dependency across the interface below, nondim.
302  logical :: is_unstable !< The return value, true if the profile is statically unstable
303  !! around the interface in question.
304 
305  ! These expressions for the local stability are long, but they have been carefully
306  ! grouped for accuracy even when the mixing rates are huge or tiny, and common
307  ! positive definite factors that would appear in the final expression for the
308  ! locally referenced potential density difference across an interface have been omitted.
309  is_unstable = (drho_dt * ((h_a * h_b * (t_b - t_a) + &
310  mix_a*mix_b * (d_a*te_bb - d_b*te_aa)) + &
311  (h_a*mix_b * (te_bb - d_b*t_a) + &
312  h_b*mix_a * (d_a*t_b - te_aa)) ) + &
313  drho_ds * ((h_a * h_b * (s_b - s_a) + &
314  mix_a*mix_b * (d_a*se_bb - d_b*se_aa)) + &
315  (h_a*mix_b * (se_bb - d_b*s_a) + &
316  h_b*mix_a * (d_a*s_b - se_aa)) ) < 0.0)

◆ smoothed_drdt_drds()

subroutine mom_full_convection::smoothed_drdt_drds ( real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h,
type(thermo_var_ptrs), intent(in)  tv,
real, intent(in)  Kddt,
real, dimension( g %isd: g %ied, g %ke+1), intent(out)  dR_dT,
real, dimension( g %isd: g %ied, g %ke+1), intent(out)  dR_dS,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
integer, intent(in)  j,
real, dimension(:,:), pointer  p_surf,
integer, intent(in), optional  halo 
)
private

Returns the partial derivatives of locally referenced potential density with temperature and salinity after the properties have been smoothed with a small 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]tvA structure pointing to various thermodynamic variables
[in]kddtA diffusivity times a time increment [H2 ~> m2 or kg2 m-4].
[out]dr_dtDerivative of locally referenced
[out]dr_dsDerivative of locally referenced
[in]usA dimensional unit scaling type
[in]jThe j-point to work on.
p_surfThe pressure at the ocean surface [R L2 T-2 ~> Pa].
[in]haloHalo width over which to compute

Definition at line 323 of file MOM_full_convection.F90.

323  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
324  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
325  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
326  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
327  type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
328  !! thermodynamic variables
329  real, intent(in) :: Kddt !< A diffusivity times a time increment [H2 ~> m2 or kg2 m-4].
330  real, dimension(SZI_(G),SZK_(G)+1), &
331  intent(out) :: dR_dT !< Derivative of locally referenced
332  !! potential density with temperature [R degC-1 ~> kg m-3 degC-1]
333  real, dimension(SZI_(G),SZK_(G)+1), &
334  intent(out) :: dR_dS !< Derivative of locally referenced
335  !! potential density with salinity [R degC-1 ~> kg m-3 ppt-1]
336  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
337  integer, intent(in) :: j !< The j-point to work on.
338  real, dimension(:,:), pointer :: p_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa].
339  integer, optional, intent(in) :: halo !< Halo width over which to compute
340 
341  ! Local variables
342  real :: mix(SZI_(G),SZK_(G)+1) ! The diffusive mixing length (kappa*dt)/dz
343  ! between layers within in a timestep [H ~> m or kg m-2].
344  real :: b1(SZI_(G)), d1(SZI_(G)) ! b1, c1, and d1 are variables used by the
345  real :: c1(SZI_(G),SZK_(G)) ! tridiagonal solver.
346  real :: T_f(SZI_(G),SZK_(G)) ! Filtered temperatures [degC]
347  real :: S_f(SZI_(G),SZK_(G)) ! Filtered salinities [ppt]
348  real :: pres(SZI_(G)) ! Interface pressures [R L2 T-2 ~> Pa].
349  real :: T_EOS(SZI_(G)) ! Filtered and vertically averaged temperatures [degC]
350  real :: S_EOS(SZI_(G)) ! Filtered and vertically averaged salinities [ppt]
351  real :: kap_dt_x2 ! The product of 2*kappa*dt [H2 ~> m2 or kg2 m-4].
352  real :: h_neglect, h0 ! Negligible thicknesses to allow for zero thicknesses,
353  ! [H ~> m or kg m-2].
354  real :: h_tr ! The thickness at tracer points, plus h_neglect [H ~> m or kg m-2].
355  integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
356  integer :: i, k, is, ie, nz
357 
358  if (present(halo)) then
359  is = g%isc-halo ; ie = g%iec+halo
360  else
361  is = g%isc ; ie = g%iec
362  endif
363  nz = g%ke
364 
365  h_neglect = gv%H_subroundoff
366  kap_dt_x2 = 2.0*kddt
367 
368  if (kddt <= 0.0) then
369  do k=1,nz ; do i=is,ie
370  t_f(i,k) = tv%T(i,j,k) ; s_f(i,k) = tv%S(i,j,k)
371  enddo ; enddo
372  else
373  h0 = 1.0e-16*sqrt(kddt) + h_neglect
374  do i=is,ie
375  mix(i,2) = kap_dt_x2 / ((h(i,j,1)+h(i,j,2)) + h0)
376 
377  h_tr = h(i,j,1) + h_neglect
378  b1(i) = 1.0 / (h_tr + mix(i,2))
379  d1(i) = b1(i) * h(i,j,1)
380  t_f(i,1) = (b1(i)*h_tr)*tv%T(i,j,1)
381  s_f(i,1) = (b1(i)*h_tr)*tv%S(i,j,1)
382  enddo
383  do k=2,nz-1 ; do i=is,ie
384  mix(i,k+1) = kap_dt_x2 / ((h(i,j,k)+h(i,j,k+1)) + h0)
385 
386  c1(i,k) = mix(i,k) * b1(i)
387  h_tr = h(i,j,k) + h_neglect
388  b1(i) = 1.0 / ((h_tr + d1(i)*mix(i,k)) + mix(i,k+1))
389  d1(i) = b1(i) * (h_tr + d1(i)*mix(i,k))
390  t_f(i,k) = b1(i) * (h_tr*tv%T(i,j,k) + mix(i,k)*t_f(i,k-1))
391  s_f(i,k) = b1(i) * (h_tr*tv%S(i,j,k) + mix(i,k)*s_f(i,k-1))
392  enddo ; enddo
393  do i=is,ie
394  c1(i,nz) = mix(i,nz) * b1(i)
395  h_tr = h(i,j,nz) + h_neglect
396  b1(i) = 1.0 / (h_tr + d1(i)*mix(i,nz))
397  t_f(i,nz) = b1(i) * (h_tr*tv%T(i,j,nz) + mix(i,nz)*t_f(i,nz-1))
398  s_f(i,nz) = b1(i) * (h_tr*tv%S(i,j,nz) + mix(i,nz)*s_f(i,nz-1))
399  enddo
400  do k=nz-1,1,-1 ; do i=is,ie
401  t_f(i,k) = t_f(i,k) + c1(i,k+1)*t_f(i,k+1)
402  s_f(i,k) = s_f(i,k) + c1(i,k+1)*s_f(i,k+1)
403  enddo ; enddo
404  endif
405 
406  if (associated(p_surf)) then
407  do i=is,ie ; pres(i) = p_surf(i,j) ; enddo
408  else
409  do i=is,ie ; pres(i) = 0.0 ; enddo
410  endif
411  eosdom(:) = eos_domain(g%HI, halo)
412  call calculate_density_derivs(t_f(:,1), s_f(:,1), pres, dr_dt(:,1), dr_ds(:,1), tv%eqn_of_state, eosdom)
413  do i=is,ie ; pres(i) = pres(i) + h(i,j,1)*(gv%H_to_RZ*gv%g_Earth) ; enddo
414  do k=2,nz
415  do i=is,ie
416  t_eos(i) = 0.5*(t_f(i,k-1) + t_f(i,k))
417  s_eos(i) = 0.5*(s_f(i,k-1) + s_f(i,k))
418  enddo
419  call calculate_density_derivs(t_eos, s_eos, pres, dr_dt(:,k), dr_ds(:,k), tv%eqn_of_state, eosdom)
420  do i=is,ie ; pres(i) = pres(i) + h(i,j,k)*(gv%H_to_RZ*gv%g_Earth) ; enddo
421  enddo
422  call calculate_density_derivs(t_f(:,nz), s_f(:,nz), pres, dr_dt(:,nz+1), dr_ds(:,nz+1), &
423  tv%eqn_of_state, eosdom)
424