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