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, 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, 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,
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]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 [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 22 of file MOM_full_convection.F90.

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

References is_unstable(), and smoothed_drdt_drds().

Referenced by mom_set_diffusivity::set_diffusivity().

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

◆ 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 [kg m-3 degC-1]
[in]drho_dsThe derivative of in situ density with salinity [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 284 of file MOM_full_convection.F90.

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

Referenced by full_convection().

Here is the caller graph for this function:

◆ 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,
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]jThe j-point to work on.
p_surfThe pressure at the ocean surface [Pa].
[in]haloHalo width over which to compute

Definition at line 321 of file MOM_full_convection.F90.

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

Referenced by full_convection().

Here is the caller graph for this function: