MOM_isopycnal_slopes.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!> Calculations of isoneutral slopes and stratification.
7
8use mom_debugging, only : hchksum, uvchksum
9use mom_error_handler, only : mom_error, fatal
10use mom_grid, only : ocean_grid_type
14use mom_eos, only : calculate_density_derivs, calculate_density_second_derivs, eos_domain
15use mom_open_boundary, only : ocean_obc_type, obc_none
16use mom_open_boundary, only : obc_direction_e, obc_direction_w, obc_direction_n, obc_direction_s
17
18implicit none ; private
19
20#include <MOM_memory.h>
21
23
24! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
25! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
26! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
27! vary with the Boussinesq approximation, the Boussinesq variant is given first.
28
29contains
30
31!> Calculate isopycnal slopes, and optionally return other stratification dependent functions such as N^2
32!! and dz*S^2*g-prime used, or calculable from factors used, during the calculation.
33subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stanley, slope_x, slope_y, &
34 N2_u, N2_v, dzu, dzv, dzSxN, dzSyN, halo, OBC, OBC_N2, &
35 drdx_u, drdy_v, drdz_u, drdz_v)
36 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
37 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
38 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
39 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
40 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: e !< Interface heights [Z ~> m]
41 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
42 !! thermodynamic variables
43 real, intent(in) :: dt_kappa_smooth !< A smoothing vertical
44 !! diffusivity times a smoothing
45 !! timescale [H Z ~> m2 or kg m-1]
46 logical, intent(in) :: use_stanley !< turn on stanley param in slope
47 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: slope_x !< Isopycnal slope in i-dir [Z L-1 ~> nondim]
48 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(inout) :: slope_y !< Isopycnal slope in j-dir [Z L-1 ~> nondim]
49 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), &
50 optional, intent(inout) :: n2_u !< Brunt-Vaisala frequency squared at
51 !! interfaces between u-points [L2 Z-2 T-2 ~> s-2]
52 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), &
53 optional, intent(inout) :: n2_v !< Brunt-Vaisala frequency squared at
54 !! interfaces between v-points [L2 Z-2 T-2 ~> s-2]
55 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), &
56 optional, intent(inout) :: dzu !< Z-thickness at u-points [Z ~> m]
57 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), &
58 optional, intent(inout) :: dzv !< Z-thickness at v-points [Z ~> m]
59 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), &
60 optional, intent(inout) :: dzsxn !< Z-thickness times zonal slope contribution to
61 !! Eady growth rate at u-points. [Z T-1 ~> m s-1]
62 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), &
63 optional, intent(inout) :: dzsyn !< Z-thickness times meridional slope contrib. to
64 !! Eady growth rate at v-points. [Z T-1 ~> m s-1]
65 integer, optional, intent(in) :: halo !< Halo width over which to compute
66 type(ocean_obc_type), optional, pointer :: obc !< Open boundaries control structure.
67 logical, optional, intent(in) :: obc_n2 !< If present and true, use interior data
68 !! to calculate stratification at open boundary
69 !! condition faces.
70 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), &
71 optional, intent(inout) :: drdx_u !< Zonal density gradient at u
72 !! along surfaces of constant z
73 !! (not along isopycnals or
74 !! model interfaces) [R L-1 ~> kg m-4]
75 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), &
76 optional, intent(inout) :: drdy_v !< Meridional density gradient at v
77 !! along surfaces of constant z
78 !! (not along isopycnals or
79 !! model interfaces) [R L-1 ~> kg m-4]
80 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), &
81 optional, intent(inout) :: drdz_u !< Vertical density gradient
82 !! at u [R Z-1 ~> kg m-4]
83 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), &
84 optional, intent(inout) :: drdz_v !< Vertical density gradient
85 !! at v [R Z-1 ~> kg m-4]
86
87 ! Local variables
88 real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: &
89 t, & ! The temperature [C ~> degC], with the values in
90 ! in massless layers filled vertically by diffusion.
91 s ! The filled salinity [S ~> ppt], with the values in
92 ! in massless layers filled vertically by diffusion.
93 real, dimension(SZI_(G), SZJ_(G),SZK_(GV)+1) :: &
94 pres ! The pressure at an interface [R L2 T-2 ~> Pa].
95 real, dimension(SZI_(G)) :: scrap ! An array to pass to calculate_density_second_derivs() that is
96 ! set there but will be ignored, it is used simultaneously with four different
97 ! inconsistent units of [R S-1 C-1 ~> kg m-3 degC-1 ppt-1], [R S-2 ~> kg m-3 ppt-2],
98 ! [T2 S-1 L-2 ~> kg m-3 ppt-1 Pa-1] and [T2 C-1 L-2 ~> kg m-3 degC-1 Pa-1].
99 real, dimension(SZIB_(G)) :: &
100 drho_dt_u, & ! The derivative of density with temperature at u points [R C-1 ~> kg m-3 degC-1].
101 drho_ds_u ! The derivative of density with salinity at u points [R S-1 ~> kg m-3 ppt-1].
102 real, dimension(SZI_(G)) :: &
103 drho_dt_v, & ! The derivative of density with temperature at v points [R C-1 ~> kg m-3 degC-1].
104 drho_ds_v, & ! The derivative of density with salinity at v points [R S-1 ~> kg m-3 ppt-1].
105 drho_dt_dt_h, & ! The second derivative of density with temperature at h points [R C-2 ~> kg m-3 degC-2]
106 drho_dt_dt_hr ! The second derivative of density with temperature at h (+1) points [R C-2 ~> kg m-3 degC-2]
107 real, dimension(SZIB_(G)) :: &
108 t_u, & ! Temperature on the interface at the u-point [C ~> degC].
109 s_u, & ! Salinity on the interface at the u-point [S ~> ppt].
110 gxspv_u, & ! Gravitiational acceleration times the specific volume at an interface
111 ! at the u-points [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1]
112 pres_u ! Pressure on the interface at the u-point [R L2 T-2 ~> Pa].
113 real, dimension(SZI_(G)) :: &
114 t_v, & ! Temperature on the interface at the v-point [C ~> degC].
115 s_v, & ! Salinity on the interface at the v-point [S ~> ppt].
116 gxspv_v, & ! Gravitiational acceleration times the specific volume at an interface
117 ! at the v-points [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1]
118 pres_v ! Pressure on the interface at the v-point [R L2 T-2 ~> Pa].
119 real, dimension(SZI_(G)) :: &
120 t_h, & ! Temperature on the interface at the h-point [C ~> degC].
121 s_h, & ! Salinity on the interface at the h-point [S ~> ppt]
122 pres_h, & ! Pressure on the interface at the h-point [R L2 T-2 ~> Pa].
123 t_hr, & ! Temperature on the interface at the h (+1) point [C ~> degC].
124 s_hr, & ! Salinity on the interface at the h (+1) point [S ~> ppt]
125 pres_hr ! Pressure on the interface at the h (+1) point [R L2 T-2 ~> Pa].
126 real :: drdia, drdib ! Along layer zonal potential density gradients in the layers above (A)
127 ! and below (B) the interface times the grid spacing [R ~> kg m-3].
128 real :: drdja, drdjb ! Along layer meridional potential density gradients in the layers above (A)
129 ! and below (B) the interface times the grid spacing [R ~> kg m-3].
130 real :: drdkl, drdkr ! Vertical density differences across an interface [R ~> kg m-3].
131 real :: hg2a, hg2b ! Squares of geometric mean thicknesses [H2 ~> m2 or kg2 m-4].
132 real :: hg2l, hg2r ! Squares of geometric mean thicknesses [H2 ~> m2 or kg2 m-4].
133 real :: haa, hab, hal, har ! Arithmetic mean thicknesses [H ~> m or kg m-2].
134 real :: dzal, dzar ! Temporary thicknesses in eta units [Z ~> m].
135 real :: wta, wtb ! Unnormalized weights of the slopes above and below [H3 ~> m3 or kg3 m-6]
136 real :: wtl, wtr ! Unnormalized weights of the slopes to the left and right [H3 Z ~> m4 or kg3 m-5]
137 real :: drdx, drdy ! Zonal and meridional density gradients [R L-1 ~> kg m-4].
138 real :: drdz ! Vertical density gradient [R Z-1 ~> kg m-4].
139 real :: slope ! The slope of density surfaces, calculated in a way
140 ! that is always between -1 and 1. [Z L-1 ~> nondim]
141 real :: mag_grad2 ! The squared magnitude of the 3-d density gradient [R2 Z-2 ~> kg2 m-8].
142 real :: h_neglect ! A thickness that is so small it is usually lost
143 ! in roundoff and can be neglected [H ~> m or kg m-2].
144 real :: h_neglect2 ! h_neglect^2 [H2 ~> m2 or kg2 m-4].
145 real :: dz_neglect ! A change in interface heights that is so small it is usually lost
146 ! in roundoff and can be neglected [Z ~> m].
147 logical :: use_eos ! If true, density is calculated from T & S using an equation of state.
148 real :: g_rho0 ! The gravitational acceleration divided by density [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1]
149
150 logical :: present_n2_u, present_n2_v
151 logical :: present_drdx_u, present_drdy_v
152 logical :: local_open_u_bc, local_open_v_bc ! True if u- or v-face OBCs exist anywhere in the global domain.
153 logical :: obc_friendly ! If true, open boundary conditions are in use and only interior data should
154 ! be used to calculate N2 at OBC faces.
155 integer, dimension(2) :: eosdom_u ! The shifted I-computational domain to use for equation of
156 ! state calculations at u-points.
157 integer, dimension(2) :: eosdom_v ! The shifted i-computational domain to use for equation of
158 ! state calculations at v-points.
159 integer, dimension(2) :: eosdom_h1 ! The shifted i-computational domain to use for equation of
160 ! state calculations at h points with 1 extra halo point
161 integer :: is, ie, js, je, nz, isdb
162 integer :: i, j, k
163
164 if (present(halo)) then
165 is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
166 eosdom_h1(:) = eos_domain(g%HI, halo=halo+1)
167 else
168 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
169 eosdom_h1(:) = eos_domain(g%HI, halo=1)
170 endif
171 eosdom_u(1) = is-1 - (g%IsdB-1) ; eosdom_u(2) = ie - (g%IsdB-1)
172 eosdom_v(:) = eos_domain(g%HI, halo=halo)
173
174 nz = gv%ke ; isdb = g%IsdB
175
176
177 h_neglect = gv%H_subroundoff ; h_neglect2 = h_neglect**2
178 dz_neglect = gv%dZ_subroundoff
179
180 local_open_u_bc = .false.
181 local_open_v_bc = .false.
182 obc_friendly = .false.
183 if (present(obc)) then ; if (associated(obc)) then
184 local_open_u_bc = obc%open_u_BCs_exist_globally
185 local_open_v_bc = obc%open_v_BCs_exist_globally
186 if (present(obc_n2)) obc_friendly = obc_n2
187 endif ; endif
188
189 use_eos = associated(tv%eqn_of_state)
190
191 present_drdx_u = PRESENT(drdx_u)
192 present_drdy_v = PRESENT(drdy_v)
193 present_n2_u = PRESENT(n2_u)
194 present_n2_v = PRESENT(n2_v)
195 g_rho0 = gv%g_Earth / gv%Rho0
196 if (present_n2_u) then
197 do j=js,je ; do i=is-1,ie
198 n2_u(i,j,1) = 0.
199 n2_u(i,j,nz+1) = 0.
200 enddo ; enddo
201 endif
202 if (present_n2_v) then
203 do j=js-1,je ; do i=is,ie
204 n2_v(i,j,1) = 0.
205 n2_v(i,j,nz+1) = 0.
206 enddo ; enddo
207 endif
208 if (present(dzu)) then
209 do j=js,je ; do i=is-1,ie
210 dzu(i,j,1) = 0.
211 dzu(i,j,nz+1) = 0.
212 enddo ; enddo
213 endif
214 if (present(dzv)) then
215 do j=js-1,je ; do i=is,ie
216 dzv(i,j,1) = 0.
217 dzv(i,j,nz+1) = 0.
218 enddo ; enddo
219 endif
220 if (present(dzsxn)) then
221 do j=js,je ; do i=is-1,ie
222 dzsxn(i,j,1) = 0.
223 dzsxn(i,j,nz+1) = 0.
224 enddo ; enddo
225 endif
226 if (present(dzsyn)) then
227 do j=js-1,je ; do i=is,ie
228 dzsyn(i,j,1) = 0.
229 dzsyn(i,j,nz+1) = 0.
230 enddo ; enddo
231 endif
232 ! Set boundary values to zero, since they will be at places
233 ! where streamfunction would be zero.
234 if (present_drdx_u) then
235 do j=js,je ; do i=is-1,ie
236 drdx_u(i,j,1) = 0.
237 drdx_u(i,j,nz+1) = 0.
238 drdz_u(i,j,1) = 0.
239 drdz_u(i,j,nz+1) = 0.
240 enddo ; enddo
241 endif
242 if (present_drdy_v) then
243 do j=js-1,je ; do i=is,ie
244 drdy_v(i,j,1) = 0.
245 drdy_v(i,j,nz+1) = 0.
246 drdz_v(i,j,1) = 0.
247 drdz_v(i,j,nz+1) = 0.
248 enddo ; enddo
249 endif
250
251 if (use_eos) then
252 if (present(halo)) then
253 call vert_fill_ts(h, tv%T, tv%S, dt_kappa_smooth, t, s, g, gv, us, halo+1)
254 else
255 call vert_fill_ts(h, tv%T, tv%S, dt_kappa_smooth, t, s, g, gv, us, 1)
256 endif
257 endif
258
259 if ((use_eos .and. allocated(tv%SpV_avg) .and. (tv%valid_SpV_halo < 1)) .and. &
260 (present_n2_u .or. present(dzsxn) .or. present_n2_v .or. present(dzsyn))) then
261 if (tv%valid_SpV_halo < 0) then
262 call mom_error(fatal, "calc_isoneutral_slopes called in fully non-Boussinesq mode "//&
263 "with invalid values of SpV_avg.")
264 else
265 call mom_error(fatal, "calc_isoneutral_slopes called in fully non-Boussinesq mode "//&
266 "with insufficiently large SpV_avg halos of width 0 but 1 is needed.")
267 endif
268 endif
269
270 ! Find the maximum and minimum permitted streamfunction.
271 if (associated(tv%p_surf)) then
272 !$OMP parallel do default(shared)
273 do j=js-1,je+1 ; do i=is-1,ie+1
274 pres(i,j,1) = tv%p_surf(i,j)
275 enddo ; enddo
276 else
277 !$OMP parallel do default(shared)
278 do j=js-1,je+1 ; do i=is-1,ie+1
279 pres(i,j,1) = 0.0
280 enddo ; enddo
281 endif
282 !$OMP parallel do default(shared)
283 do j=js-1,je+1
284 do k=1,nz ; do i=is-1,ie+1
285 pres(i,j,k+1) = pres(i,j,k) + gv%g_Earth * gv%H_to_RZ * h(i,j,k)
286 enddo ; enddo
287 enddo
288
289 do i=is-1,ie
290 gxspv_u(i) = g_rho0 ! This will be changed if both use_EOS and allocated(tv%SpV_avg) are true
291 enddo
292 !$OMP parallel do default(none) shared(nz,is,ie,js,je,IsdB,use_EOS,G,GV,US,pres,T,S,tv,h,e, &
293 !$OMP h_neglect,dz_neglect,h_neglect2, &
294 !$OMP present_N2_u,G_Rho0,N2_u,slope_x,dzSxN,EOSdom_u,EOSdom_h1, &
295 !$OMP present_drdx_u,drdx_u,drdz_u, &
296 !$OMP local_open_u_BC,dzu,OBC,use_stanley,OBC_friendly) &
297 !$OMP private(drdiA,drdiB,drdkL,drdkR,pres_u,T_u,S_u, &
298 !$OMP drho_dT_u,drho_dS_u,hg2A,hg2B,hg2L,hg2R,haA, &
299 !$OMP drho_dT_dT_h,scrap,pres_h,T_h,S_h, &
300 !$OMP haB,haL,haR,dzaL,dzaR,wtA,wtB,wtL,wtR,drdz, &
301 !$OMP drdx,mag_grad2,slope) &
302 !$OMP firstprivate(GxSpV_u)
303 do j=js,je ; do k=nz,2,-1
304 if (.not.(use_eos)) then
305 drdia = 0.0 ; drdib = 0.0
306 drdkl = gv%Rlay(k)-gv%Rlay(k-1) ; drdkr = gv%Rlay(k)-gv%Rlay(k-1)
307 endif
308
309 ! Calculate the zonal isopycnal slope.
310 if (use_eos) then
311 do i=is-1,ie
312 pres_u(i) = 0.5*(pres(i,j,k) + pres(i+1,j,k))
313 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)))
314 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)))
315 enddo
316 if (obc_friendly) then
317 if (obc%u_E_OBCs_on_PE .and. (j>=obc%js_u_E_obc) .and. (j<=obc%je_u_E_obc)) then
318 do i = max(is-1, obc%Is_u_E_obc), min(ie, obc%Ie_u_E_obc)
319 if (obc%segnum_u(i,j) > 0) then ! OBC_DIRECTION_E
320 pres_u(i) = pres(i,j,k)
321 t_u(i) = 0.5*(t(i,j,k) + t(i,j,k-1))
322 s_u(i) = 0.5*(s(i,j,k) + s(i,j,k-1))
323 endif
324 enddo
325 endif
326 if (obc%u_W_OBCs_on_PE .and. (j>=obc%js_u_W_obc) .and. (j<=obc%je_u_W_obc)) then
327 do i = max(is-1, obc%Is_u_W_obc), min(ie, obc%Ie_u_W_obc)
328 if (obc%segnum_u(i,j) < 0) then ! OBC_DIRECTION_W
329 pres_u(i) = pres(i+1,j,k)
330 t_u(i) = 0.5*(t(i+1,j,k) + t(i+1,j,k-1))
331 s_u(i) = 0.5*(s(i+1,j,k) + s(i+1,j,k-1))
332 endif
333 enddo
334 endif
335 endif
336 call calculate_density_derivs(t_u, s_u, pres_u, drho_dt_u, drho_ds_u, &
337 tv%eqn_of_state, eosdom_u)
338 if (present_n2_u .or. (present(dzsxn))) then
339 if (allocated(tv%SpV_avg)) then
340 do i=is-1,ie
341 gxspv_u(i) = gv%g_Earth * 0.25* ((tv%SpV_avg(i,j,k) + tv%SpV_avg(i+1,j,k)) + &
342 (tv%SpV_avg(i,j,k-1) + tv%SpV_avg(i+1,j,k-1)))
343 enddo
344 endif
345 endif
346 endif
347
348 if (use_stanley) then
349 do i=is-1,ie+1
350 pres_h(i) = pres(i,j,k)
351 t_h(i) = 0.5*(t(i,j,k) + t(i,j,k-1))
352 s_h(i) = 0.5*(s(i,j,k) + s(i,j,k-1))
353 enddo
354 ! The second line below would correspond to arguments
355 ! drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, drho_dT_dP, &
356 call calculate_density_second_derivs(t_h, s_h, pres_h, &
357 scrap, scrap, drho_dt_dt_h, scrap, scrap, &
358 tv%eqn_of_state, dom=eosdom_h1)
359 endif
360
361 do i=is-1,ie
362 if (use_eos) then
363 ! Estimate the horizontal density gradients along layers.
364 drdia = drho_dt_u(i) * (t(i+1,j,k-1)-t(i,j,k-1)) + &
365 drho_ds_u(i) * (s(i+1,j,k-1)-s(i,j,k-1))
366 drdib = drho_dt_u(i) * (t(i+1,j,k)-t(i,j,k)) + &
367 drho_ds_u(i) * (s(i+1,j,k)-s(i,j,k))
368
369 ! Estimate the vertical density gradients times the grid spacing.
370 drdkl = (drho_dt_u(i) * (t(i,j,k)-t(i,j,k-1)) + &
371 drho_ds_u(i) * (s(i,j,k)-s(i,j,k-1)))
372 drdkr = (drho_dt_u(i) * (t(i+1,j,k)-t(i+1,j,k-1)) + &
373 drho_ds_u(i) * (s(i+1,j,k)-s(i+1,j,k-1)))
374 endif
375 if (use_stanley) then
376 ! Correction to the horizontal density gradient due to nonlinearity in
377 ! the EOS rectifying SGS temperature anomalies
378 drdia = drdia + 0.5 * ((drho_dt_dt_h(i+1) * tv%varT(i+1,j,k-1)) - &
379 (drho_dt_dt_h(i) * tv%varT(i,j,k-1)) )
380 drdib = drdib + 0.5 * ((drho_dt_dt_h(i+1) * tv%varT(i+1,j,k)) - &
381 (drho_dt_dt_h(i) * tv%varT(i,j,k)) )
382 endif
383
384 hg2a = h(i,j,k-1)*h(i+1,j,k-1) + h_neglect2
385 hg2b = h(i,j,k)*h(i+1,j,k) + h_neglect2
386 hg2l = h(i,j,k-1)*h(i,j,k) + h_neglect2
387 hg2r = h(i+1,j,k-1)*h(i+1,j,k) + h_neglect2
388 haa = 0.5*(h(i,j,k-1) + h(i+1,j,k-1)) + h_neglect
389 hab = 0.5*(h(i,j,k) + h(i+1,j,k)) + h_neglect
390 hal = 0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect
391 har = 0.5*(h(i+1,j,k-1) + h(i+1,j,k)) + h_neglect
392 if (gv%Boussinesq) then
393 dzal = hal * gv%H_to_Z ; dzar = har * gv%H_to_Z
394 else
395 dzal = 0.5*(e(i,j,k-1) - e(i,j,k+1)) + dz_neglect
396 dzar = 0.5*(e(i+1,j,k-1) - e(i+1,j,k+1)) + dz_neglect
397 endif
398 if (present(dzu)) dzu(i,j,k) = 0.5*( dzal + dzar )
399 ! Use the harmonic mean thicknesses to weight the horizontal gradients.
400 ! These unnormalized weights have been rearranged to minimize divisions.
401 wta = hg2a*hab ; wtb = hg2b*haa
402 wtl = hg2l*(har*dzar) ; wtr = hg2r*(hal*dzal)
403
404 drdz = ((wtl * drdkl) + (wtr * drdkr)) / ((dzal*wtl) + (dzar*wtr))
405 ! The expression for drdz above is mathematically equivalent to:
406 ! drdz = ((hg2L/haL) * drdkL/dzaL + (hg2R/haR) * drdkR/dzaR) / &
407 ! ((hg2L/haL) + (hg2R/haR))
408 ! which is an estimate of the gradient of density across geopotentials.
409 if (present_n2_u) then
410 if (obc_friendly) then ; if (obc%segnum_u(i,j) /= 0) then
411 if (obc%segnum_u(i,j) > 0) then ! OBC_DIRECTION_E
412 drdz = drdkl / dzal ! Note that drdz is not used for slopes at OBC faces.
413 if (use_eos .and. allocated(tv%SpV_avg)) &
414 gxspv_u(i) = gv%g_Earth * 0.5 * (tv%SpV_avg(i,j,k) + tv%SpV_avg(i,j,k-1))
415 elseif (obc%segnum_u(i,j) < 0) then ! OBC_DIRECTION_W
416 drdz = drdkr / dzar
417 if (use_eos .and. allocated(tv%SpV_avg)) &
418 gxspv_u(i) = gv%g_Earth * 0.5 * (tv%SpV_avg(i+1,j,k) + tv%SpV_avg(i+1,j,k-1))
419 endif
420 endif ; endif
421
422 n2_u(i,j,k) = gxspv_u(i) * drdz * g%mask2dCu(i,j) ! Square of buoyancy freq. [L2 Z-2 T-2 ~> s-2]
423 endif
424
425 if (use_eos) then
426 drdx = ((wta * drdia + wtb * drdib) / (wta + wtb) - &
427 drdz * (e(i,j,k)-e(i+1,j,k))) * g%IdxCu(i,j)
428
429 ! This estimate of slope is accurate for small slopes, but bounded
430 ! to be between -1 and 1.
431 mag_grad2 = (us%Z_to_L*drdx)**2 + drdz**2
432 if (mag_grad2 > 0.0) then
433 slope = drdx / sqrt(mag_grad2)
434 else ! Just in case mag_grad2 = 0 ever.
435 slope = 0.0
436 endif
437 else ! With .not.use_EOS, the layers are constant density.
438 slope = (e(i+1,j,k)-e(i,j,k)) * g%IdxCu(i,j)
439 endif
440
441 if (local_open_u_bc) then
442 if (obc%segnum_u(i,j) /= 0) then
443 if (obc%segment(abs(obc%segnum_u(i,j)))%open) then
444 slope = 0.
445 ! This and/or the masking code below is to make slopes match inside
446 ! land mask. Might not be necessary except for DEBUG output.
447! if (OBC%segnum_u(I,j) > 0) then ! OBC_DIRECTION_E
448! slope_x(I+1,j,K) = 0.
449! else
450! slope_x(I-1,j,K) = 0.
451! endif
452 endif
453 endif
454 slope = slope * max(g%mask2dT(i,j), g%mask2dT(i+1,j))
455 endif
456
457 slope_x(i,j,k) = slope
458 if (present_drdx_u) then
459 drdx_u(i,j,k) = drdx
460 drdz_u(i,j,k) = drdz
461 endif
462 if (present(dzsxn)) &
463 dzsxn(i,j,k) = sqrt( gxspv_u(i) * max(0., (wtl * ( dzal * drdkl )) &
464 + (wtr * ( dzar * drdkr ))) / (wtl + wtr) ) & ! dz * N
465 * abs(slope) * g%mask2dCu(i,j) ! x-direction contribution to S^2
466
467 enddo ! I
468 enddo ; enddo ! end of j-loop
469
470 do i=is,ie
471 gxspv_v(i) = g_rho0 !This will be changed if both use_EOS and allocated(tv%SpV_avg) are true
472 enddo
473 ! Calculate the meridional isopycnal slope.
474 !$OMP parallel do default(none) shared(nz,is,ie,js,je,IsdB,use_EOS,G,GV,US,pres,T,S,tv, &
475 !$OMP h,h_neglect,e,dz_neglect, &
476 !$OMP h_neglect2,present_N2_v,G_Rho0,N2_v,slope_y,dzSyN,EOSdom_v, &
477 !$OMP present_drdy_v,drdy_v,drdz_v, &
478 !$OMP dzv,local_open_v_BC,OBC,use_stanley,OBC_friendly) &
479 !$OMP private(drdjA,drdjB,drdkL,drdkR,pres_v,T_v,S_v, &
480 !$OMP drho_dT_v,drho_dS_v,hg2A,hg2B,hg2L,hg2R,haA, &
481 !$OMP drho_dT_dT_h,scrap,pres_h,T_h,S_h, &
482 !$OMP drho_dT_dT_hr,pres_hr,T_hr,S_hr, &
483 !$OMP haB,haL,haR,dzaL,dzaR,wtA,wtB,wtL,wtR,drdz, &
484 !$OMP drdy,mag_grad2,slope) &
485 !$OMP firstprivate(GxSpV_v)
486 do j=js-1,je ; do k=nz,2,-1
487 if (.not.(use_eos)) then
488 drdja = 0.0 ; drdjb = 0.0
489 drdkl = gv%Rlay(k)-gv%Rlay(k-1) ; drdkr = gv%Rlay(k)-gv%Rlay(k-1)
490 endif
491
492 if (use_eos) then
493 do i=is,ie
494 pres_v(i) = 0.5*(pres(i,j,k) + pres(i,j+1,k))
495 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)))
496 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)))
497 enddo
498 if (obc_friendly) then
499 if (obc%v_N_OBCs_on_PE .and. (j>=obc%Js_v_N_obc) .and. (j<=obc%Je_v_N_obc)) then
500 do i = max(is, obc%is_v_N_obc), min(ie, obc%ie_v_N_obc)
501 if (obc%segnum_v(i,j) > 0) then ! OBC_DIRECTION_N
502 pres_v(i) = pres(i,j,k)
503 t_v(i) = 0.5*(t(i,j,k) + t(i,j,k-1))
504 s_v(i) = 0.5*(s(i,j,k) + s(i,j,k-1))
505 endif
506 enddo
507 endif
508 if (obc%v_S_OBCs_on_PE .and. (j>=obc%Js_v_S_obc) .and. (j<=obc%Je_v_S_obc)) then
509 do i = max(is, obc%is_v_S_obc), min(ie, obc%ie_v_S_obc)
510 if (obc%segnum_v(i,j) < 0) then ! OBC_DIRECTION_S
511 pres_v(i) = pres(i,j+1,k)
512 t_v(i) = 0.5*(t(i,j+1,k) + t(i,j+1,k-1))
513 s_v(i) = 0.5*(s(i,j+1,k) + s(i,j+1,k-1))
514 endif
515 enddo
516 endif
517 endif
518 call calculate_density_derivs(t_v, s_v, pres_v, drho_dt_v, drho_ds_v, &
519 tv%eqn_of_state, eosdom_v)
520
521 if ((present_n2_v) .or. (present(dzsyn))) then
522 if (allocated(tv%SpV_avg)) then
523 do i=is,ie
524 gxspv_v(i) = gv%g_Earth * 0.25* ((tv%SpV_avg(i,j,k) + tv%SpV_avg(i,j+1,k)) + &
525 (tv%SpV_avg(i,j,k-1) + tv%SpV_avg(i,j+1,k-1)))
526 enddo
527 endif
528 endif
529 endif
530
531 if (use_stanley) then
532 do i=is,ie
533 pres_h(i) = pres(i,j,k)
534 t_h(i) = 0.5*(t(i,j,k) + t(i,j,k-1))
535 s_h(i) = 0.5*(s(i,j,k) + s(i,j,k-1))
536
537 pres_hr(i) = pres(i,j+1,k)
538 t_hr(i) = 0.5*(t(i,j+1,k) + t(i,j+1,k-1))
539 s_hr(i) = 0.5*(s(i,j+1,k) + s(i,j+1,k-1))
540 enddo
541 ! The second line below would correspond to arguments
542 ! drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, drho_dT_dP, &
543 call calculate_density_second_derivs(t_h, s_h, pres_h, &
544 scrap, scrap, drho_dt_dt_h, scrap, scrap, &
545 tv%eqn_of_state, dom=eosdom_v)
546 call calculate_density_second_derivs(t_hr, s_hr, pres_hr, &
547 scrap, scrap, drho_dt_dt_hr, scrap, scrap, &
548 tv%eqn_of_state, dom=eosdom_v)
549 endif
550 do i=is,ie
551 if (use_eos) then
552 ! Estimate the horizontal density gradients along layers.
553 drdja = drho_dt_v(i) * (t(i,j+1,k-1)-t(i,j,k-1)) + &
554 drho_ds_v(i) * (s(i,j+1,k-1)-s(i,j,k-1))
555 drdjb = drho_dt_v(i) * (t(i,j+1,k)-t(i,j,k)) + &
556 drho_ds_v(i) * (s(i,j+1,k)-s(i,j,k))
557
558 ! Estimate the vertical density gradients times the grid spacing.
559 drdkl = (drho_dt_v(i) * (t(i,j,k)-t(i,j,k-1)) + &
560 drho_ds_v(i) * (s(i,j,k)-s(i,j,k-1)))
561 drdkr = (drho_dt_v(i) * (t(i,j+1,k)-t(i,j+1,k-1)) + &
562 drho_ds_v(i) * (s(i,j+1,k)-s(i,j+1,k-1)))
563 endif
564 if (use_stanley) then
565 ! Correction to the horizontal density gradient due to nonlinearity in
566 ! the EOS rectifying SGS temperature anomalies
567 drdja = drdja + 0.5 * ((drho_dt_dt_hr(i) * tv%varT(i,j+1,k-1)) - &
568 (drho_dt_dt_h(i) * tv%varT(i,j,k-1)) )
569 drdjb = drdjb + 0.5 * ((drho_dt_dt_hr(i) * tv%varT(i,j+1,k)) - &
570 (drho_dt_dt_h(i) * tv%varT(i,j,k)) )
571 endif
572
573 hg2a = h(i,j,k-1)*h(i,j+1,k-1) + h_neglect2
574 hg2b = h(i,j,k)*h(i,j+1,k) + h_neglect2
575 hg2l = h(i,j,k-1)*h(i,j,k) + h_neglect2
576 hg2r = h(i,j+1,k-1)*h(i,j+1,k) + h_neglect2
577 haa = 0.5*(h(i,j,k-1) + h(i,j+1,k-1)) + h_neglect
578 hab = 0.5*(h(i,j,k) + h(i,j+1,k)) + h_neglect
579 hal = 0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect
580 har = 0.5*(h(i,j+1,k-1) + h(i,j+1,k)) + h_neglect
581 if (gv%Boussinesq) then
582 dzal = hal * gv%H_to_Z ; dzar = har * gv%H_to_Z
583 else
584 dzal = 0.5*(e(i,j,k-1) - e(i,j,k+1)) + dz_neglect
585 dzar = 0.5*(e(i,j+1,k-1) - e(i,j+1,k+1)) + dz_neglect
586 endif
587 if (present(dzv)) dzv(i,j,k) = 0.5*( dzal + dzar )
588 ! Use the harmonic mean thicknesses to weight the horizontal gradients.
589 ! These unnormalized weights have been rearranged to minimize divisions.
590 wta = hg2a*hab ; wtb = hg2b*haa
591 wtl = hg2l*(har*dzar) ; wtr = hg2r*(hal*dzal)
592
593 drdz = ((wtl * drdkl) + (wtr * drdkr)) / ((dzal*wtl) + (dzar*wtr))
594 ! The expression for drdz above is mathematically equivalent to:
595 ! drdz = ((hg2L/haL) * drdkL/dzaL + (hg2R/haR) * drdkR/dzaR) / &
596 ! ((hg2L/haL) + (hg2R/haR))
597 ! which is an estimate of the gradient of density across geopotentials.
598 if (present_n2_v) then
599 if (obc_friendly) then ; if (obc%segnum_v(i,j) /= 0) then
600 if (obc%segnum_v(i,j) > 0) then ! OBC_DIRECTION_N
601 drdz = drdkl / dzal ! Note that drdz is not used for slopes at OBC faces.
602 if (use_eos .and. allocated(tv%SpV_avg)) &
603 gxspv_v(i) = gv%g_Earth * 0.5 * (tv%SpV_avg(i,j,k) + tv%SpV_avg(i,j,k-1))
604 elseif (obc%segnum_v(i,j) < 0) then ! OBC_DIRECTION_S
605 drdz = drdkl / dzal
606 if (use_eos .and. allocated(tv%SpV_avg)) &
607 gxspv_v(i) = gv%g_Earth * 0.5 * (tv%SpV_avg(i,j+1,k) + tv%SpV_avg(i,j+1,k-1))
608 endif
609 endif ; endif
610
611 n2_v(i,j,k) = gxspv_v(i) * drdz * g%mask2dCv(i,j) ! Square of buoyancy freq. [L2 Z-2 T-2 ~> s-2]
612 endif
613
614 if (use_eos) then
615 drdy = ((wta * drdja + wtb * drdjb) / (wta + wtb) - &
616 drdz * (e(i,j,k)-e(i,j+1,k))) * g%IdyCv(i,j)
617
618 ! This estimate of slope is accurate for small slopes, but bounded
619 ! to be between -1 and 1.
620 mag_grad2 = (us%Z_to_L*drdy)**2 + drdz**2
621 if (mag_grad2 > 0.0) then
622 slope = drdy / sqrt(mag_grad2)
623 else ! Just in case mag_grad2 = 0 ever.
624 slope = 0.0
625 endif
626 else ! With .not.use_EOS, the layers are constant density.
627 slope = (e(i,j+1,k)-e(i,j,k)) * g%IdyCv(i,j)
628 endif
629
630 if (local_open_v_bc) then
631 if (obc%segnum_v(i,j) /= 0) then
632 if (obc%segment(abs(obc%segnum_v(i,j)))%open) then
633 slope = 0.
634 ! This and/or the masking code below is to make slopes match inside
635 ! land mask. Might not be necessary except for DEBUG output.
636! if (OBC%segnum_v(i,J)) > 0) then ! OBC_DIRECTION_N
637! slope_y(i,J+1,K) = 0.
638! else
639! slope_y(i,J-1,K) = 0.
640! endif
641 endif
642 endif
643 slope = slope * max(g%mask2dT(i,j), g%mask2dT(i,j+1))
644 endif
645 slope_y(i,j,k) = slope
646 if (present_drdy_v) then
647 drdy_v(i,j,k) = drdy
648 drdz_v(i,j,k) = drdz
649 endif
650 if (present(dzsyn)) &
651 dzsyn(i,j,k) = sqrt( gxspv_v(i) * max(0., (wtl * ( dzal * drdkl )) &
652 + (wtr * ( dzar * drdkr ))) / (wtl + wtr) ) & ! dz * N
653 * abs(slope) * g%mask2dCv(i,j) ! x-direction contribution to S^2
654
655 enddo ! i
656 enddo ; enddo ! end of j-loop
657
658end subroutine calc_isoneutral_slopes
659
660!> Returns tracer arrays (nominally T and S) with massless layers filled with
661!! sensible values, by diffusing vertically with a small but constant diffusivity.
662subroutine vert_fill_ts(h, T_in, S_in, kappa_dt, T_f, S_f, G, GV, US, halo_here, larger_h_denom)
663 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
664 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
665 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
666 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
667 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: t_in !< Input temperature [C ~> degC]
668 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: s_in !< Input salinity [S ~> ppt]
669 real, intent(in) :: kappa_dt !< A vertical diffusivity to use for smoothing
670 !! times a smoothing timescale [H Z ~> m2 or kg m-1]
671 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: t_f !< Filled temperature [C ~> degC]
672 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: s_f !< Filled salinity [S ~> ppt]
673 integer, optional, intent(in) :: halo_here !< Number of halo points to work on,
674 !! 0 by default
675 logical, optional, intent(in) :: larger_h_denom !< Present and true, add a large
676 !! enough minimal thickness in the denominator of
677 !! the flux calculations so that the fluxes are
678 !! never so large as eliminate the transmission
679 !! of information across groups of massless layers.
680 ! Local variables
681 real :: ent(szi_(g),szk_(gv)+1) ! The diffusive entrainment (kappa*dt)/dz
682 ! between layers in a timestep [H ~> m or kg m-2].
683 real :: b1(szi_(g)) ! A variable used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1]
684 real :: d1(szi_(g)) ! A variable used by the tridiagonal solver [nondim], d1 = 1 - c1.
685 real :: c1(szi_(g),szk_(gv)) ! A variable used by the tridiagonal solver [nondim].
686 real :: kap_dt_x2 ! The 2*kappa_dt converted to H units [H2 ~> m2 or kg2 m-4].
687 real :: h_neglect ! A negligible thickness [H ~> m or kg m-2], to allow for zero thicknesses.
688 real :: h0 ! A negligible thickness to allow for zero thickness layers without
689 ! completely decoupling groups of layers [H ~> m or kg m-2].
690 ! Often 0 < h_neglect << h0.
691 real :: h_tr ! h_tr is h at tracer points with a tiny thickness
692 ! added to ensure positive definiteness [H ~> m or kg m-2].
693 integer :: i, j, k, is, ie, js, je, nz, halo
694
695 halo=0 ; if (present(halo_here)) halo = max(halo_here,0)
696
697 is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo ; nz = gv%ke
698
699 h_neglect = gv%H_subroundoff
700 ! The use of the fixed rescaling factor in the next line avoids an extra call to thickness_to_dz()
701 ! and the use of an extra 3-d array of vertical distnaces across layers (dz). This would be more
702 ! physically consistent, but it would also be more expensive, and given that this routine applies
703 ! a small (but arbitrary) amount of mixing to clean up the properties of nearly massless layers,
704 ! the added expense is hard to justify.
705 kap_dt_x2 = (2.0*kappa_dt) * (us%Z_to_m*gv%m_to_H) ! Usually the latter term is GV%Z_to_H.
706 h0 = h_neglect
707 if (present(larger_h_denom)) then
708 if (larger_h_denom) h0 = 1.0e-16*sqrt(0.5*kap_dt_x2)
709 endif
710
711 if (kap_dt_x2 <= 0.0) then
712 !$OMP parallel do default(shared)
713 do k=1,nz ; do j=js,je ; do i=is,ie
714 t_f(i,j,k) = t_in(i,j,k) ; s_f(i,j,k) = s_in(i,j,k)
715 enddo ; enddo ; enddo
716 else
717 !$OMP parallel do default(shared) private(ent,b1,d1,c1,h_tr)
718 do j=js,je
719 do i=is,ie
720 ent(i,2) = kap_dt_x2 / ((h(i,j,1)+h(i,j,2)) + h0)
721 h_tr = h(i,j,1) + h_neglect
722 b1(i) = 1.0 / (h_tr + ent(i,2))
723 d1(i) = b1(i) * h_tr
724 t_f(i,j,1) = (b1(i)*h_tr)*t_in(i,j,1)
725 s_f(i,j,1) = (b1(i)*h_tr)*s_in(i,j,1)
726 enddo
727 do k=2,nz-1 ; do i=is,ie
728 ent(i,k+1) = kap_dt_x2 / ((h(i,j,k)+h(i,j,k+1)) + h0)
729 h_tr = h(i,j,k) + h_neglect
730 c1(i,k) = ent(i,k) * b1(i)
731 b1(i) = 1.0 / ((h_tr + d1(i)*ent(i,k)) + ent(i,k+1))
732 d1(i) = b1(i) * (h_tr + d1(i)*ent(i,k))
733 t_f(i,j,k) = b1(i) * (h_tr*t_in(i,j,k) + ent(i,k)*t_f(i,j,k-1))
734 s_f(i,j,k) = b1(i) * (h_tr*s_in(i,j,k) + ent(i,k)*s_f(i,j,k-1))
735 enddo ; enddo
736 do i=is,ie
737 c1(i,nz) = ent(i,nz) * b1(i)
738 h_tr = h(i,j,nz) + h_neglect
739 b1(i) = 1.0 / (h_tr + d1(i)*ent(i,nz))
740 t_f(i,j,nz) = b1(i) * (h_tr*t_in(i,j,nz) + ent(i,nz)*t_f(i,j,nz-1))
741 s_f(i,j,nz) = b1(i) * (h_tr*s_in(i,j,nz) + ent(i,nz)*s_f(i,j,nz-1))
742 enddo
743 do k=nz-1,1,-1 ; do i=is,ie
744 t_f(i,j,k) = t_f(i,j,k) + c1(i,k+1)*t_f(i,j,k+1)
745 s_f(i,j,k) = s_f(i,j,k) + c1(i,k+1)*s_f(i,j,k+1)
746 enddo ; enddo
747 enddo
748 endif
749
750end subroutine vert_fill_ts
751
752end module mom_isopycnal_slopes