portedportable, not yet portedexecuted, not portableexecutable, not hit by this run
| 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. | |
| 6 | module MOM_isopycnal_slopes | |
| 7 | ||
| 8 | use MOM_debugging, only : hchksum, uvchksum | |
| 9 | use MOM_error_handler, only : MOM_error, FATAL | |
| 10 | use MOM_grid, only : ocean_grid_type | |
| 11 | use MOM_unit_scaling, only : unit_scale_type | |
| 12 | use MOM_variables, only : thermo_var_ptrs | |
| 13 | use MOM_verticalGrid, only : verticalGrid_type | |
| 14 | use MOM_EOS, only : calculate_density_derivs, calculate_density_second_derivs, EOS_domain | |
| 15 | use MOM_open_boundary, only : ocean_OBC_type, OBC_NONE | |
| 16 | use MOM_open_boundary, only : OBC_DIRECTION_E, OBC_DIRECTION_W, OBC_DIRECTION_N, OBC_DIRECTION_S | |
| 17 | ||
| 18 | implicit none ; private | |
| 19 | ||
| 20 | #include <MOM_memory.h> | |
| 21 | #include "do_concurrent_compat.h" | |
| 22 | ||
| 23 | public calc_isoneutral_slopes, vert_fill_TS | |
| 24 | ||
| 25 | ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional | |
| 26 | ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with | |
| 27 | ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units | |
| 28 | ! vary with the Boussinesq approximation, the Boussinesq variant is given first. | |
| 29 | ||
| 30 | contains | |
| 31 | ||
| 32 | !> Calculate isopycnal slopes, and optionally return other stratification dependent functions such as N^2 | |
| 33 | !! and dz*S^2*g-prime used, or calculable from factors used, during the calculation. | |
| 34 | 24 | subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stanley, slope_x, slope_y, & |
| 35 | 24 | niblock, njblock, nkblock, N2_u, N2_v, dzu, dzv, dzSxN, dzSyN, halo, & |
| 36 | OBC, OBC_N2) | |
| 37 | integer, intent(in) :: niblock !< Blocksize in i direction | |
| 38 | integer, intent(in) :: njblock !< Blocksize in j direction | |
| 39 | integer, intent(in) :: nkblock !< Blocksize in k direction | |
| 40 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure | |
| 41 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure | |
| 42 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 43 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] | |
| 44 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: e !< Interface heights [Z ~> m] | |
| 45 | type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various | |
| 46 | !! thermodynamic variables | |
| 47 | real, intent(in) :: dt_kappa_smooth !< A smoothing vertical | |
| 48 | !! diffusivity times a smoothing | |
| 49 | !! timescale [H Z ~> m2 or kg m-1] | |
| 50 | logical, intent(in) :: use_stanley !< turn on stanley param in slope | |
| 51 | real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: slope_x !< Isopycnal slope in i-dir [Z L-1 ~> nondim] | |
| 52 | real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(inout) :: slope_y !< Isopycnal slope in j-dir [Z L-1 ~> nondim] | |
| 53 | real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), & | |
| 54 | optional, intent(inout) :: N2_u !< Brunt-Vaisala frequency squared at | |
| 55 | !! interfaces between u-points [L2 Z-2 T-2 ~> s-2] | |
| 56 | real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), & | |
| 57 | optional, intent(inout) :: N2_v !< Brunt-Vaisala frequency squared at | |
| 58 | !! interfaces between v-points [L2 Z-2 T-2 ~> s-2] | |
| 59 | real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), & | |
| 60 | optional, intent(inout) :: dzu !< Z-thickness at u-points [Z ~> m] | |
| 61 | real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), & | |
| 62 | optional, intent(inout) :: dzv !< Z-thickness at v-points [Z ~> m] | |
| 63 | real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), & | |
| 64 | optional, intent(inout) :: dzSxN !< Z-thickness times zonal slope contribution to | |
| 65 | !! Eady growth rate at u-points. [Z T-1 ~> m s-1] | |
| 66 | real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), & | |
| 67 | optional, intent(inout) :: dzSyN !< Z-thickness times meridional slope contrib. to | |
| 68 | !! Eady growth rate at v-points. [Z T-1 ~> m s-1] | |
| 69 | integer, optional, intent(in) :: halo !< Halo width over which to compute | |
| 70 | type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure. | |
| 71 | logical, optional, intent(in) :: OBC_N2 !< If present and true, use interior data | |
| 72 | !! to calculate stratification at open boundary | |
| 73 | !! condition faces. | |
| 74 | ||
| 75 | ! Local variables | |
| 76 | real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: & | |
| 77 | 24 | T, & ! The temperature [C ~> degC], with the values in |
| 78 | ! in massless layers filled vertically by diffusion. | |
| 79 | 24 | S ! The filled salinity [S ~> ppt], with the values in |
| 80 | ! in massless layers filled vertically by diffusion. | |
| 81 | real, dimension(SZI_(G), SZJ_(G),SZK_(GV)+1) :: & | |
| 82 | 24 | pres ! The pressure at an interface [R L2 T-2 ~> Pa]. |
| 83 | ||
| 84 | ! Blocked work arrays. | |
| 85 | real, dimension(niblock, njblock, nkblock) :: & | |
| 86 | 24 | drho_dT, & ! The derivative of density with temperature [R C-1 ~> kg m-3 degC-1]. |
| 87 | 24 | drho_dS, & ! The derivative of density with salinity [R S-1 ~> kg m-3 ppt-1]. |
| 88 | 24 | drho_dT_dT_h, & ! The second derivative of density with temperature at h points [R C-2 ~> kg m-3 degC-2] |
| 89 | 24 | T_uvh, & ! Temperature at the u, v or h-points for derivative calculations [C ~> degC]. |
| 90 | 24 | S_uvh, & ! Salinity at the u, v or h-points for derivative calculations [S ~> ppt]. |
| 91 | 24 | GxSpV_uvh, & ! g * specific volume at u/v-point interfaces [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1] |
| 92 | 24 | pres_uvh ! Pressure at the u, v or h-points [R L2 T-2 ~> Pa]. |
| 93 | 24 | real, dimension(niblock) :: scrap ! Ignored output from calculate_density_second_derivs() |
| 94 | ||
| 95 | real :: drdiA, drdiB ! Along layer zonal potential density gradients in the layers above (A) | |
| 96 | ! and below (B) the interface times the grid spacing [R ~> kg m-3]. | |
| 97 | real :: drdjA, drdjB ! Along layer meridional potential density gradients in the layers above (A) | |
| 98 | ! and below (B) the interface times the grid spacing [R ~> kg m-3]. | |
| 99 | real :: drdkL, drdkR ! Vertical density differences across an interface [R ~> kg m-3]. | |
| 100 | real :: hg2A, hg2B ! Squares of geometric mean thicknesses [H2 ~> m2 or kg2 m-4]. | |
| 101 | real :: hg2L, hg2R ! Squares of geometric mean thicknesses [H2 ~> m2 or kg2 m-4]. | |
| 102 | real :: haA, haB, haL, haR ! Arithmetic mean thicknesses [H ~> m or kg m-2]. | |
| 103 | real :: dzaL, dzaR ! Temporary thicknesses in eta units [Z ~> m]. | |
| 104 | real :: wtA, wtB ! Unnormalized weights of the slopes above and below [H3 ~> m3 or kg3 m-6] | |
| 105 | real :: wtL, wtR ! Unnormalized weights of the slopes to the left and right [H3 Z ~> m4 or kg3 m-5] | |
| 106 | real :: drdx, drdy ! Zonal and meridional density gradients [R L-1 ~> kg m-4]. | |
| 107 | real :: drdz ! Vertical density gradient [R Z-1 ~> kg m-4]. | |
| 108 | real :: slope ! The slope of density surfaces, calculated in a way | |
| 109 | ! that is always between -1 and 1. [Z L-1 ~> nondim] | |
| 110 | real :: mag_grad2 ! The squared magnitude of the 3-d density gradient [R2 Z-2 ~> kg2 m-8]. | |
| 111 | real :: h_neglect ! A thickness that is so small it is usually lost | |
| 112 | ! in roundoff and can be neglected [H ~> m or kg m-2]. | |
| 113 | real :: h_neglect2 ! h_neglect^2 [H2 ~> m2 or kg2 m-4]. | |
| 114 | real :: dz_neglect ! A change in interface heights that is so small it is usually lost | |
| 115 | ! in roundoff and can be neglected [Z ~> m]. | |
| 116 | logical :: use_EOS ! If true, density is calculated from T & S using an equation of state. | |
| 117 | real :: G_Rho0 ! The gravitational acceleration divided by density [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1] | |
| 118 | ||
| 119 | logical :: present_N2_u, present_N2_v | |
| 120 | logical :: local_open_u_BC, local_open_v_BC ! True if u- or v-face OBCs exist anywhere in the global domain. | |
| 121 | logical :: OBC_friendly ! If true, open boundary conditions are in use and only interior data should | |
| 122 | ! be used to calculate N2 at OBC faces. | |
| 123 | integer, dimension(3,2) :: EOSdom_block !< 1-based EOS domain for the current block [nondim]. | |
| 124 | integer, dimension(2) :: EOSdom_block_h1 !< 1-based EOS domain for h-point fills (one extra column) [nondim]. | |
| 125 | integer :: is, ie, js, je, nz, IsdB | |
| 126 | integer :: i, j, k | |
| 127 | integer :: istart, iend !< First and last global i (or I) indices of the current block. | |
| 128 | integer :: jstart, jend !< First and last global j (or J) indices of the current block. | |
| 129 | integer :: kstart, kend !< First and last global K (or K) indices of the current block. | |
| 130 | integer :: ii, jj, kk !< Block-local 1-based i, j, and k indices. | |
| 131 | integer :: delta_i, delta_j | |
| 132 | integer :: ie_read, je_read !< Read-only extent of the h-point block used to supply the | |
| 133 | !! ii+1 (or jj+1) access needed by the Stanley stencil; equal | |
| 134 | !! to ie (or je) plus one column/row when use_stanley is true, | |
| 135 | !! and otherwise equal to ie (or je). | |
| 136 | integer :: iend_stanley, jend_stanley !< Last i (or j) index filled into the Stanley h-point block, | |
| 137 | !! as opposed to iend/jend, which bound what is actually written | |
| 138 | !! to the output arrays this block. | |
| 139 | ||
| 140 | ! Allocate locals on device | |
| 141 | !$omp target enter data map(alloc: T, S, pres, T_uvh, S_uvh, pres_uvh, GxSpV_uvh, & | |
| 142 | !$omp & scrap, drho_dT, drho_dS, drho_dT_dT_h) | |
| 143 | ||
| 144 | 12 | if (present(halo)) then |
| 145 | 12 | is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo |
| 146 | else | |
| 147 | 0 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec |
| 148 | endif | |
| 149 | ||
| 150 | 12 | nz = GV%ke ; IsdB = G%IsdB |
| 151 | ||
| 152 | 12 | h_neglect = GV%H_subroundoff ; h_neglect2 = h_neglect**2 |
| 153 | 12 | dz_neglect = GV%dZ_subroundoff |
| 154 | ||
| 155 | 12 | local_open_u_BC = .false. |
| 156 | 12 | local_open_v_BC = .false. |
| 157 | 12 | OBC_friendly = .false. |
| 158 | 12 | if (present(OBC)) then ; if (associated(OBC)) then |
| 159 | 0 | local_open_u_BC = OBC%open_u_BCs_exist_globally |
| 160 | 0 | local_open_v_BC = OBC%open_v_BCs_exist_globally |
| 161 | 0 | if (present(OBC_N2)) OBC_friendly = OBC_N2 |
| 162 | endif ; endif | |
| 163 | ||
| 164 | 12 | use_EOS = associated(tv%eqn_of_state) |
| 165 | ||
| 166 | 12 | present_N2_u = PRESENT(N2_u) |
| 167 | 12 | present_N2_v = PRESENT(N2_v) |
| 168 | 12 | G_Rho0 = GV%g_Earth / GV%Rho0 |
| 169 | 12 | if (present_N2_u) then |
| 170 | 12 | do concurrent (j=js:je, I=is-1:ie) |
| 171 | 91512 | N2_u(I,j,1) = 0. |
| 172 | 93000 | N2_u(I,j,nz+1) = 0. |
| 173 | enddo | |
| 174 | endif | |
| 175 | 12 | if (present_N2_v) then |
| 176 | 12 | do concurrent (J=js-1:je, i=is:ie) |
| 177 | 92232 | N2_v(i,J,1) = 0. |
| 178 | 93708 | N2_v(i,J,nz+1) = 0. |
| 179 | enddo | |
| 180 | endif | |
| 181 | 12 | if (present(dzu)) then |
| 182 | 0 | do concurrent (j=js:je, I=is-1:ie) |
| 183 | 0 | dzu(I,j,1) = 0. |
| 184 | 0 | dzu(I,j,nz+1) = 0. |
| 185 | enddo | |
| 186 | endif | |
| 187 | 12 | if (present(dzv)) then |
| 188 | 0 | do concurrent (J=js-1:je, i=is:ie) |
| 189 | 0 | dzv(i,J,1) = 0. |
| 190 | 0 | dzv(i,J,nz+1) = 0. |
| 191 | enddo | |
| 192 | endif | |
| 193 | 12 | if (present(dzSxN)) then |
| 194 | 0 | do concurrent (j=js:je, I=is-1:ie) |
| 195 | 0 | dzSxN(I,j,1) = 0. |
| 196 | 0 | dzSxN(I,j,nz+1) = 0. |
| 197 | enddo | |
| 198 | endif | |
| 199 | 12 | if (present(dzSyN)) then |
| 200 | 0 | do concurrent (J=js-1:je, i=is:ie) |
| 201 | 0 | dzSyN(i,J,1) = 0. |
| 202 | 0 | dzSyN(i,J,nz+1) = 0. |
| 203 | enddo | |
| 204 | endif | |
| 205 | ||
| 206 | 12 | if (use_EOS) then |
| 207 | 12 | if (present(halo)) then |
| 208 | 12 | call vert_fill_TS(h, tv%T, tv%S, dt_kappa_smooth, T, S, G, GV, US, halo+1) |
| 209 | else | |
| 210 | 0 | call vert_fill_TS(h, tv%T, tv%S, dt_kappa_smooth, T, S, G, GV, US, 1) |
| 211 | endif | |
| 212 | endif | |
| 213 | ||
| 214 | 12 | if ((use_EOS .and. allocated(tv%SpV_avg) .and. (tv%valid_SpV_halo < 1)) .and. & |
| 215 | (present_N2_u .or. present(dzSxN) .or. present_N2_v .or. present(dzSyN))) then | |
| 216 | 0 | if (tv%valid_SpV_halo < 0) then |
| 217 | call MOM_error(FATAL, "calc_isoneutral_slopes called in fully non-Boussinesq mode "//& | |
| 218 | 0 | "with invalid values of SpV_avg.") |
| 219 | else | |
| 220 | call MOM_error(FATAL, "calc_isoneutral_slopes called in fully non-Boussinesq mode "//& | |
| 221 | 0 | "with insufficiently large SpV_avg halos of width 0 but 1 is needed.") |
| 222 | endif | |
| 223 | endif | |
| 224 | ||
| 225 | ! Find the maximum and minimum permitted streamfunction. | |
| 226 | 12 | if (associated(tv%p_surf)) then |
| 227 | 1500 | do concurrent (j=js-1:je+1, i=is-1:ie+1) |
| 228 | 96732 | pres(i,j,1) = tv%p_surf(i,j) |
| 229 | enddo | |
| 230 | else | |
| 231 | 0 | do concurrent (j=js-1:je+1, i=is-1:ie+1) |
| 232 | 0 | pres(i,j,1) = 0.0 |
| 233 | enddo | |
| 234 | endif | |
| 235 | ||
| 236 | 12 | do concurrent (j=js-1:je+1) |
| 237 | 58380 | do k=1,nz |
| 238 | 768 | do concurrent (i=is-1:ie+1) |
| 239 | 7200000 | pres(i,j,K+1) = pres(i,j,K) + GV%g_Earth * GV%H_to_RZ * h(i,j,k) |
| 240 | enddo | |
| 241 | enddo | |
| 242 | enddo | |
| 243 | ||
| 244 | 12 | do concurrent(kk=1:nkblock, jj=1:njblock, ii=1:niblock) |
| 245 | 4332 | GxSpV_uvh(ii,jj,kk) = G_Rho0 ! This will be changed if both use_EOS and allocated(tv%SpV_avg) are true |
| 246 | enddo | |
| 247 | ||
| 248 | ! Stanley param needs access to h-point element ii+1, so when using stanley param, | |
| 249 | ! iterate in chunks of size block - 1 but read chunks of size block to ensure blocks | |
| 250 | ! always included needed elements. Allow access to element ie+1 so element ie is filled | |
| 251 | ! Block sizes are validated to be > 0 and not equal to 1 if use_stanley at initialization | |
| 252 | 12 | if (use_stanley) then |
| 253 | 0 | delta_i = niblock - 1 |
| 254 | 0 | ie_read = ie + 1 |
| 255 | else | |
| 256 | 12 | delta_i = niblock |
| 257 | 12 | ie_read = ie |
| 258 | endif | |
| 259 | ||
| 260 | ! Calculate the zonal isopycnal slope. | |
| 261 | 166956 | do kstart=2,nz,nkblock ; do jstart=js,je,njblock ; do istart=is-1,ie,delta_i |
| 262 | 110112 | iend = min(istart + delta_i - 1, ie) |
| 263 | 110112 | iend_stanley = min(istart + niblock - 1, ie_read) |
| 264 | 110112 | jend = min(jstart + njblock - 1, je) |
| 265 | 110112 | kend = min(kstart + nkblock - 1, nz) |
| 266 | 110112 | EOSdom_block(1,1) = 1 ; EOSdom_block(1,2) = iend - istart + 1 |
| 267 | 110112 | EOSdom_block(2,1) = 1 ; EOSdom_block(2,2) = jend - jstart + 1 |
| 268 | 110112 | EOSdom_block(3,1) = 1 ; EOSdom_block(3,2) = kend - kstart + 1 |
| 269 | ||
| 270 | 110112 | if (use_EOS) then |
| 271 | ! Fill block T_uvh/S_uvh/pres_uvh at u-points | |
| 272 | do concurrent(kk=1:kend-kstart+1, jj=1:jend-jstart+1, ii=1:iend-istart+1) & | |
| 273 | 13653888 | DO_LOCALITY(local(i,j,k)) |
| 274 | 6771888 | i = istart + ii - 1 |
| 275 | 6771888 | j = jstart + jj - 1 |
| 276 | 6771888 | k = kstart + kk - 1 |
| 277 | 6771888 | pres_uvh(ii,jj,kk) = 0.5*(pres(i,j,K) + pres(i+1,j,K)) |
| 278 | 6771888 | T_uvh(ii,jj,kk) = 0.25*((T(i,j,k) + T(i+1,j,k)) + (T(i,j,k-1) + T(i+1,j,k-1))) |
| 279 | 20425776 | S_uvh(ii,jj,kk) = 0.25*((S(i,j,k) + S(i+1,j,k)) + (S(i,j,k-1) + S(i+1,j,k-1))) |
| 280 | enddo | |
| 281 | ||
| 282 | 110112 | if (OBC_friendly) then |
| 283 | !$omp target update from(T_uvh, S_uvh, pres_uvh) | |
| 284 | ||
| 285 | 0 | if (OBC%u_E_OBCs_on_PE) then |
| 286 | 0 | do k = kstart,kend |
| 287 | 0 | kk = k - kstart + 1 |
| 288 | 0 | do j = max(jstart, OBC%js_u_E_obc), min(jend, OBC%je_u_E_obc) |
| 289 | 0 | jj = j - jstart + 1 |
| 290 | 0 | do i = max(istart, OBC%Is_u_E_obc), min(iend, OBC%Ie_u_E_obc) |
| 291 | 0 | ii = i - istart + 1 |
| 292 | 0 | if (OBC%segnum_u(i,j) > 0) then ! OBC_DIRECTION_E |
| 293 | 0 | pres_uvh(ii,jj,kk) = pres(i,j,k) |
| 294 | 0 | T_uvh(ii,jj,kk) = 0.5*(T(i,j,k) + T(i,j,k-1)) |
| 295 | 0 | S_uvh(ii,jj,kk) = 0.5*(S(i,j,k) + S(i,j,k-1)) |
| 296 | endif | |
| 297 | enddo | |
| 298 | enddo | |
| 299 | enddo | |
| 300 | endif | |
| 301 | 0 | if (OBC%u_W_OBCs_on_PE) then |
| 302 | 0 | do k = kstart, kend |
| 303 | 0 | kk = k - kstart + 1 |
| 304 | 0 | do j = max(jstart, OBC%js_u_W_obc), min(jend, OBC%je_u_W_obc) |
| 305 | 0 | jj = j - jstart + 1 |
| 306 | 0 | do i = max(istart, OBC%Is_u_W_obc), min(iend, OBC%Ie_u_W_obc) |
| 307 | 0 | ii = i - istart + 1 |
| 308 | 0 | if (OBC%segnum_u(i,j) < 0) then ! OBC_DIRECTION_W |
| 309 | 0 | pres_uvh(ii,jj,kk) = pres(i+1,j,k) |
| 310 | 0 | T_uvh(ii,jj,kk) = 0.5*(T(i+1,j,k) + T(i+1,j,k-1)) |
| 311 | 0 | S_uvh(ii,jj,kk) = 0.5*(S(i+1,j,k) + S(i+1,j,k-1)) |
| 312 | endif | |
| 313 | enddo | |
| 314 | enddo | |
| 315 | enddo | |
| 316 | endif | |
| 317 | ||
| 318 | !$omp target update to(T_uvh, S_uvh, pres_uvh) | |
| 319 | endif | |
| 320 | ||
| 321 | ! Pre-fill GxSpV at u-points: four-cell SpV average over the u-face and the two | |
| 322 | ! adjacent layers. Individual OBC faces may override this inside the slope loop. | |
| 323 | 110112 | if (present_N2_u .or. present(dzSxN)) then |
| 324 | 110112 | if (allocated(tv%SpV_avg)) then |
| 325 | do concurrent(kk=1:kend-kstart+1, jj=1:jend-jstart+1, II=1:iend-istart+1) & | |
| 326 | 0 | DO_LOCALITY(local(I,j,k)) |
| 327 | 0 | I = istart + II - 1 |
| 328 | 0 | j = jstart + jj - 1 |
| 329 | 0 | k = kstart + kk - 1 |
| 330 | GxSpV_uvh(II,jj,kk) = GV%g_Earth * 0.25 * ((tv%SpV_avg(I,j,k) + tv%SpV_avg(I+1,j,k)) + & | |
| 331 | 0 | (tv%SpV_avg(I,j,k-1) + tv%SpV_avg(I+1,j,k-1))) |
| 332 | enddo | |
| 333 | endif | |
| 334 | endif | |
| 335 | ||
| 336 | call calculate_density_derivs(T_uvh, S_uvh, pres_uvh, & | |
| 337 | drho_dT, drho_dS, & | |
| 338 | 110112 | tv%eqn_of_state, EOSdom_block) |
| 339 | ||
| 340 | 110112 | if (use_stanley) then |
| 341 | ! Reset T_uvh/S_uvh/pres_uvh at h-points for the Stanley second-derivative calculation. | |
| 342 | ! This loop fills all niblock elements of the _uvh arrays and can access index ie+1 of T | |
| 343 | ! and S to ensure there always exists an ii+1 element for the compute loop below. | |
| 344 | 0 | EOSdom_block_h1(1) = 1 ; EOSdom_block_h1(2) = iend_stanley - istart + 1 |
| 345 | do concurrent(kk=1:kend-kstart+1, jj=1:jend-jstart+1, ii=1:iend_stanley-istart+1) & | |
| 346 | 0 | DO_LOCALITY(local(i,j,k)) |
| 347 | 0 | i = istart + ii - 1 |
| 348 | 0 | j = jstart + jj - 1 |
| 349 | 0 | k = kstart + kk - 1 |
| 350 | 0 | pres_uvh(ii,jj,KK) = pres(i,j,K) |
| 351 | 0 | T_uvh(ii,jj,kk) = 0.5*(T(i,j,k) + T(i,j,k-1)) |
| 352 | 0 | S_uvh(ii,jj,kk) = 0.5*(S(i,j,k) + S(i,j,k-1)) |
| 353 | enddo | |
| 354 | ||
| 355 | ! TODO: Make 3D version of this | |
| 356 | 0 | do kk=1,kend-kstart+1 ; do jj=1,jend-jstart+1 |
| 357 | ! The second line below would correspond to arguments | |
| 358 | ! drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, drho_dT_dP, & | |
| 359 | call calculate_density_second_derivs(T_uvh(:,jj,kk), S_uvh(:,jj,kk), pres_uvh(:,jj,kk), & | |
| 360 | scrap, scrap, drho_dT_dT_h(:,jj,kk), scrap, scrap, & | |
| 361 | 0 | tv%eqn_of_state, dom=EOSdom_block_h1) |
| 362 | enddo ; enddo | |
| 363 | endif ! end use_stanley | |
| 364 | ||
| 365 | endif ! end use_EOS for zonal block | |
| 366 | ||
| 367 | ! Zonal slope compute over the block | |
| 368 | do concurrent(kk=1:kend-kstart+1, jj=1:jend-jstart+1, II=1:iend-istart+1) & | |
| 369 | 55056 | DO_LOCALITY(local(drdkL, drdkR, drdiA, drdiB, I, j)) |
| 370 | 6771888 | I = istart + II - 1 |
| 371 | 6771888 | j = jstart + jj - 1 |
| 372 | 6771888 | k = kstart + kk - 1 |
| 373 | ||
| 374 | 6771888 | if (use_EOS) then |
| 375 | ! Estimate the horizontal density gradients along layers. | |
| 376 | drdiA = drho_dT(ii,jj,kk) * (T(i+1,j,k-1)-T(i,j,k-1)) + & | |
| 377 | 6771888 | drho_dS(ii,jj,kk) * (S(i+1,j,k-1)-S(i,j,k-1)) |
| 378 | drdiB = drho_dT(ii,jj,kk) * (T(i+1,j,k)-T(i,j,k)) + & | |
| 379 | 6771888 | drho_dS(ii,jj,kk) * (S(i+1,j,k)-S(i,j,k)) |
| 380 | ||
| 381 | ! Estimate the vertical density gradients times the grid spacing. | |
| 382 | drdkL = (drho_dT(ii,jj,kk) * (T(i,j,k)-T(i,j,k-1)) + & | |
| 383 | 6771888 | drho_dS(ii,jj,kk) * (S(i,j,k)-S(i,j,k-1))) |
| 384 | drdkR = (drho_dT(ii,jj,kk) * (T(i+1,j,k)-T(i+1,j,k-1)) + & | |
| 385 | 6771888 | drho_dS(ii,jj,kk) * (S(i+1,j,k)-S(i+1,j,k-1))) |
| 386 | 6771888 | if (use_stanley) then |
| 387 | ! Correction to the horizontal density gradient due to nonlinearity in | |
| 388 | ! the EOS rectifying SGS temperature anomalies | |
| 389 | drdiA = drdiA + 0.5 * ((drho_dT_dT_h(ii+1,jj,kk) * tv%varT(i+1,j,K-1)) - & | |
| 390 | 0 | (drho_dT_dT_h(ii,jj,kk) * tv%varT(i,j,K-1)) ) |
| 391 | drdiB = drdiB + 0.5 * ((drho_dT_dT_h(ii+1,jj,kk) * tv%varT(i+1,j,K)) - & | |
| 392 | 0 | (drho_dT_dT_h(ii,jj,kk) * tv%varT(i,j,K)) ) |
| 393 | endif | |
| 394 | else | |
| 395 | 0 | drdiA = 0.0 ; drdiB = 0.0 |
| 396 | 0 | drdkL = GV%Rlay(K)-GV%Rlay(K-1) ; drdkR = drdkL |
| 397 | endif | |
| 398 | ||
| 399 | 6771888 | hg2A = h(i,j,k-1)*h(i+1,j,k-1) + h_neglect2 |
| 400 | 6771888 | hg2B = h(i,j,k)*h(i+1,j,k) + h_neglect2 |
| 401 | 6771888 | hg2L = h(i,j,k-1)*h(i,j,k) + h_neglect2 |
| 402 | 6771888 | hg2R = h(i+1,j,k-1)*h(i+1,j,k) + h_neglect2 |
| 403 | 6771888 | haA = 0.5*(h(i,j,k-1) + h(i+1,j,k-1)) + h_neglect |
| 404 | 6771888 | haB = 0.5*(h(i,j,k) + h(i+1,j,k)) + h_neglect |
| 405 | 6771888 | haL = 0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect |
| 406 | 6771888 | haR = 0.5*(h(i+1,j,k-1) + h(i+1,j,k)) + h_neglect |
| 407 | 6771888 | if (GV%Boussinesq) then |
| 408 | 6771888 | dzaL = haL * GV%H_to_Z ; dzaR = haR * GV%H_to_Z |
| 409 | else | |
| 410 | 0 | dzaL = 0.5*(e(i,j,K-1) - e(i,j,K+1)) + dz_neglect |
| 411 | 0 | dzaR = 0.5*(e(i+1,j,K-1) - e(i+1,j,K+1)) + dz_neglect |
| 412 | endif | |
| 413 | 6771888 | if (present(dzu)) dzu(I,j,K) = 0.5*( dzaL + dzaR ) |
| 414 | ! Use the harmonic mean thicknesses to weight the horizontal gradients. | |
| 415 | ! These unnormalized weights have been rearranged to minimize divisions. | |
| 416 | 6771888 | wtA = hg2A*haB ; wtB = hg2B*haA |
| 417 | 6771888 | wtL = hg2L*(haR*dzaR) ; wtR = hg2R*(haL*dzaL) |
| 418 | ||
| 419 | 6771888 | drdz = ((wtL * drdkL) + (wtR * drdkR)) / ((dzaL*wtL) + (dzaR*wtR)) |
| 420 | ! The expression for drdz above is mathematically equivalent to: | |
| 421 | ! drdz = ((hg2L/haL) * drdkL/dzaL + (hg2R/haR) * drdkR/dzaR) / & | |
| 422 | ! ((hg2L/haL) + (hg2R/haR)) | |
| 423 | ! which is an estimate of the gradient of density across geopotentials. | |
| 424 | 6771888 | if (present_N2_u) then |
| 425 | 6771888 | if (OBC_friendly) then ; if (OBC%segnum_u(I,j) /= 0) then |
| 426 | 0 | if (OBC%segnum_u(I,j) > 0) then ! OBC_DIRECTION_E |
| 427 | 0 | drdz = drdkL / dzaL ! Note that drdz is not used for slopes at OBC faces. |
| 428 | 0 | if (use_EOS .and. allocated(tv%SpV_avg)) & |
| 429 | 0 | GxSpV_uvh(ii,jj,kk) = GV%g_Earth * 0.5 * (tv%SpV_avg(i,j,k) + tv%SpV_avg(i,j,k-1)) |
| 430 | 0 | elseif (OBC%segnum_u(I,j) < 0) then ! OBC_DIRECTION_W |
| 431 | 0 | drdz = drdkR / dzaR |
| 432 | 0 | if (use_EOS .and. allocated(tv%SpV_avg)) & |
| 433 | 0 | GxSpV_uvh(ii,jj,kk) = GV%g_Earth * 0.5 * (tv%SpV_avg(i+1,j,k) + tv%SpV_avg(i+1,j,k-1)) |
| 434 | endif | |
| 435 | endif ; endif | |
| 436 | 6771888 | N2_u(I,j,k) = GxSpV_uvh(ii,jj,kk) * drdz * G%mask2dCu(I,j) ! Square of buoyancy freq. [L2 Z-2 T-2 ~> s-2] |
| 437 | endif | |
| 438 | ||
| 439 | 6771888 | if (use_EOS) then |
| 440 | drdx = ((wtA * drdiA + wtB * drdiB) / (wtA + wtB) - & | |
| 441 | 6771888 | drdz * (e(i,j,K)-e(i+1,j,K))) * G%IdxCu(I,j) |
| 442 | ||
| 443 | ! This estimate of slope is accurate for small slopes, but bounded | |
| 444 | ! to be between -1 and 1. | |
| 445 | 6771888 | mag_grad2 = (US%Z_to_L*drdx)**2 + drdz**2 |
| 446 | 6771888 | if (mag_grad2 > 0.0) then |
| 447 | 6468111 | slope = drdx / sqrt(mag_grad2) |
| 448 | else ! Just in case mag_grad2 = 0 ever. | |
| 449 | 303777 | slope = 0.0 |
| 450 | endif | |
| 451 | else ! With .not.use_EOS, the layers are constant density. | |
| 452 | 0 | slope = (e(i+1,j,K)-e(i,j,K)) * G%IdxCu(I,j) |
| 453 | endif | |
| 454 | ||
| 455 | 6771888 | if (local_open_u_BC) then |
| 456 | 0 | if (OBC%segnum_u(I,j) /= 0) then |
| 457 | 0 | if (OBC%segment(abs(OBC%segnum_u(I,j)))%open) then |
| 458 | 0 | slope = 0. |
| 459 | ! This and/or the masking code below is to make slopes match inside | |
| 460 | ! land mask. Might not be necessary except for DEBUG output. | |
| 461 | ! if (OBC%segnum_u(I,j) > 0) then ! OBC_DIRECTION_E | |
| 462 | ! slope_x(I+1,j,K) = 0. | |
| 463 | ! else | |
| 464 | ! slope_x(I-1,j,K) = 0. | |
| 465 | ! endif | |
| 466 | endif | |
| 467 | endif | |
| 468 | 0 | slope = slope * max(G%mask2dT(i,j), G%mask2dT(i+1,j)) |
| 469 | endif | |
| 470 | ||
| 471 | 6771888 | slope_x(I,j,K) = slope |
| 472 | 6771888 | if (present(dzSxN)) & |
| 473 | dzSxN(I,j,K) = sqrt( GxSpV_uvh(ii,jj,kk) * max(0., (wtL * ( dzaL * drdkL )) & | |
| 474 | + (wtR * ( dzaR * drdkR ))) / (wtL + wtR) ) & | |
| 475 | 20425776 | * abs(slope) * G%mask2dCu(I,j) |
| 476 | enddo ! end zonal block do concurrent | |
| 477 | enddo ; enddo ; enddo ! end zonal outer block loops | |
| 478 | ||
| 479 | 12 | if (use_stanley) then |
| 480 | 0 | delta_j = njblock - 1 |
| 481 | 0 | je_read = je + 1 |
| 482 | else | |
| 483 | 12 | delta_j = njblock |
| 484 | 12 | je_read = je |
| 485 | endif | |
| 486 | ||
| 487 | ! Blocked meridional loop | |
| 488 | 169620 | do kstart=2,nz, nkblock ; do jstart=js-1,je,delta_j ; do istart=is,ie,niblock |
| 489 | 111888 | iend = min(istart + niblock - 1, ie) |
| 490 | 111888 | jend = min(jstart + delta_j - 1, je) |
| 491 | 111888 | jend_stanley = min(jstart + njblock - 1, je_read) |
| 492 | 111888 | kend = min(kstart + nkblock - 1, nz) |
| 493 | 111888 | EOSdom_block(1,1) = 1 ; EOSdom_block(1,2) = iend - istart + 1 |
| 494 | 111888 | EOSdom_block(2,1) = 1 ; EOSdom_block(2,2) = jend - jstart + 1 |
| 495 | 111888 | EOSdom_block(3,1) = 1 ; EOSdom_block(3,2) = kend - kstart + 1 |
| 496 | ||
| 497 | 111888 | if (use_EOS) then |
| 498 | ! Fill block T_uvh/S_uvh/pres_uvh at v-points | |
| 499 | do concurrent(kk=1:kend-kstart+1, jj=1:jend-jstart+1, ii=1:iend-istart+1) & | |
| 500 | 13762224 | DO_LOCALITY(local(i,j,k)) |
| 501 | 6825168 | i = istart + ii - 1 |
| 502 | 6825168 | j = jstart + jj - 1 |
| 503 | 6825168 | k = kstart + kk - 1 |
| 504 | 6825168 | pres_uvh(ii,jj,kk) = 0.5*(pres(i,j,K) + pres(i,j+1,K)) |
| 505 | 6825168 | T_uvh(ii,jj,kk) = 0.25*((T(i,j,K) + T(i,j+1,K)) + (T(i,j,K-1) + T(i,j+1,K-1))) |
| 506 | 20587392 | S_uvh(ii,jj,kk) = 0.25*((S(i,j,K) + S(i,j+1,K)) + (S(i,j,K-1) + S(i,j+1,K-1))) |
| 507 | enddo | |
| 508 | ||
| 509 | 111888 | if (OBC_friendly) then |
| 510 | !$omp target update from(T_uvh, S_uvh, pres_uvh) | |
| 511 | ||
| 512 | 0 | if (OBC%v_N_OBCs_on_PE) then |
| 513 | 0 | do k = kstart, kend |
| 514 | 0 | kk = k - kstart + 1 |
| 515 | 0 | do j = max(jstart, OBC%Js_v_N_obc), min(jend, OBC%Je_v_N_obc) |
| 516 | 0 | jj = j - jstart + 1 |
| 517 | 0 | do i = max(istart, OBC%is_v_N_obc), min(iend, OBC%ie_v_N_obc) |
| 518 | 0 | ii = i - istart + 1 |
| 519 | 0 | if (OBC%segnum_v(i,j) > 0) then ! OBC_DIRECTION_N |
| 520 | 0 | pres_uvh(ii,jj,kk) = pres(i,j,K) |
| 521 | 0 | T_uvh(ii,jj,kk) = 0.5*(T(i,j,K) + T(i,j,K-1)) |
| 522 | 0 | S_uvh(ii,jj,kk) = 0.5*(S(i,j,K) + S(i,j,K-1)) |
| 523 | endif | |
| 524 | enddo | |
| 525 | enddo | |
| 526 | enddo | |
| 527 | endif | |
| 528 | 0 | if (OBC%v_S_OBCs_on_PE) then |
| 529 | 0 | do k = kstart, kend |
| 530 | 0 | kk = k - kstart + 1 |
| 531 | 0 | do j = max(jstart, OBC%Js_v_S_obc), min(jend, OBC%Je_v_S_obc) |
| 532 | 0 | jj = j - jstart + 1 |
| 533 | 0 | do i = max(istart, OBC%is_v_S_obc), min(iend, OBC%ie_v_S_obc) |
| 534 | 0 | ii = i - istart + 1 |
| 535 | 0 | if (OBC%segnum_v(i,j) < 0) then ! OBC_DIRECTION_S |
| 536 | 0 | pres_uvh(ii,jj,kk) = pres(i,j+1,K) |
| 537 | 0 | T_uvh(ii,jj,kk) = 0.5*(T(i,j+1,K) + T(i,j+1,K-1)) |
| 538 | 0 | S_uvh(ii,jj,kk) = 0.5*(S(i,j+1,K) + S(i,j+1,K-1)) |
| 539 | endif | |
| 540 | enddo | |
| 541 | enddo | |
| 542 | enddo | |
| 543 | endif | |
| 544 | ||
| 545 | !$omp target update to(T_uvh, S_uvh, pres_uvh) | |
| 546 | endif | |
| 547 | ||
| 548 | 111888 | if (present_N2_v .or. present(dzSyN)) then |
| 549 | 111888 | if (allocated(tv%SpV_avg)) then |
| 550 | do concurrent(kk=1:kend-kstart+1, jj=1:jend-jstart+1, ii=1:iend-istart+1) & | |
| 551 | 0 | DO_LOCALITY(local(i,j,k)) |
| 552 | 0 | i = istart + ii - 1 |
| 553 | 0 | j = jstart + jj - 1 |
| 554 | 0 | k = kstart + kk - 1 |
| 555 | GxSpV_uvh(ii,jj,kk) = GV%g_Earth * 0.25 * ((tv%SpV_avg(i,j,K) + tv%SpV_avg(i,j+1,K)) + & | |
| 556 | 0 | (tv%SpV_avg(i,j,K-1) + tv%SpV_avg(i,j+1,K-1))) |
| 557 | enddo | |
| 558 | endif | |
| 559 | endif | |
| 560 | ||
| 561 | call calculate_density_derivs(T_uvh, S_uvh, pres_uvh, & | |
| 562 | drho_dT, drho_dS, & | |
| 563 | 111888 | tv%eqn_of_state, EOSdom_block) |
| 564 | ||
| 565 | 111888 | if (use_stanley) then |
| 566 | ! Reset at h-points for Stanley second derivatives. | |
| 567 | 0 | EOSdom_block_h1(1) = 1 ; EOSdom_block_h1(2) = iend - istart + 1 |
| 568 | do concurrent(kk=1:kend-kstart+1, jj=1:jend_stanley-jstart+1, ii=1:iend-istart+1) & | |
| 569 | 0 | DO_LOCALITY(local(i,j,k)) |
| 570 | 0 | i = istart + ii - 1 |
| 571 | 0 | j = jstart + jj - 1 |
| 572 | 0 | k = kstart + kk - 1 |
| 573 | 0 | pres_uvh(ii,jj,kk) = pres(i,j,K) |
| 574 | 0 | T_uvh(ii,jj,kk) = 0.5*(T(i,j,K) + T(i,j,K-1)) |
| 575 | 0 | S_uvh(ii,jj,kk) = 0.5*(S(i,j,K) + S(i,j,K-1)) |
| 576 | enddo | |
| 577 | ||
| 578 | ! TODO: Make 3D version of this | |
| 579 | 0 | do kk=1,kend-kstart+1 ; do jj=1,jend_stanley-jstart+1 |
| 580 | ! The second line below would correspond to arguments | |
| 581 | ! drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, drho_dT_dP, & | |
| 582 | call calculate_density_second_derivs(T_uvh(:,jj,kk), S_uvh(:,jj,kk), pres_uvh(:,jj,kk), & | |
| 583 | scrap, scrap, drho_dT_dT_h(:,jj,kk), scrap, scrap, & | |
| 584 | 0 | tv%eqn_of_state, dom=EOSdom_block_h1) |
| 585 | enddo ; enddo | |
| 586 | endif ! end use_stanley | |
| 587 | ||
| 588 | endif ! end use_EOS for meridional block | |
| 589 | ||
| 590 | ! Meridional slope compute over the block | |
| 591 | do concurrent(kk=1:kend-kstart+1, jj=1:jend-jstart+1, ii=1:iend-istart+1) & | |
| 592 | 55944 | DO_LOCALITY(local(drdkL, drdkR, drdjA, drdjB, i, J)) |
| 593 | 6825168 | i = istart + ii - 1 |
| 594 | 6825168 | J = jstart + jj - 1 |
| 595 | 6825168 | K = kstart + kk - 1 |
| 596 | ||
| 597 | 6825168 | if (use_EOS) then |
| 598 | ! Estimate the horizontal density gradients along layers. | |
| 599 | drdjA = drho_dT(ii,jj,kk) * (T(i,j+1,k-1)-T(i,j,k-1)) + & | |
| 600 | 6825168 | drho_dS(ii,jj,kk) * (S(i,j+1,k-1)-S(i,j,k-1)) |
| 601 | drdjB = drho_dT(ii,jj,kk) * (T(i,j+1,k)-T(i,j,k)) + & | |
| 602 | 6825168 | drho_dS(ii,jj,kk) * (S(i,j+1,k)-S(i,j,k)) |
| 603 | ||
| 604 | ! Estimate the vertical density gradients times the grid spacing. | |
| 605 | drdkL = (drho_dT(ii,jj,kk) * (T(i,j,K)-T(i,j,K-1)) + & | |
| 606 | 6825168 | drho_dS(ii,jj,kk) * (S(i,j,K)-S(i,j,K-1))) |
| 607 | drdkR = (drho_dT(ii,jj,kk) * (T(i,j+1,K)-T(i,j+1,K-1)) + & | |
| 608 | 6825168 | drho_dS(ii,jj,kk) * (S(i,j+1,K)-S(i,j+1,K-1))) |
| 609 | 6825168 | if (use_stanley) then |
| 610 | ! Correction to the horizontal density gradient due to nonlinearity in | |
| 611 | ! the EOS rectifying SGS temperature anomalies | |
| 612 | drdjA = drdjA + 0.5 * ((drho_dT_dT_h(ii,jj+1,kk) * tv%varT(i,j+1,K-1)) - & | |
| 613 | 0 | (drho_dT_dT_h(ii,jj,kk) * tv%varT(i,j,K-1)) ) |
| 614 | drdjB = drdjB + 0.5 * ((drho_dT_dT_h(ii,jj+1,kk) * tv%varT(i,j+1,K)) - & | |
| 615 | 0 | (drho_dT_dT_h(ii,jj,kk) * tv%varT(i,j,K)) ) |
| 616 | endif | |
| 617 | else | |
| 618 | 0 | drdjA = 0.0 ; drdjB = 0.0 |
| 619 | 0 | drdkL = GV%Rlay(K)-GV%Rlay(K-1) ; drdkR = drdkL |
| 620 | endif | |
| 621 | ||
| 622 | 6825168 | hg2A = h(i,j,k-1)*h(i,j+1,k-1) + h_neglect2 |
| 623 | 6825168 | hg2B = h(i,j,k)*h(i,j+1,k) + h_neglect2 |
| 624 | 6825168 | hg2L = h(i,j,k-1)*h(i,j,k) + h_neglect2 |
| 625 | 6825168 | hg2R = h(i,j+1,k-1)*h(i,j+1,k) + h_neglect2 |
| 626 | 6825168 | haA = 0.5*(h(i,j,k-1) + h(i,j+1,k-1)) + h_neglect |
| 627 | 6825168 | haB = 0.5*(h(i,j,k) + h(i,j+1,k)) + h_neglect |
| 628 | 6825168 | haL = 0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect |
| 629 | 6825168 | haR = 0.5*(h(i,j+1,k-1) + h(i,j+1,k)) + h_neglect |
| 630 | 6825168 | if (GV%Boussinesq) then |
| 631 | 6825168 | dzaL = haL * GV%H_to_Z ; dzaR = haR * GV%H_to_Z |
| 632 | else | |
| 633 | 0 | dzaL = 0.5*(e(i,j,K-1) - e(i,j,K+1)) + dz_neglect |
| 634 | 0 | dzaR = 0.5*(e(i,j+1,K-1) - e(i,j+1,K+1)) + dz_neglect |
| 635 | endif | |
| 636 | 6825168 | if (present(dzv)) dzv(i,J,K) = 0.5*( dzaL + dzaR ) |
| 637 | ! Use the harmonic mean thicknesses to weight the horizontal gradients. | |
| 638 | ! These unnormalized weights have been rearranged to minimize divisions. | |
| 639 | 6825168 | wtA = hg2A*haB ; wtB = hg2B*haA |
| 640 | 6825168 | wtL = hg2L*(haR*dzaR) ; wtR = hg2R*(haL*dzaL) |
| 641 | ||
| 642 | 6825168 | drdz = ((wtL * drdkL) + (wtR * drdkR)) / ((dzaL*wtL) + (dzaR*wtR)) |
| 643 | ! The expression for drdz above is mathematically equivalent to: | |
| 644 | ! drdz = ((hg2L/haL) * drdkL/dzaL + (hg2R/haR) * drdkR/dzaR) / & | |
| 645 | ! ((hg2L/haL) + (hg2R/haR)) | |
| 646 | ! which is an estimate of the gradient of density across geopotentials. | |
| 647 | 6825168 | if (present_N2_v) then |
| 648 | 6825168 | if (OBC_friendly) then ; if (OBC%segnum_v(i,J) /= 0) then |
| 649 | 0 | if (OBC%segnum_v(i,J) > 0) then ! OBC_DIRECTION_N |
| 650 | 0 | drdz = drdkL / dzaL ! Note that drdz is not used for slopes at OBC faces. |
| 651 | 0 | if (use_EOS .and. allocated(tv%SpV_avg)) & |
| 652 | 0 | GxSpV_uvh(ii,jj,kk) = GV%g_Earth * 0.5 * (tv%SpV_avg(i,j,k) + tv%SpV_avg(i,j,k-1)) |
| 653 | 0 | elseif (OBC%segnum_v(i,J) < 0) then ! OBC_DIRECTION_S |
| 654 | 0 | drdz = drdkL / dzaL |
| 655 | 0 | if (use_EOS .and. allocated(tv%SpV_avg)) & |
| 656 | 0 | GxSpV_uvh(ii,jj,kk) = GV%g_Earth * 0.5 * (tv%SpV_avg(i,j+1,k) + tv%SpV_avg(i,j+1,k-1)) |
| 657 | endif | |
| 658 | endif ; endif | |
| 659 | 6825168 | N2_v(i,J,K) = GxSpV_uvh(ii,jj,kk) * drdz * G%mask2dCv(i,J) ! Square of buoyancy freq. [L2 Z-2 T-2 ~> s-2] |
| 660 | endif | |
| 661 | ||
| 662 | 6825168 | if (use_EOS) then |
| 663 | drdy = ((wtA * drdjA + wtB * drdjB) / (wtA + wtB) - & | |
| 664 | 6825168 | drdz * (e(i,j,K)-e(i,j+1,K))) * G%IdyCv(i,J) |
| 665 | ||
| 666 | ! This estimate of slope is accurate for small slopes, but bounded | |
| 667 | ! to be between -1 and 1. | |
| 668 | 6825168 | mag_grad2 = (US%Z_to_L*drdy)**2 + drdz**2 |
| 669 | 6825168 | if (mag_grad2 > 0.0) then |
| 670 | 6601099 | slope = drdy / sqrt(mag_grad2) |
| 671 | else ! Just in case mag_grad2 = 0 ever. | |
| 672 | 224069 | slope = 0.0 |
| 673 | endif | |
| 674 | else ! With .not.use_EOS, the layers are constant density. | |
| 675 | 0 | slope = (e(i,j+1,K)-e(i,j,K)) * G%IdyCv(i,J) |
| 676 | endif | |
| 677 | ||
| 678 | 6825168 | if (local_open_v_BC) then |
| 679 | 0 | if (OBC%segnum_v(i,J) /= 0) then |
| 680 | 0 | if (OBC%segment(abs(OBC%segnum_v(i,J)))%open) then |
| 681 | 0 | slope = 0. |
| 682 | ! This and/or the masking code below is to make slopes match inside | |
| 683 | ! land mask. Might not be necessary except for DEBUG output. | |
| 684 | ! if (OBC%segnum_v(i,J)) > 0) then ! OBC_DIRECTION_N | |
| 685 | ! slope_y(i,J+1,K) = 0. | |
| 686 | ! else | |
| 687 | ! slope_y(i,J-1,K) = 0. | |
| 688 | ! endif | |
| 689 | endif | |
| 690 | endif | |
| 691 | 0 | slope = slope * max(G%mask2dT(i,j), G%mask2dT(i,j+1)) |
| 692 | endif | |
| 693 | 6825168 | slope_y(i,J,K) = slope |
| 694 | 6825168 | if (present(dzSyN)) & |
| 695 | dzSyN(i,J,K) = sqrt( GxSpV_uvh(ii,jj,kk) * max(0., (wtL * ( dzaL * drdkL )) & | |
| 696 | + (wtR * ( dzaR * drdkR ))) / (wtL + wtR) ) & | |
| 697 | 20587392 | * abs(slope) * G%mask2dCv(i,J) |
| 698 | enddo ! end meridional block do concurrent | |
| 699 | enddo ; enddo ; enddo ! end meridional outer block loops | |
| 700 | ||
| 701 | ! Delete all block and field arrays from device | |
| 702 | !$omp target exit data map(delete: T, S, pres, T_uvh, S_uvh, pres_uvh, GxSpV_uvh, & | |
| 703 | !$omp scrap, drho_dT, drho_dS, drho_dT_dT_h) | |
| 704 | ||
| 705 | 36 | end subroutine calc_isoneutral_slopes |
| 706 | ||
| 707 | !> Returns tracer arrays (nominally T and S) with massless layers filled with | |
| 708 | !! sensible values, by diffusing vertically with a small but constant diffusivity. | |
| 709 | 36 | subroutine vert_fill_TS(h, T_in, S_in, kappa_dt, T_f, S_f, G, GV, US, halo_here, larger_h_denom) |
| 710 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure | |
| 711 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure | |
| 712 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 713 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] | |
| 714 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: T_in !< Input temperature [C ~> degC] | |
| 715 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: S_in !< Input salinity [S ~> ppt] | |
| 716 | real, intent(in) :: kappa_dt !< A vertical diffusivity to use for smoothing | |
| 717 | !! times a smoothing timescale [H Z ~> m2 or kg m-1] | |
| 718 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T_f !< Filled temperature [C ~> degC] | |
| 719 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S_f !< Filled salinity [S ~> ppt] | |
| 720 | integer, optional, intent(in) :: halo_here !< Number of halo points to work on, | |
| 721 | !! 0 by default | |
| 722 | logical, optional, intent(in) :: larger_h_denom !< Present and true, add a large | |
| 723 | !! enough minimal thickness in the denominator of | |
| 724 | !! the flux calculations so that the fluxes are | |
| 725 | !! never so large as eliminate the transmission | |
| 726 | !! of information across groups of massless layers. | |
| 727 | ! Local variables | |
| 728 | 72 | real :: ent(SZI_(G),SZK_(GV)+1) ! The diffusive entrainment (kappa*dt)/dz |
| 729 | ! between layers in a timestep [H ~> m or kg m-2]. | |
| 730 | 72 | real :: b1(SZI_(G)) ! A variable used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1] |
| 731 | 72 | real :: d1(SZI_(G)) ! A variable used by the tridiagonal solver [nondim], d1 = 1 - c1. |
| 732 | 72 | real :: c1(SZI_(G),SZK_(GV)) ! A variable used by the tridiagonal solver [nondim]. |
| 733 | real :: kap_dt_x2 ! The 2*kappa_dt converted to H units [H2 ~> m2 or kg2 m-4]. | |
| 734 | real :: h_neglect ! A negligible thickness [H ~> m or kg m-2], to allow for zero thicknesses. | |
| 735 | real :: h0 ! A negligible thickness to allow for zero thickness layers without | |
| 736 | ! completely decoupling groups of layers [H ~> m or kg m-2]. | |
| 737 | ! Often 0 < h_neglect << h0. | |
| 738 | real :: h_tr ! h_tr is h at tracer points with a tiny thickness | |
| 739 | ! added to ensure positive definiteness [H ~> m or kg m-2]. | |
| 740 | integer :: i, j, k, is, ie, js, je, nz, halo | |
| 741 | ||
| 742 | 36 | halo=0 ; if (present(halo_here)) halo = max(halo_here,0) |
| 743 | ||
| 744 | 36 | is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo ; nz = GV%ke |
| 745 | ||
| 746 | 36 | h_neglect = GV%H_subroundoff |
| 747 | ! The use of the fixed rescaling factor in the next line avoids an extra call to thickness_to_dz() | |
| 748 | ! and the use of an extra 3-d array of vertical distnaces across layers (dz). This would be more | |
| 749 | ! physically consistent, but it would also be more expensive, and given that this routine applies | |
| 750 | ! a small (but arbitrary) amount of mixing to clean up the properties of nearly massless layers, | |
| 751 | ! the added expense is hard to justify. | |
| 752 | 36 | kap_dt_x2 = (2.0*kappa_dt) * (US%Z_to_m*GV%m_to_H) ! Usually the latter term is GV%Z_to_H. |
| 753 | 36 | h0 = h_neglect |
| 754 | 36 | if (present(larger_h_denom)) then |
| 755 | 24 | if (larger_h_denom) h0 = 1.0e-16*sqrt(0.5*kap_dt_x2) |
| 756 | endif | |
| 757 | ||
| 758 | !$omp target enter data map(alloc: ent, b1, d1, c1) | |
| 759 | ||
| 760 | 36 | if (kap_dt_x2 <= 0.0) then |
| 761 | 0 | do concurrent (k=1:nz , j=js:je , i=is:ie) |
| 762 | 0 | T_f(i,j,k) = T_in(i,j,k) ; S_f(i,j,k) = S_in(i,j,k) |
| 763 | enddo | |
| 764 | else | |
| 765 | !$omp target teams loop private( ent, b1, d1, c1, h_tr ) | |
| 766 | 2268 | do j=js,je |
| 767 | 2232 | do concurrent( i=is:ie ) |
| 768 | 272400 | ent(i,2) = kap_dt_x2 / ((h(i,j,1)+h(i,j,2)) + h0) |
| 769 | 272400 | h_tr = h(i,j,1) + h_neglect |
| 770 | 272400 | b1(i) = 1.0 / (h_tr + ent(i,2)) |
| 771 | 272400 | d1(i) = b1(i) * h_tr |
| 772 | 272400 | T_f(i,j,1) = (b1(i)*h_tr)*T_in(i,j,1) |
| 773 | 274632 | S_f(i,j,1) = (b1(i)*h_tr)*S_in(i,j,1) |
| 774 | enddo | |
| 775 | 165168 | do k=2,nz-1 ; do concurrent( i=is:ie ) |
| 776 | 19885200 | ent(i,K+1) = kap_dt_x2 / ((h(i,j,k)+h(i,j,k+1)) + h0) |
| 777 | 19885200 | h_tr = h(i,j,k) + h_neglect |
| 778 | 19885200 | c1(i,k) = ent(i,K) * b1(i) |
| 779 | 19885200 | b1(i) = 1.0 / ((h_tr + d1(i)*ent(i,K)) + ent(i,K+1)) |
| 780 | 19885200 | d1(i) = b1(i) * (h_tr + d1(i)*ent(i,K)) |
| 781 | 19885200 | T_f(i,j,k) = b1(i) * (h_tr*T_in(i,j,k) + ent(i,K)*T_f(i,j,k-1)) |
| 782 | 20048136 | S_f(i,j,k) = b1(i) * (h_tr*S_in(i,j,k) + ent(i,K)*S_f(i,j,k-1)) |
| 783 | enddo ; enddo | |
| 784 | 2232 | do concurrent( i=is:ie ) |
| 785 | 272400 | c1(i,nz) = ent(i,nz) * b1(i) |
| 786 | 272400 | h_tr = h(i,j,nz) + h_neglect |
| 787 | 272400 | b1(i) = 1.0 / (h_tr + d1(i)*ent(i,nz)) |
| 788 | 272400 | T_f(i,j,nz) = b1(i) * (h_tr*T_in(i,j,nz) + ent(i,nz)*T_f(i,j,nz-1)) |
| 789 | 274632 | S_f(i,j,nz) = b1(i) * (h_tr*S_in(i,j,nz) + ent(i,nz)*S_f(i,j,nz-1)) |
| 790 | enddo | |
| 791 | 167436 | do k=nz-1,1,-1 ; do concurrent( i=is:ie ) |
| 792 | 20157600 | T_f(i,j,k) = T_f(i,j,k) + c1(i,k+1)*T_f(i,j,k+1) |
| 793 | 20322768 | S_f(i,j,k) = S_f(i,j,k) + c1(i,k+1)*S_f(i,j,k+1) |
| 794 | enddo ; enddo | |
| 795 | enddo | |
| 796 | endif | |
| 797 | ||
| 798 | !$omp target exit data map(delete: ent, b1, d1, c1) | |
| 799 | ||
| 800 | 36 | end subroutine vert_fill_TS |
| 801 | ||
| 802 | 0 | end module MOM_isopycnal_slopes |