MOM_full_convection.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Does full convective adjustment of unstable regions via a strong diffusivity.
7
9use mom_interface_heights, only : thickness_to_dz
13use mom_eos, only : calculate_density_derivs, eos_domain
14
15implicit none ; private
16
17#include <MOM_memory.h>
18
19public full_convection
20
21contains
22
23!> Calculate new temperatures and salinities that have been subject to full convective mixing.
24subroutine full_convection(G, GV, US, h, tv, T_adj, S_adj, p_surf, Kddt_smooth, halo)
25 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
26 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
27 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
28 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
29 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
30 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
31 !! thermodynamic variables
32 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
33 intent(out) :: t_adj !< Adjusted potential temperature [C ~> degC].
34 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
35 intent(out) :: s_adj !< Adjusted salinity [S ~> ppt].
36 real, dimension(:,:), pointer :: p_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa] (or NULL).
37 real, intent(in) :: kddt_smooth !< A smoothing vertical diffusivity
38 !! times a timestep [H Z ~> m2 or kg m-1].
39 integer, intent(in) :: halo !< Halo width over which to compute
40
41 ! Local variables
42 real, dimension(SZI_(G),SZK_(GV)+1) :: &
43 drho_dt, & ! The derivative of density with temperature [R C-1 ~> kg m-3 degC-1]
44 drho_ds ! The derivative of density with salinity [R S-1 ~> kg m-3 ppt-1].
45 real :: dz(szi_(g),szk_(gv)) ! Height change across layers [Z ~> m]
46 real :: h_neglect ! A thickness that is so small it is usually lost
47 ! in roundoff and can be neglected [H ~> m or kg m-2].
48! logical :: use_EOS ! If true, density is calculated from T & S using an equation of state.
49 real, dimension(SZI_(G),SZK0_(G)) :: &
50 te_a, & ! A partially updated temperature estimate including the influence from
51 ! mixing with layers above rescaled by a factor of d_a [C ~> degC].
52 ! This array is discretized on tracer cells, but contains an extra
53 ! layer at the top for algorithmic convenience.
54 se_a ! A partially updated salinity estimate including the influence from
55 ! mixing with layers above rescaled by a factor of d_a [S ~> ppt].
56 ! This array is discretized on tracer cells, but contains an extra
57 ! layer at the top for algorithmic convenience.
58 real, dimension(SZI_(G),SZK_(GV)+1) :: &
59 te_b, & ! A partially updated temperature estimate including the influence from
60 ! mixing with layers below rescaled by a factor of d_b [C ~> degC].
61 ! This array is discretized on tracer cells, but contains an extra
62 ! layer at the bottom for algorithmic convenience.
63 se_b ! A partially updated salinity estimate including the influence from
64 ! mixing with layers below rescaled by a factor of d_b [S ~> ppt].
65 ! This array is discretized on tracer cells, but contains an extra
66 ! layer at the bottom for algorithmic convenience.
67 real, dimension(SZI_(G),SZK_(GV)+1) :: &
68 c_a, & ! The fractional influence of the properties of the layer below
69 ! in the final properties with a downward-first solver [nondim]
70 d_a, & ! The fractional influence of the properties of the layer in question
71 ! and layers above in the final properties with a downward-first solver [nondim]
72 ! d_a = 1.0 - c_a
73 c_b, & ! The fractional influence of the properties of the layer above
74 ! in the final properties with a upward-first solver [nondim]
75 d_b ! The fractional influence of the properties of the layer in question
76 ! and layers below in the final properties with a upward-first solver [nondim]
77 ! d_b = 1.0 - c_b
78 real, dimension(SZI_(G),SZK_(GV)+1) :: &
79 mix !< The amount of mixing across the interface between layers [H ~> m or kg m-2].
80 real :: mix_len ! The length-scale of mixing, when it is active [H ~> m or kg m-2]
81 real :: h_b, h_a ! The thicknesses of the layers above and below an interface [H ~> m or kg m-2]
82 real :: b_b, b_a ! Inverse pivots used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1].
83
84 logical, dimension(SZI_(G)) :: do_i ! Do more work on this column.
85 logical, dimension(SZI_(G)) :: last_down ! The last setup pass was downward.
86 integer, dimension(SZI_(G)) :: change_ct ! The number of interfaces where the
87 ! mixing has changed this iteration.
88 integer :: changed_col ! The number of columns whose mixing changed.
89 integer :: i, j, k, is, ie, js, je, nz, itt
90
91 is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
92 nz = gv%ke
93
94 if (.not.associated(tv%eqn_of_state)) return
95
96 h_neglect = gv%H_subroundoff
97 mix_len = (1.0e20 * nz) * (g%max_depth * us%Z_to_m * gv%m_to_H)
98
99 do j=js,je
100 mix(:,:) = 0.0 ; d_b(:,:) = 1.0
101 ! These would be Te_b(:,:) = tv%T(:,j,:), etc., but the values are not used
102 te_b(:,:) = 0.0 ; se_b(:,:) = 0.0
103
104 ! Find the vertical distances across layers.
105 call thickness_to_dz(h, tv, dz, j, g, gv, halo_size=halo)
106
107 call smoothed_drdt_drds(h, dz, tv, kddt_smooth, drho_dt, drho_ds, g, gv, us, j, p_surf, halo)
108
109 do i=is,ie
110 do_i(i) = (g%mask2dT(i,j) > 0.0)
111
112 d_a(i,1) = 1.0
113 last_down(i) = .true. ! This is set for debuggers.
114 ! These are extra values are used for convenience in the stability test
115 te_a(i,0) = 0.0 ; se_a(i,0) = 0.0
116 enddo
117
118 do itt=1,nz ! At least 2 interfaces will change with each full pass, or the
119 ! iterations stop, so the maximum count of nz is very conservative.
120
121 do i=is,ie ; change_ct(i) = 0 ; enddo
122 ! Move down the water column, finding unstable interfaces, and building up the
123 ! temporary arrays for the tridiagonal solver.
124 do k=2,nz ; do i=is,ie ; if (do_i(i)) then
125
126 h_a = h(i,j,k-1) + h_neglect ; h_b = h(i,j,k) + h_neglect
127 if (mix(i,k) <= 0.0) then
128 if (is_unstable(drho_dt(i,k), drho_ds(i,k), h_a, h_b, mix(i,k-1), mix(i,k+1), &
129 tv%T(i,j,k-1), tv%T(i,j,k), tv%S(i,j,k-1), tv%S(i,j,k), &
130 te_a(i,k-2), te_b(i,k+1), se_a(i,k-2), se_b(i,k+1), &
131 d_a(i,k-1), d_b(i,k+1))) then
132 mix(i,k) = mix_len
133 change_ct(i) = change_ct(i) + 1
134 endif
135 endif
136
137 b_a = 1.0 / ((h_a + d_a(i,k-1)*mix(i,k-1)) + mix(i,k))
138 if (mix(i,k) <= 0.0) then
139 c_a(i,k) = 0.0 ; d_a(i,k) = 1.0
140 else
141 d_a(i,k) = b_a * (h_a + d_a(i,k-1)*mix(i,k-1)) ! = 1.0-c_a(i,K)
142 c_a(i,k) = 1.0 ; if (d_a(i,k) > epsilon(b_a)) c_a(i,k) = b_a * mix(i,k)
143 endif
144
145 if (k>2) then
146 te_a(i,k-1) = b_a * (h_a*tv%T(i,j,k-1) + mix(i,k-1)*te_a(i,k-2))
147 se_a(i,k-1) = b_a * (h_a*tv%S(i,j,k-1) + mix(i,k-1)*se_a(i,k-2))
148 else
149 te_a(i,k-1) = b_a * (h_a*tv%T(i,j,k-1))
150 se_a(i,k-1) = b_a * (h_a*tv%S(i,j,k-1))
151 endif
152 endif ; enddo ; enddo
153
154 ! Determine which columns might have further instabilities.
155 changed_col = 0
156 do i=is,ie ; if (do_i(i)) then
157 if (change_ct(i) == 0) then
158 last_down(i) = .true. ; do_i(i) = .false.
159 else
160 changed_col = changed_col + 1 ; change_ct(i) = 0
161 endif
162 endif ; enddo
163 if (changed_col == 0) exit ! No more columns are unstable.
164
165 ! This is the same as above, but with the direction reversed (bottom to top)
166 do k=nz,2,-1 ; do i=is,ie ; if (do_i(i)) then
167
168 h_a = h(i,j,k-1) + h_neglect ; h_b = h(i,j,k) + h_neglect
169 if (mix(i,k) <= 0.0) then
170 if (is_unstable(drho_dt(i,k), drho_ds(i,k), h_a, h_b, mix(i,k-1), mix(i,k+1), &
171 tv%T(i,j,k-1), tv%T(i,j,k), tv%S(i,j,k-1), tv%S(i,j,k), &
172 te_a(i,k-2), te_b(i,k+1), se_a(i,k-2), se_b(i,k+1), &
173 d_a(i,k-1), d_b(i,k+1))) then
174 mix(i,k) = mix_len
175 change_ct(i) = change_ct(i) + 1
176 endif
177 endif
178
179 b_b = 1.0 / ((h_b + d_b(i,k+1)*mix(i,k+1)) + mix(i,k))
180 if (mix(i,k) <= 0.0) then
181 c_b(i,k) = 0.0 ; d_b(i,k) = 1.0
182 else
183 d_b(i,k) = b_b * (h_b + d_b(i,k+1)*mix(i,k+1)) ! = 1.0-c_b(i,K)
184 c_b(i,k) = 1.0 ; if (d_b(i,k) > epsilon(b_b)) c_b(i,k) = b_b * mix(i,k)
185 endif
186
187 if (k<nz) then
188 te_b(i,k) = b_b * (h_b*tv%T(i,j,k) + mix(i,k+1)*te_b(i,k+1))
189 se_b(i,k) = b_b * (h_b*tv%S(i,j,k) + mix(i,k+1)*se_b(i,k+1))
190 else
191 te_b(i,k) = b_b * (h_b*tv%T(i,j,k))
192 se_b(i,k) = b_b * (h_b*tv%S(i,j,k))
193 endif
194 endif ; enddo ; enddo
195
196 ! Determine which columns might have further instabilities.
197 changed_col = 0
198 do i=is,ie ; if (do_i(i)) then
199 if (change_ct(i) == 0) then
200 last_down(i) = .false. ; do_i(i) = .false.
201 else
202 changed_col = changed_col + 1 ; change_ct(i) = 0
203 endif
204 endif ; enddo
205 if (changed_col == 0) exit ! No more columns are unstable.
206
207 enddo ! End of iterations, all columns are now stable.
208
209 ! Do the final return pass on the columns where the penultimate pass was downward.
210 do i=is,ie ; do_i(i) = ((g%mask2dT(i,j) > 0.0) .and. last_down(i)) ; enddo
211 do i=is,ie ; if (do_i(i)) then
212 h_a = h(i,j,nz) + h_neglect
213 b_a = 1.0 / (h_a + d_a(i,nz)*mix(i,nz))
214 t_adj(i,j,nz) = b_a * (h_a*tv%T(i,j,nz) + mix(i,nz)*te_a(i,nz-1))
215 s_adj(i,j,nz) = b_a * (h_a*tv%S(i,j,nz) + mix(i,nz)*se_a(i,nz-1))
216 endif ; enddo
217 do k=nz-1,1,-1 ; do i=is,ie ; if (do_i(i)) then
218 t_adj(i,j,k) = te_a(i,k) + c_a(i,k+1)*t_adj(i,j,k+1)
219 s_adj(i,j,k) = se_a(i,k) + c_a(i,k+1)*s_adj(i,j,k+1)
220 endif ; enddo ; enddo
221
222 do i=is,ie ; if (do_i(i)) then
223 k = 1 ! A hook for debugging.
224 endif ; enddo
225
226 ! Do the final return pass on the columns where the penultimate pass was upward.
227 ! Also do a simple copy of T & S values on land points.
228 do i=is,ie
229 do_i(i) = ((g%mask2dT(i,j) > 0.0) .and. .not.last_down(i))
230 if (do_i(i)) then
231 h_b = h(i,j,1) + h_neglect
232 b_b = 1.0 / (h_b + d_b(i,2)*mix(i,2))
233 t_adj(i,j,1) = b_b * (h_b*tv%T(i,j,1) + mix(i,2)*te_b(i,2))
234 s_adj(i,j,1) = b_b * (h_b*tv%S(i,j,1) + mix(i,2)*se_b(i,2))
235 elseif (g%mask2dT(i,j) <= 0.0) then
236 t_adj(i,j,1) = tv%T(i,j,1) ; s_adj(i,j,1) = tv%S(i,j,1)
237 endif
238 enddo
239 do k=2,nz ; do i=is,ie
240 if (do_i(i)) then
241 t_adj(i,j,k) = te_b(i,k) + c_b(i,k)*t_adj(i,j,k-1)
242 s_adj(i,j,k) = se_b(i,k) + c_b(i,k)*s_adj(i,j,k-1)
243 elseif (g%mask2dT(i,j) <= 0.0) then
244 t_adj(i,j,k) = tv%T(i,j,k) ; s_adj(i,j,k) = tv%S(i,j,k)
245 endif
246 enddo ; enddo
247
248 do i=is,ie ; if (do_i(i)) then
249 k = 1 ! A hook for debugging.
250 endif ; enddo
251
252 enddo ! j-loop
253
254 k = 1 ! A hook for debugging.
255
256 ! The following set of expressions for the final values are derived from the partial
257 ! updates for the estimated temperatures and salinities around an interface, then directly
258 ! solving for the final temperatures and salinities. They are here for later reference
259 ! and to document an intermediate step in the stability calculation.
260 ! hp_a = (h_a + d_a(i,K-1)*mix(i,K-1))
261 ! hp_b = (h_b + d_b(i,K+1)*mix(i,K+1))
262 ! b2_c = 1.0 / (hp_a*hp_b + (hp_a + hp_b) * mix(i,K))
263 ! Th_a = h_a*tv%T(i,j,k-1) + mix(i,K-1)*Te_a(i,k-2)
264 ! Th_b = h_b*tv%T(i,j,k) + mix(i,K+1)*Te_b(i,k+1)
265 ! T_fin(i,k) = ( (hp_a + mix(i,K)) * Th_b + Th_a * mix(i,K) ) * b2_c
266 ! T_fin(i,k-1) = ( (hp_b + mix(i,K)) * Th_a + Th_b * mix(i,K) ) * b2_c
267 ! Sh_a = h_a*tv%S(i,j,k-1) + mix(i,K-1)*Se_a(i,k-2)
268 ! Sh_b = h_b*tv%S(i,j,k) + mix(i,K+1)*Se_b(i,k+1)
269 ! S_fin(i,k) = ( (hp_a + mix(i,K)) * Sh_b + Sh_a * mix(i,K) ) * b2_c
270 ! S_fin(i,k-1) = ( (hp_b + mix(i,K)) * Sh_a + Sh_b * mix(i,K) ) * b2_c
271
272end subroutine full_convection
273
274!> This function returns True if the profiles around the given interface will be
275!! statically unstable after mixing above below. The arguments are the ocean state
276!! above and below, including partial calculations from a tridiagonal solver.
277function is_unstable(dRho_dT, dRho_dS, h_a, h_b, mix_A, mix_B, T_a, T_b, S_a, S_b, &
278 Te_aa, Te_bb, Se_aa, Se_bb, d_A, d_B)
279 real, intent(in) :: drho_dt !< The derivative of in situ density with temperature [R C-1 ~> kg m-3 degC-1]
280 real, intent(in) :: drho_ds !< The derivative of in situ density with salinity [R S-1 ~> kg m-3 ppt-1]
281 real, intent(in) :: h_a !< The thickness of the layer above [H ~> m or kg m-2]
282 real, intent(in) :: h_b !< The thickness of the layer below [H ~> m or kg m-2]
283 real, intent(in) :: mix_a !< The time integrated mixing rate of the interface above [H ~> m or kg m-2]
284 real, intent(in) :: mix_b !< The time integrated mixing rate of the interface below [H ~> m or kg m-2]
285 real, intent(in) :: t_a !< The initial temperature of the layer above [C ~> degC]
286 real, intent(in) :: t_b !< The initial temperature of the layer below [C ~> degC]
287 real, intent(in) :: s_a !< The initial salinity of the layer below [S ~> ppt]
288 real, intent(in) :: s_b !< The initial salinity of the layer below [S ~> ppt]
289 real, intent(in) :: te_aa !< The estimated temperature two layers above rescaled by d_A [C ~> degC]
290 real, intent(in) :: te_bb !< The estimated temperature two layers below rescaled by d_B [C ~> degC]
291 real, intent(in) :: se_aa !< The estimated salinity two layers above rescaled by d_A [S ~> ppt]
292 real, intent(in) :: se_bb !< The estimated salinity two layers below rescaled by d_B [S ~> ppt]
293 real, intent(in) :: d_a !< The rescaling dependency across the interface above [nondim]
294 real, intent(in) :: d_b !< The rescaling dependency across the interface below [nondim]
295 logical :: is_unstable !< The return value, true if the profile is statically unstable
296 !! around the interface in question.
297
298 ! These expressions for the local stability are long, but they have been carefully
299 ! grouped for accuracy even when the mixing rates are huge or tiny, and common
300 ! positive definite factors that would appear in the final expression for the
301 ! locally referenced potential density difference across an interface have been omitted.
302 is_unstable = (drho_dt * ((h_a * h_b * (t_b - t_a) + &
303 mix_a*mix_b * (d_a*te_bb - d_b*te_aa)) + &
304 (h_a*mix_b * (te_bb - d_b*t_a) + &
305 h_b*mix_a * (d_a*t_b - te_aa)) ) + &
306 drho_ds * ((h_a * h_b * (s_b - s_a) + &
307 mix_a*mix_b * (d_a*se_bb - d_b*se_aa)) + &
308 (h_a*mix_b * (se_bb - d_b*s_a) + &
309 h_b*mix_a * (d_a*s_b - se_aa)) ) < 0.0)
310end function is_unstable
311
312!> Returns the partial derivatives of locally referenced potential density with
313!! temperature and salinity after the properties have been smoothed with a small
314!! constant diffusivity.
315subroutine smoothed_drdt_drds(h, dz, tv, Kddt, dR_dT, dR_dS, G, GV, US, j, p_surf, halo)
316 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
317 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
318 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
319 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
320 real, dimension(SZI_(G),SZK_(GV)), &
321 intent(in) :: dz !< Height change across layers [Z ~> m]
322 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
323 !! thermodynamic variables
324 real, intent(in) :: Kddt !< A diffusivity times a time increment [H Z ~> m2 or kg m-1].
325 real, dimension(SZI_(G),SZK_(GV)+1), &
326 intent(out) :: dR_dT !< Derivative of locally referenced
327 !! potential density with temperature [R C-1 ~> kg m-3 degC-1]
328 real, dimension(SZI_(G),SZK_(GV)+1), &
329 intent(out) :: dR_dS !< Derivative of locally referenced
330 !! potential density with salinity [R S-1 ~> kg m-3 ppt-1]
331 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
332 integer, intent(in) :: j !< The j-point to work on.
333 real, dimension(:,:), pointer :: p_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa].
334 integer, intent(in) :: halo !< Halo width over which to compute
335
336 ! Local variables
337 real :: mix(SZI_(G),SZK_(GV)+1) ! The diffusive mixing length (kappa*dt)/dz
338 ! between layers within in a timestep [H ~> m or kg m-2].
339 real :: b1(SZI_(G)) ! A tridiagonal solver variable [H-1 ~> m-1 or m2 kg-1]
340 real :: d1(SZI_(G)) ! A tridiagonal solver variable [nondim]
341 real :: c1(SZI_(G),SZK_(GV)) ! A tridiagonal solver variable [nondim]
342 real :: T_f(SZI_(G),SZK_(GV)) ! Filtered temperatures [C ~> degC]
343 real :: S_f(SZI_(G),SZK_(GV)) ! Filtered salinities [S ~> ppt]
344 real :: pres(SZI_(G)) ! Interface pressures [R L2 T-2 ~> Pa].
345 real :: T_EOS(SZI_(G)) ! Filtered and vertically averaged temperatures [C ~> degC]
346 real :: S_EOS(SZI_(G)) ! Filtered and vertically averaged salinities [S ~> ppt]
347 real :: kap_dt_x2 ! The product of 2*kappa*dt [H Z ~> m2 or kg m-1].
348 real :: dz_neglect, h0 ! A negligible vertical distances [Z ~> m]
349 real :: h_neglect ! A negligible thickness 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, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
353 integer :: i, k, is, ie, nz
354
355 is = g%isc-halo ; ie = g%iec+halo
356 nz = gv%ke
357
358 h_neglect = gv%H_subroundoff
359 dz_neglect = gv%dz_subroundoff
360 kap_dt_x2 = 2.0*kddt
361
362 if (kddt <= 0.0) then
363 do k=1,nz ; do i=is,ie
364 t_f(i,k) = tv%T(i,j,k) ; s_f(i,k) = tv%S(i,j,k)
365 enddo ; enddo
366 else
367 h0 = 1.0e-16*sqrt(gv%H_to_m*us%m_to_Z*kddt) + dz_neglect
368 do i=is,ie
369 mix(i,2) = kap_dt_x2 / ((dz(i,1)+dz(i,2)) + h0)
370
371 h_tr = h(i,j,1) + h_neglect
372 b1(i) = 1.0 / (h_tr + mix(i,2))
373 d1(i) = b1(i) * h(i,j,1)
374 t_f(i,1) = (b1(i)*h_tr)*tv%T(i,j,1)
375 s_f(i,1) = (b1(i)*h_tr)*tv%S(i,j,1)
376 enddo
377 do k=2,nz-1 ; do i=is,ie
378 mix(i,k+1) = kap_dt_x2 / ((dz(i,k)+dz(i,k+1)) + h0)
379
380 c1(i,k) = mix(i,k) * b1(i)
381 h_tr = h(i,j,k) + h_neglect
382 b1(i) = 1.0 / ((h_tr + d1(i)*mix(i,k)) + mix(i,k+1))
383 d1(i) = b1(i) * (h_tr + d1(i)*mix(i,k))
384 t_f(i,k) = b1(i) * (h_tr*tv%T(i,j,k) + mix(i,k)*t_f(i,k-1))
385 s_f(i,k) = b1(i) * (h_tr*tv%S(i,j,k) + mix(i,k)*s_f(i,k-1))
386 enddo ; enddo
387 do i=is,ie
388 c1(i,nz) = mix(i,nz) * b1(i)
389 h_tr = h(i,j,nz) + h_neglect
390 b1(i) = 1.0 / (h_tr + d1(i)*mix(i,nz))
391 t_f(i,nz) = b1(i) * (h_tr*tv%T(i,j,nz) + mix(i,nz)*t_f(i,nz-1))
392 s_f(i,nz) = b1(i) * (h_tr*tv%S(i,j,nz) + mix(i,nz)*s_f(i,nz-1))
393 enddo
394 do k=nz-1,1,-1 ; do i=is,ie
395 t_f(i,k) = t_f(i,k) + c1(i,k+1)*t_f(i,k+1)
396 s_f(i,k) = s_f(i,k) + c1(i,k+1)*s_f(i,k+1)
397 enddo ; enddo
398 endif
399
400 if (associated(p_surf)) then
401 do i=is,ie ; pres(i) = p_surf(i,j) ; enddo
402 else
403 do i=is,ie ; pres(i) = 0.0 ; enddo
404 endif
405 eosdom(:) = eos_domain(g%HI, halo)
406 call calculate_density_derivs(t_f(:,1), s_f(:,1), pres, dr_dt(:,1), dr_ds(:,1), tv%eqn_of_state, eosdom)
407 do i=is,ie ; pres(i) = pres(i) + h(i,j,1)*(gv%H_to_RZ*gv%g_Earth) ; enddo
408 do k=2,nz
409 do i=is,ie
410 t_eos(i) = 0.5*(t_f(i,k-1) + t_f(i,k))
411 s_eos(i) = 0.5*(s_f(i,k-1) + s_f(i,k))
412 enddo
413 call calculate_density_derivs(t_eos, s_eos, pres, dr_dt(:,k), dr_ds(:,k), tv%eqn_of_state, eosdom)
414 do i=is,ie ; pres(i) = pres(i) + h(i,j,k)*(gv%H_to_RZ*gv%g_Earth) ; enddo
415 enddo
416 call calculate_density_derivs(t_f(:,nz), s_f(:,nz), pres, dr_dt(:,nz+1), dr_ds(:,nz+1), &
417 tv%eqn_of_state, eosdom)
418
419end subroutine smoothed_drdt_drds
420
421end module mom_full_convection