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 | !> Functions for calculating interface heights, including free surface height. | |
| 6 | module MOM_interface_heights | |
| 7 | ||
| 8 | use MOM_density_integrals, only : int_specific_vol_dp, avg_specific_vol, int_density_dz | |
| 9 | use MOM_debugging, only : hchksum | |
| 10 | use MOM_error_handler, only : MOM_error, FATAL | |
| 11 | use MOM_EOS, only : calculate_density, average_specific_vol, EOS_type, EOS_domain | |
| 12 | use MOM_file_parser, only : log_version | |
| 13 | use MOM_grid, only : ocean_grid_type | |
| 14 | use MOM_unit_scaling, only : unit_scale_type | |
| 15 | use MOM_variables, only : thermo_var_ptrs | |
| 16 | use MOM_verticalGrid, only : verticalGrid_type | |
| 17 | ||
| 18 | implicit none ; private | |
| 19 | ||
| 20 | #include <MOM_memory.h> | |
| 21 | ||
| 22 | public find_eta, find_dz_for_eta, dz_to_thickness, thickness_to_dz, dz_to_thickness_simple | |
| 23 | public calc_derived_thermo | |
| 24 | public convert_MLD_to_ML_thickness | |
| 25 | public find_rho_bottom, find_col_avg_SpV, find_col_mass | |
| 26 | ||
| 27 | !> Calculates the heights of the free surface or all interfaces from layer thicknesses. | |
| 28 | interface find_eta | |
| 29 | module procedure find_eta_2d, find_eta_3d | |
| 30 | end interface find_eta | |
| 31 | ||
| 32 | !> Calculates layer thickness in thickness units from geometric distance between the | |
| 33 | !! interfaces around that layer in height units. | |
| 34 | interface dz_to_thickness | |
| 35 | module procedure dz_to_thickness_tv, dz_to_thickness_EoS | |
| 36 | end interface dz_to_thickness | |
| 37 | ||
| 38 | !> Converts layer thickness in thickness units into the vertical distance between the | |
| 39 | !! interfaces around a layer in height units. | |
| 40 | interface thickness_to_dz | |
| 41 | module procedure thickness_to_dz_3d, thickness_to_dz_jslice | |
| 42 | end interface thickness_to_dz | |
| 43 | ||
| 44 | contains | |
| 45 | ||
| 46 | !> Calculates the change in height across layers, using the appropriate form for | |
| 47 | !! consistency with the calculation of the pressure gradient forces. | |
| 48 | 0 | subroutine find_dz_for_eta(h, tv, G, GV, US, dz_lay, halo_size) |
| 49 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. | |
| 50 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. | |
| 51 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 52 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] | |
| 53 | type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various | |
| 54 | !! thermodynamic variables. | |
| 55 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: dz_lay !< Height change across layers [Z ~> m] | |
| 56 | integer, optional, intent(in) :: halo_size !< width of halo points on | |
| 57 | !! which to calculate eta. | |
| 58 | ||
| 59 | ! Local variables | |
| 60 | 0 | real :: p(SZI_(G),SZJ_(G),SZK_(GV)+1) ! Hydrostatic pressure at each interface [R L2 T-2 ~> Pa] |
| 61 | 0 | real :: dz_geo(SZI_(G),SZJ_(G)) ! The change in geopotential height across a layer [L2 T-2 ~> m2 s-2] |
| 62 | real :: SpV_lay_conv(SZK_(GV)) ! The prescribed layer specific volume times a conversion factor from | |
| 63 | ! the units of thickness to layer mass [Z H-1 ~> nondim or m3 kg-1] | |
| 64 | real :: I_gEarth ! The inverse of the gravitational acceleration times the | |
| 65 | ! rescaling factor derived from eta_to_m [T2 Z L-2 ~> s2 m-1] | |
| 66 | integer :: i, j, k, isv, iev, jsv, jev, nz, halo | |
| 67 | ||
| 68 | 0 | halo = 0 ; if (present(halo_size)) halo = max(0,halo_size) |
| 69 | ||
| 70 | 0 | isv = G%isc-halo ; iev = G%iec+halo ; jsv = G%jsc-halo ; jev = G%jec+halo |
| 71 | 0 | nz = GV%ke |
| 72 | ||
| 73 | 0 | if ((isv<G%isd) .or. (iev>G%ied) .or. (jsv<G%jsd) .or. (jev>G%jed)) & |
| 74 | 0 | call MOM_error(FATAL,"find_dz_for_eta called with an overly large halo_size.") |
| 75 | ||
| 76 | 0 | if (GV%Boussinesq) then |
| 77 | 0 | do k=1,nz ; do j=jsv,jev ; do i=isv,iev |
| 78 | 0 | dz_lay(i,j,K) = h(i,j,k)*GV%H_to_Z |
| 79 | enddo ; enddo ; enddo | |
| 80 | 0 | elseif (associated(tv%eqn_of_state)) then |
| 81 | 0 | I_gEarth = 1.0 / GV%g_Earth |
| 82 | !$OMP parallel do default(shared) | |
| 83 | 0 | do j=jsv,jev |
| 84 | 0 | if (associated(tv%p_surf)) then |
| 85 | 0 | do i=isv,iev ; p(i,j,1) = tv%p_surf(i,j) ; enddo |
| 86 | else | |
| 87 | 0 | do i=isv,iev ; p(i,j,1) = 0.0 ; enddo |
| 88 | endif | |
| 89 | 0 | do k=1,nz ; do i=isv,iev |
| 90 | 0 | p(i,j,K+1) = p(i,j,K) + GV%g_Earth*GV%H_to_RZ*h(i,j,k) |
| 91 | enddo ; enddo | |
| 92 | enddo | |
| 93 | !$OMP parallel do default(shared) private(dz_geo) | |
| 94 | 0 | do k=1,nz |
| 95 | call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p(:,:,K), p(:,:,K+1), & | |
| 96 | 0 | 0.0, G%HI, tv%eqn_of_state, US, dz_geo, halo_size=halo) |
| 97 | 0 | do j=jsv,jev ; do i=isv,iev |
| 98 | 0 | dz_lay(i,j,K) = I_gEarth * dz_geo(i,j) |
| 99 | enddo ; enddo | |
| 100 | enddo | |
| 101 | else ! non-Boussinesq but with no equation of state | |
| 102 | 0 | do k=1,nz ; do j=jsv,jev ; do i=isv,iev |
| 103 | 0 | dz_lay(i,j,K) = GV%H_to_RZ*h(i,j,k) / GV%Rlay(k) |
| 104 | enddo ; enddo ; enddo | |
| 105 | ! This would be faster but could change answers. | |
| 106 | ! do k=1,nz ; SpV_lay_conv(k) = GV%H_to_RZ / GV%Rlay(k) ; enddo | |
| 107 | ! do k=1,nz ; do j=jsv,jev ; do i=isv,iev | |
| 108 | ! dz_lay(i,j,K) = h(i,j,k) * SpV_lay_conv(k) | |
| 109 | ! enddo ; enddo ; enddo | |
| 110 | endif | |
| 111 | ||
| 112 | ! To find eta, do the following: | |
| 113 | ! do j=jsv,jev ; do i=isv,iev ; eta(i,j,nz+1) = -(G%bathyT(i,j) + dZ_ref) ; enddo ; enddo | |
| 114 | ! do k=nz,1,-1 ; do j=jsv,jev ; do i=isv,iev | |
| 115 | ! eta(i,j,K) = eta(i,j,K+1) + dz_lay(i,j,K) | |
| 116 | ! enddo ; enddo ; enddo | |
| 117 | ||
| 118 | 0 | end subroutine find_dz_for_eta |
| 119 | ||
| 120 | !> Calculates the heights of all interfaces between layers, using the appropriate | |
| 121 | !! form for consistency with the calculation of the pressure gradient forces. | |
| 122 | !! Additionally, these height may be dilated for consistency with the | |
| 123 | !! corresponding time-average quantity from the barotropic calculation. | |
| 124 | 0 | subroutine find_eta_3d(h, tv, G, GV, US, eta, eta_bt, halo_size, dZref) |
| 125 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. | |
| 126 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. | |
| 127 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 128 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] | |
| 129 | type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various | |
| 130 | !! thermodynamic variables. | |
| 131 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(out) :: eta !< layer interface heights [Z ~> m] | |
| 132 | real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: eta_bt !< optional barotropic variable | |
| 133 | !! that gives the "correct" free surface height (Boussinesq) or total water | |
| 134 | !! column mass per unit area (non-Boussinesq). This is used to dilate the layer | |
| 135 | !! thicknesses when calculating interface heights [H ~> m or kg m-2]. | |
| 136 | !! In Boussinesq mode, eta_bt and G%bathyT use the same reference height. | |
| 137 | integer, optional, intent(in) :: halo_size !< width of halo points on | |
| 138 | !! which to calculate eta. | |
| 139 | real, optional, intent(in) :: dZref !< The difference in the | |
| 140 | !! reference height between G%bathyT and eta [Z ~> m]. The default is 0. | |
| 141 | ||
| 142 | ! Local variables | |
| 143 | 50 | real :: dz_lay(SZI_(G),SZJ_(G),SZK_(GV)) ! The change in height across a layer [Z ~> m] |
| 144 | 50 | real :: dilate(SZI_(G)) ! A non-dimensional dilation factor [nondim] |
| 145 | 50 | real :: htot(SZI_(G)) ! total thickness [H ~> m or kg m-2] |
| 146 | real :: dZ_ref ! The difference in the reference height between G%bathyT and eta [Z ~> m]. | |
| 147 | ! dZ_ref is 0 unless the optional argument dZref is present. | |
| 148 | integer :: i, j, k, isv, iev, jsv, jev, nz, halo | |
| 149 | ||
| 150 | 25 | halo = 0 ; if (present(halo_size)) halo = max(0,halo_size) |
| 151 | ||
| 152 | 25 | isv = G%isc-halo ; iev = G%iec+halo ; jsv = G%jsc-halo ; jev = G%jec+halo |
| 153 | 25 | nz = GV%ke |
| 154 | ||
| 155 | 25 | if ((isv<G%isd) .or. (iev>G%ied) .or. (jsv<G%jsd) .or. (jev>G%jed)) & |
| 156 | 0 | call MOM_error(FATAL,"find_eta called with an overly large halo_size.") |
| 157 | ||
| 158 | 25 | dZ_ref = 0.0 ; if (present(dZref)) dZ_ref = dZref |
| 159 | ||
| 160 | 25 | if (GV%Boussinesq) then |
| 161 | 25 | do concurrent (j=jsv:jev, i=isv:iev) |
| 162 | 196297 | eta(i,j,nz+1) = -(G%bathyT(i,j) + dZ_ref) |
| 163 | enddo | |
| 164 | ||
| 165 | 25 | do concurrent (j=jsv:jev, i=isv:iev) |
| 166 | 14686297 | do k=nz,1,-1 |
| 167 | 14683200 | eta(i,j,K) = eta(i,j,K+1) + h(i,j,k)*GV%H_to_Z |
| 168 | enddo | |
| 169 | enddo | |
| 170 | ||
| 171 | 25 | if (present(eta_bt)) then |
| 172 | ! Dilate the water column to agree with the free surface height | |
| 173 | ! that is used for the dynamics. | |
| 174 | !$omp target update from(eta) | |
| 175 | 0 | do j=jsv,jev |
| 176 | 0 | do i=isv,iev |
| 177 | dilate(i) = (eta_bt(i,j)*GV%H_to_Z + G%bathyT(i,j)) / & | |
| 178 | 0 | (eta(i,j,1) + (G%bathyT(i,j) + dZ_ref)) |
| 179 | enddo | |
| 180 | 0 | do k=1,nz ; do i=isv,iev |
| 181 | eta(i,j,K) = dilate(i) * (eta(i,j,K) + (G%bathyT(i,j) + dZ_ref)) - & | |
| 182 | 0 | (G%bathyT(i,j) + dZ_ref) |
| 183 | enddo ; enddo | |
| 184 | enddo | |
| 185 | !$omp target update to(eta) | |
| 186 | endif | |
| 187 | else | |
| 188 | !$omp target update from(eta) | |
| 189 | 0 | call find_dz_for_eta(h, tv, G, GV, US, dz_lay, halo_size) |
| 190 | ||
| 191 | 0 | do j=jsv,jev |
| 192 | 0 | do i=isv,iev ; eta(i,j,nz+1) = -(G%bathyT(i,j) + dZ_ref) ; enddo |
| 193 | 0 | do k=nz,1,-1 ; do i=isv,iev |
| 194 | 0 | eta(i,j,K) = eta(i,j,K+1) + dz_lay(i,j,k) |
| 195 | enddo ; enddo | |
| 196 | enddo | |
| 197 | ||
| 198 | 0 | if (present(eta_bt)) then |
| 199 | ! Dilate the water column to agree with the free surface height | |
| 200 | ! from the time-averaged barotropic solution. | |
| 201 | 0 | do j=jsv,jev |
| 202 | 0 | do i=isv,iev ; htot(i) = GV%H_subroundoff ; enddo |
| 203 | 0 | do k=1,nz ; do i=isv,iev ; htot(i) = htot(i) + h(i,j,k) ; enddo ; enddo |
| 204 | 0 | do i=isv,iev ; dilate(i) = eta_bt(i,j) / htot(i) ; enddo |
| 205 | 0 | do k=1,nz ; do i=isv,iev |
| 206 | eta(i,j,K) = dilate(i) * (eta(i,j,K) + (G%bathyT(i,j) + dZ_ref)) - & | |
| 207 | 0 | (G%bathyT(i,j) + dZ_ref) |
| 208 | enddo ; enddo | |
| 209 | enddo | |
| 210 | endif | |
| 211 | !$omp target update to(eta) | |
| 212 | endif | |
| 213 | 25 | end subroutine find_eta_3d |
| 214 | ||
| 215 | !> Calculates the free surface height, using the appropriate form for consistency | |
| 216 | !! with the calculation of the pressure gradient forces. Additionally, the sea | |
| 217 | !! surface height may be adjusted for consistency with the corresponding | |
| 218 | !! time-average quantity from the barotropic calculation. | |
| 219 | 25 | subroutine find_eta_2d(h, tv, G, GV, US, eta, eta_bt, halo_size, dZref) |
| 220 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. | |
| 221 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. | |
| 222 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 223 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] | |
| 224 | type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various | |
| 225 | !! thermodynamic variables. | |
| 226 | real, dimension(SZI_(G),SZJ_(G)), intent(out) :: eta !< free surface height relative to | |
| 227 | !! mean sea level (z=0) often [Z ~> m]. | |
| 228 | real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: eta_bt !< optional barotropic | |
| 229 | !! variable that gives the "correct" free surface height (Boussinesq) or total | |
| 230 | !! water column mass per unit area (non-Boussinesq) [H ~> m or kg m-2]. | |
| 231 | !! In Boussinesq mode, eta_bt and G%bathyT use the same reference height. | |
| 232 | integer, optional, intent(in) :: halo_size !< width of halo points on | |
| 233 | !! which to calculate eta. | |
| 234 | real, optional, intent(in) :: dZref !< The difference in the | |
| 235 | !! reference height between G%bathyT and eta [Z ~> m]. The default is 0. | |
| 236 | ||
| 237 | ! Local variables | |
| 238 | 50 | real :: dz_lay(SZI_(G),SZJ_(G),SZK_(GV)) ! The change in height across a layer [Z ~> m] |
| 239 | 50 | real :: htot(SZI_(G)) ! The sum of all layers' thicknesses [H ~> m or kg m-2]. |
| 240 | real :: dZ_ref ! The difference in the reference height between G%bathyT and eta [Z ~> m]. | |
| 241 | ! dZ_ref is 0 unless the optional argument dZref is present. | |
| 242 | integer :: i, j, k, is, ie, js, je, nz, halo | |
| 243 | ||
| 244 | 25 | halo = 0 ; if (present(halo_size)) halo = max(0,halo_size) |
| 245 | 25 | is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo |
| 246 | 25 | nz = GV%ke |
| 247 | ||
| 248 | 25 | dZ_ref = 0.0 ; if (present(dZref)) dZ_ref = dZref |
| 249 | ||
| 250 | 25 | if (GV%Boussinesq) then |
| 251 | 25 | if (present(eta_bt)) then |
| 252 | 25 | do concurrent (j=js:je, i=is:ie) |
| 253 | 183025 | eta(i,j) = GV%H_to_Z*eta_bt(i,j) - dZ_ref |
| 254 | enddo | |
| 255 | else | |
| 256 | 0 | do concurrent (j=js:je) |
| 257 | 0 | do concurrent (i=is:ie) |
| 258 | 0 | eta(i,j) = -G%bathyT(i,j) + dZ_ref |
| 259 | enddo | |
| 260 | ||
| 261 | 0 | do k=1,nz ; do concurrent (i=is:ie) |
| 262 | 0 | eta(i,j) = eta(i,j) + h(i,j,k)*GV%H_to_Z |
| 263 | enddo ; enddo | |
| 264 | enddo | |
| 265 | endif | |
| 266 | else | |
| 267 | !$omp target update from(eta) | |
| 268 | 0 | call find_dz_for_eta(h, tv, G, GV, US, dz_lay, halo_size) |
| 269 | ||
| 270 | 0 | do j=js,je |
| 271 | 0 | do i=is,ie ; eta(i,j) = -(G%bathyT(i,j) + dZ_ref) ; enddo |
| 272 | 0 | do k=1,nz ; do i=is,ie |
| 273 | 0 | eta(i,j) = eta(i,j) + dz_lay(i,j,k) |
| 274 | enddo ; enddo | |
| 275 | enddo | |
| 276 | ||
| 277 | 0 | if (present(eta_bt)) then |
| 278 | ! Dilate the water column to agree with the time-averaged column | |
| 279 | ! mass from the barotropic solution. | |
| 280 | 0 | do j=js,je |
| 281 | 0 | do i=is,ie ; htot(i) = GV%H_subroundoff ; enddo |
| 282 | 0 | do k=1,nz ; do i=is,ie ; htot(i) = htot(i) + h(i,j,k) ; enddo ; enddo |
| 283 | 0 | do i=is,ie |
| 284 | eta(i,j) = (eta_bt(i,j) / htot(i)) * (eta(i,j) + (G%bathyT(i,j) + dZ_ref)) - & | |
| 285 | 0 | (G%bathyT(i,j) + dZ_ref) |
| 286 | enddo | |
| 287 | enddo | |
| 288 | !$omp target update to(eta) | |
| 289 | endif | |
| 290 | endif | |
| 291 | 25 | end subroutine find_eta_2d |
| 292 | ||
| 293 | ||
| 294 | !> Calculate derived thermodynamic quantities for re-use later. | |
| 295 | 0 | subroutine calc_derived_thermo(tv, h, G, GV, US, halo, debug) |
| 296 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure | |
| 297 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure | |
| 298 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 299 | type(thermo_var_ptrs), intent(inout) :: tv !< A structure pointing to various | |
| 300 | !! thermodynamic variables, some of | |
| 301 | !! which will be set here. | |
| 302 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & | |
| 303 | intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]. | |
| 304 | integer, optional, intent(in) :: halo !< Width of halo within which to | |
| 305 | !! calculate thicknesses | |
| 306 | logical, optional, intent(in) :: debug !< If present and true, write debugging checksums | |
| 307 | ! Local variables | |
| 308 | 0 | real, dimension(SZI_(G),SZJ_(G)) :: p_t ! Hydrostatic pressure atop a layer [R L2 T-2 ~> Pa] |
| 309 | 0 | real, dimension(SZI_(G),SZJ_(G)) :: dp ! Pressure change across a layer [R L2 T-2 ~> Pa] |
| 310 | 0 | real, dimension(SZK_(GV)) :: SpV_lay ! The specific volume of each layer when no equation of |
| 311 | ! state is used [R-1 ~> m3 kg-1] | |
| 312 | logical :: do_debug ! If true, write checksums for debugging. | |
| 313 | integer :: i, j, k, is, ie, js, je, halos, nz | |
| 314 | ||
| 315 | 0 | do_debug = .false. ; if (present(debug)) do_debug = debug |
| 316 | 0 | halos = 0 ; if (present(halo)) halos = max(0,halo) |
| 317 | 0 | is = G%isc-halos ; ie = G%iec+halos ; js = G%jsc-halos ; je = G%jec+halos ; nz = GV%ke |
| 318 | ||
| 319 | 0 | if (allocated(tv%Spv_avg) .and. associated(tv%eqn_of_state)) then |
| 320 | 0 | if (associated(tv%p_surf)) then |
| 321 | 0 | do j=js,je ; do i=is,ie ; p_t(i,j) = tv%p_surf(i,j) ; enddo ; enddo |
| 322 | else | |
| 323 | 0 | do j=js,je ; do i=is,ie ; p_t(i,j) = 0.0 ; enddo ; enddo |
| 324 | endif | |
| 325 | 0 | do k=1,nz |
| 326 | 0 | do j=js,je ; do i=is,ie |
| 327 | 0 | dp(i,j) = GV%g_Earth*GV%H_to_RZ*h(i,j,k) |
| 328 | enddo ; enddo | |
| 329 | 0 | call avg_specific_vol(tv%T(:,:,k), tv%S(:,:,k), p_t, dp, G%HI, tv%eqn_of_state, tv%SpV_avg(:,:,k), halo) |
| 330 | 0 | if (k<nz) then ; do j=js,je ; do i=is,ie |
| 331 | 0 | p_t(i,j) = p_t(i,j) + dp(i,j) |
| 332 | enddo ; enddo ; endif | |
| 333 | enddo | |
| 334 | 0 | tv%valid_SpV_halo = halos |
| 335 | ||
| 336 | 0 | if (do_debug) then |
| 337 | 0 | call hchksum(h, "derived_thermo h", G%HI, haloshift=halos, unscale=GV%H_to_MKS) |
| 338 | 0 | if (associated(tv%p_surf)) call hchksum(tv%p_surf, "derived_thermo p_surf", G%HI, & |
| 339 | 0 | haloshift=halos, unscale=US%RL2_T2_to_Pa) |
| 340 | 0 | call hchksum(tv%T, "derived_thermo T", G%HI, haloshift=halos, unscale=US%C_to_degC) |
| 341 | 0 | call hchksum(tv%S, "derived_thermo S", G%HI, haloshift=halos, unscale=US%S_to_ppt) |
| 342 | endif | |
| 343 | 0 | elseif (allocated(tv%Spv_avg)) then |
| 344 | 0 | do k=1,nz ; SpV_lay(k) = 1.0 / GV%Rlay(k) ; enddo |
| 345 | 0 | do k=1,nz ; do j=js,je ; do i=is,ie |
| 346 | 0 | tv%SpV_avg(i,j,k) = SpV_lay(k) |
| 347 | enddo ; enddo ; enddo | |
| 348 | 0 | tv%valid_SpV_halo = halos |
| 349 | endif | |
| 350 | ||
| 351 | 0 | end subroutine calc_derived_thermo |
| 352 | ||
| 353 | ||
| 354 | !> Determine the column average specific volumes. | |
| 355 | 0 | subroutine find_col_avg_SpV(h, SpV_avg, tv, G, GV, US, halo_size) |
| 356 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure | |
| 357 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure | |
| 358 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 359 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & | |
| 360 | intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] | |
| 361 | real, dimension(SZI_(G),SZJ_(G)), & | |
| 362 | intent(inout) :: SpV_avg !< Column average specific volume [R-1 ~> m3 kg-1] | |
| 363 | ! SpV_avg is intent inout to retain excess halo values. | |
| 364 | type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any available | |
| 365 | !! thermodynamic fields. | |
| 366 | integer, optional, intent(in) :: halo_size !< width of halo points on which to work | |
| 367 | ||
| 368 | ! Local variables | |
| 369 | 0 | real :: h_tot(SZI_(G)) ! Sum of the layer thicknesses [H ~> m or kg m-2] |
| 370 | 0 | real :: SpV_x_h_tot(SZI_(G)) ! Vertical sum of the layer average specific volume times |
| 371 | ! the layer thicknesses [H R-1 ~> m4 kg-1 or m] | |
| 372 | real :: I_rho ! The inverse of the Boussiensq reference density [R-1 ~> m3 kg-1] | |
| 373 | 0 | real :: SpV_lay(SZK_(GV)) ! The inverse of the layer target potential densities [R-1 ~> m3 kg-1] |
| 374 | character(len=128) :: mesg ! A string for error messages | |
| 375 | integer :: i, j, k, is, ie, js, je, nz, halo | |
| 376 | ||
| 377 | 0 | halo = 0 ; if (present(halo_size)) halo = max(0,halo_size) |
| 378 | ||
| 379 | 0 | is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo |
| 380 | 0 | nz = GV%ke |
| 381 | ||
| 382 | 0 | if (GV%Boussinesq) then |
| 383 | 0 | I_rho = 1.0 / GV%Rho0 |
| 384 | 0 | do j=js,je ; do i=is,ie |
| 385 | 0 | SpV_avg(i,j) = I_rho |
| 386 | enddo ; enddo | |
| 387 | 0 | elseif (.not.allocated(tv%SpV_avg)) then |
| 388 | 0 | do k=1,nz ; Spv_lay(k) = 1.0 / GV%Rlay(k) ; enddo |
| 389 | 0 | do j=js,je |
| 390 | 0 | do i=is,ie ; SpV_x_h_tot(i) = 0.0 ; h_tot(i) = 0.0 ; enddo |
| 391 | 0 | do k=1,nz ; do i=is,ie |
| 392 | 0 | h_tot(i) = h_tot(i) + max(h(i,j,k), GV%H_subroundoff) |
| 393 | 0 | SpV_x_h_tot(i) = SpV_x_h_tot(i) + Spv_lay(k)*max(h(i,j,k), GV%H_subroundoff) |
| 394 | enddo ; enddo | |
| 395 | 0 | do i=is,ie ; SpV_avg(i,j) = SpV_x_h_tot(i) / h_tot(i) ; enddo |
| 396 | enddo | |
| 397 | else | |
| 398 | ! Check that SpV_avg has been set. | |
| 399 | 0 | if ((allocated(tv%SpV_avg)) .and. (tv%valid_SpV_halo < halo)) then |
| 400 | 0 | if (tv%valid_SpV_halo < 0) then |
| 401 | 0 | mesg = "invalid values of SpV_avg." |
| 402 | else | |
| 403 | write(mesg, '("insufficiently large SpV_avg halos of width ", i2, " but ", i2," is needed.")') & | |
| 404 | 0 | tv%valid_SpV_halo, halo |
| 405 | endif | |
| 406 | 0 | call MOM_error(FATAL, "find_col_avg_SpV called in fully non-Boussinesq mode with "//trim(mesg)) |
| 407 | endif | |
| 408 | ||
| 409 | 0 | do j=js,je |
| 410 | 0 | do i=is,ie ; SpV_x_h_tot(i) = 0.0 ; h_tot(i) = 0.0 ; enddo |
| 411 | 0 | do k=1,nz ; do i=is,ie |
| 412 | 0 | h_tot(i) = h_tot(i) + max(h(i,j,k), GV%H_subroundoff) |
| 413 | 0 | SpV_x_h_tot(i) = SpV_x_h_tot(i) + tv%SpV_avg(i,j,k)*max(h(i,j,k), GV%H_subroundoff) |
| 414 | enddo ; enddo | |
| 415 | 0 | do i=is,ie ; SpV_avg(i,j) = SpV_x_h_tot(i) / h_tot(i) ; enddo |
| 416 | enddo | |
| 417 | endif | |
| 418 | ||
| 419 | 0 | end subroutine find_col_avg_SpV |
| 420 | ||
| 421 | !> Calculate the integrated mass of the water column. | |
| 422 | 0 | subroutine find_col_mass(h, tv, G, GV, US, mass, p_bot, p_surf) |
| 423 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. | |
| 424 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. | |
| 425 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 426 | type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various | |
| 427 | !! thermodynamic variables. | |
| 428 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] | |
| 429 | real, dimension(SZI_(G),SZJ_(G)), intent(out) :: mass !< Integrated mass of the water column | |
| 430 | !! [R Z ~> kg m-2] | |
| 431 | real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: p_bot !< Bottom pressure = g * mass + psurf | |
| 432 | !! [R L2 T-2 ~> Pa] | |
| 433 | real, dimension(:,:), optional, pointer :: p_surf !< A pointer to surface pressure | |
| 434 | !! [R L2 T-2 ~> Pa] | |
| 435 | ||
| 436 | ! Local variables | |
| 437 | real :: I_gEarth ! The inverse of GV%g_Earth [T2 Z L-2 ~> s2 m-1] | |
| 438 | real, dimension(SZI_(G),SZJ_(G)) :: & | |
| 439 | 0 | z_top, & ! Height of the top of a layer [Z ~> m]. |
| 440 | 0 | z_bot, & ! Height of the bottom of a layer [Z ~> m]. |
| 441 | 0 | dp ! Change in hydrostatic pressure across a layer [R L2 T-2 ~> Pa]. |
| 442 | integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz | |
| 443 | ||
| 444 | 0 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec |
| 445 | 0 | isq = G%iscB ; ieq = G%iecB ; jsq = G%jscB ; jeq = G%jecB |
| 446 | 0 | nz = GV%ke |
| 447 | ||
| 448 | 0 | do j=js,je ; do i=is,ie ; mass(i,j) = 0.0 ; enddo ; enddo |
| 449 | 0 | if (GV%Boussinesq) then |
| 450 | 0 | if (associated(tv%eqn_of_state)) then |
| 451 | 0 | I_gEarth = 1.0 / GV%g_Earth |
| 452 | 0 | do j=jsq,jeq+1 ; do i=isq,ieq+1 ; z_bot(i,j) = 0.0 ; enddo ; enddo |
| 453 | 0 | do k=1,nz |
| 454 | ! NOTE: int_density_z expects z_top and z_bot values from [ij]sq to [ij]eq+1 | |
| 455 | 0 | do j=jsq,jeq+1 ; do i=isq,ieq+1 |
| 456 | 0 | z_top(i,j) = z_bot(i,j) |
| 457 | 0 | z_bot(i,j) = z_top(i,j) - GV%H_to_Z * h(i,j,k) |
| 458 | enddo ; enddo | |
| 459 | call int_density_dz(tv%T(:,:,k), tv%S(:,:,k), z_top, z_bot, 0.0, GV%Rho0, GV%g_Earth, & | |
| 460 | 0 | G%HI, tv%eqn_of_state, US, dp) |
| 461 | 0 | do j=js,je ; do i=is,ie |
| 462 | 0 | mass(i,j) = mass(i,j) + dp(i,j) * I_gEarth |
| 463 | enddo ; enddo | |
| 464 | enddo | |
| 465 | else | |
| 466 | 0 | do k=1,nz ; do j=js,je ; do i=is,ie |
| 467 | 0 | mass(i,j) = mass(i,j) + (GV%H_to_Z * GV%Rlay(k)) * h(i,j,k) |
| 468 | enddo ; enddo ; enddo | |
| 469 | endif | |
| 470 | else | |
| 471 | 0 | do k=1,nz ; do j=js,je ; do i=is,ie |
| 472 | 0 | mass(i,j) = mass(i,j) + GV%H_to_RZ * h(i,j,k) |
| 473 | enddo ; enddo ; enddo | |
| 474 | endif | |
| 475 | ||
| 476 | 0 | if (present(p_bot)) then |
| 477 | 0 | do j=js,je ; do i=is,ie |
| 478 | 0 | p_bot(i,j) = GV%g_Earth * mass(i,j) |
| 479 | enddo ; enddo | |
| 480 | 0 | if (present(p_surf) .and. associated(p_surf)) then ; do j=js,je ; do i=is,ie |
| 481 | 0 | p_bot(i,j) = p_bot(i,j) + p_surf(i,j) |
| 482 | enddo ; enddo ; endif | |
| 483 | endif | |
| 484 | ||
| 485 | 0 | end subroutine find_col_mass |
| 486 | ||
| 487 | !> Determine the in situ density averaged over a specified distance from the bottom, | |
| 488 | !! calculating it as the inverse of the mass-weighted average specific volume. | |
| 489 | 720 | subroutine find_rho_bottom(G, GV, US, tv, h, dz, pres_int, dz_avg, j, Rho_bot, h_bot, k_bot) |
| 490 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure | |
| 491 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure | |
| 492 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 493 | type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any available | |
| 494 | !! thermodynamic fields. | |
| 495 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & | |
| 496 | intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] | |
| 497 | real, dimension(SZI_(G),SZK_(GV)), & | |
| 498 | intent(in) :: dz !< Height change across layers [Z ~> m] | |
| 499 | real, dimension(SZI_(G),SZK_(GV)+1), & | |
| 500 | intent(in) :: pres_int !< Pressure at each interface [R L2 T-2 ~> Pa] | |
| 501 | real, dimension(SZI_(G)), intent(in) :: dz_avg !< The vertical distance over which to average [Z ~> m] | |
| 502 | integer, intent(in) :: j !< j-index of row to work on | |
| 503 | real, dimension(SZI_(G)), intent(out) :: Rho_bot !< Near-bottom density [R ~> kg m-3]. | |
| 504 | real, dimension(SZI_(G)), intent(out) :: h_bot !< Bottom boundary layer thickness [H ~> m or kg m-2] | |
| 505 | integer, dimension(SZI_(G)), intent(out) :: k_bot !< Bottom boundary layer top layer index | |
| 506 | ||
| 507 | ! Local variables | |
| 508 | 1440 | real :: hb(SZI_(G)) ! Running sum of the thickness in the bottom boundary layer [H ~> m or kg m-2] |
| 509 | 1440 | real :: SpV_h_bot(SZI_(G)) ! Running sum of the specific volume times thickness in the bottom |
| 510 | ! boundary layer [H R-1 ~> m4 kg-1 or m] | |
| 511 | 1440 | real :: dz_bbl_rem(SZI_(G)) ! Vertical extent of the boundary layer that has yet to be accounted |
| 512 | ! for [Z ~> m] | |
| 513 | 1440 | real :: h_bbl_frac(SZI_(G)) ! Thickness of the fractional layer that makes up the top of the |
| 514 | ! boundary layer [H ~> m or kg m-2] | |
| 515 | 1440 | real :: T_bbl(SZI_(G)) ! Temperature of the fractional layer that makes up the top of the |
| 516 | ! boundary layer [C ~> degC] | |
| 517 | 1440 | real :: S_bbl(SZI_(G)) ! Salinity of the fractional layer that makes up the top of the |
| 518 | ! boundary layer [S ~> ppt] | |
| 519 | 1440 | real :: P_bbl(SZI_(G)) ! Pressure the top of the boundary layer [R L2 T-2 ~> Pa] |
| 520 | 1440 | real :: dp(SZI_(G)) ! Pressure change across the fractional layer that makes up the top |
| 521 | ! of the boundary layer [R L2 T-2 ~> Pa] | |
| 522 | 1440 | real :: SpV_bbl(SZI_(G)) ! In situ specific volume of the fractional layer that makes up the |
| 523 | ! top of the boundary layer [R-1 ~> m3 kg-1] | |
| 524 | real :: frac_in ! The fraction of a layer that is within the bottom boundary layer [nondim] | |
| 525 | 1440 | logical :: do_i(SZI_(G)), do_any |
| 526 | logical :: use_EOS | |
| 527 | integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state | |
| 528 | integer :: i, k, is, ie, nz | |
| 529 | ||
| 530 | 720 | is = G%isc ; ie = G%iec ; nz = GV%ke |
| 531 | ||
| 532 | 720 | use_EOS = associated(tv%T) .and. associated(tv%S) .and. associated(tv%eqn_of_state) |
| 533 | ||
| 534 | 720 | if (GV%Boussinesq .or. GV%semi_Boussinesq .or. .not.allocated(tv%SpV_avg)) then |
| 535 | 87120 | do i=is,ie |
| 536 | 87120 | rho_bot(i) = GV%Rho0 |
| 537 | enddo | |
| 538 | ||
| 539 | ! Obtain bottom boundary layer thickness and index of top layer | |
| 540 | 87120 | do i=is,ie |
| 541 | 86400 | hb(i) = 0.0 ; h_bot(i) = 0.0 ; k_bot(i) = nz |
| 542 | 86400 | dz_bbl_rem(i) = G%mask2dT(i,j) * max(0.0, dz_avg(i)) |
| 543 | 86400 | do_i(i) = .true. |
| 544 | 87120 | if (G%mask2dT(i,j) <= 0.0) then |
| 545 | 26232 | h_bbl_frac(i) = 0.0 |
| 546 | 26232 | do_i(i) = .false. |
| 547 | endif | |
| 548 | enddo | |
| 549 | ||
| 550 | 720 | do k=nz,1,-1 |
| 551 | 720 | do_any = .false. |
| 552 | 87120 | do i=is,ie ; if (do_i(i)) then |
| 553 | 60168 | if (dz(i,k) < dz_bbl_rem(i)) then |
| 554 | ! This layer is fully within the averaging depth. | |
| 555 | 0 | dz_bbl_rem(i) = dz_bbl_rem(i) - dz(i,k) |
| 556 | 0 | hb(i) = hb(i) + h(i,j,k) |
| 557 | 0 | k_bot(i) = k |
| 558 | 0 | do_any = .true. |
| 559 | else | |
| 560 | 60168 | if (dz(i,k) > 0.0) then |
| 561 | 60168 | frac_in = dz_bbl_rem(i) / dz(i,k) |
| 562 | 60168 | if (frac_in >= 0.5) k_bot(i) = k ! update bbl top index if >= 50% of layer |
| 563 | else | |
| 564 | 0 | frac_in = 0.0 |
| 565 | endif | |
| 566 | 60168 | h_bbl_frac(i) = frac_in * h(i,j,k) |
| 567 | 60168 | dz_bbl_rem(i) = 0.0 |
| 568 | 60168 | do_i(i) = .false. |
| 569 | endif | |
| 570 | endif ; enddo | |
| 571 | 720 | if (.not.do_any) exit |
| 572 | enddo | |
| 573 | 87120 | do i=is,ie ; if (do_i(i)) then |
| 574 | ! The nominal bottom boundary layer is thicker than the water column, but layer 1 is | |
| 575 | ! already included in the averages. These values are set so that the call to find | |
| 576 | ! the layer-average specific volume will behave sensibly. | |
| 577 | 0 | h_bbl_frac(i) = 0.0 |
| 578 | endif ; enddo | |
| 579 | ||
| 580 | 87120 | do i=is,ie |
| 581 | 86400 | if (hb(i) + h_bbl_frac(i) < GV%H_subroundoff) h_bbl_frac(i) = GV%H_subroundoff |
| 582 | 87120 | h_bot(i) = hb(i) + h_bbl_frac(i) |
| 583 | enddo | |
| 584 | ||
| 585 | else | |
| 586 | ! Check that SpV_avg has been set. | |
| 587 | 0 | if (tv%valid_SpV_halo < 0) call MOM_error(FATAL, & |
| 588 | 0 | "find_rho_bottom called in fully non-Boussinesq mode with invalid values of SpV_avg.") |
| 589 | ||
| 590 | ! Set the bottom density to the inverse of the in situ specific volume averaged over the | |
| 591 | ! specified distance, with care taken to avoid having compressibility lead to an imprint | |
| 592 | ! of the layer thicknesses on this density. | |
| 593 | 0 | do i=is,ie |
| 594 | 0 | hb(i) = 0.0 ; SpV_h_bot(i) = 0.0 ; h_bot(i) = 0.0 ; k_bot(i) = nz |
| 595 | 0 | dz_bbl_rem(i) = G%mask2dT(i,j) * max(0.0, dz_avg(i)) |
| 596 | 0 | do_i(i) = .true. |
| 597 | 0 | if (G%mask2dT(i,j) <= 0.0) then |
| 598 | ! Set acceptable values for calling the equation of state over land. | |
| 599 | 0 | T_bbl(i) = 0.0 ; S_bbl(i) = 0.0 ; dp(i) = 0.0 ; P_bbl(i) = 0.0 |
| 600 | 0 | SpV_bbl(i) = 1.0 ! This value is arbitrary, provided it is non-zero. |
| 601 | 0 | h_bbl_frac(i) = 0.0 |
| 602 | 0 | do_i(i) = .false. |
| 603 | endif | |
| 604 | enddo | |
| 605 | ||
| 606 | 0 | do k=nz,1,-1 |
| 607 | 0 | do_any = .false. |
| 608 | 0 | do i=is,ie ; if (do_i(i)) then |
| 609 | 0 | if (dz(i,k) < dz_bbl_rem(i)) then |
| 610 | ! This layer is fully within the averaging depth. | |
| 611 | 0 | SpV_h_bot(i) = SpV_h_bot(i) + h(i,j,k) * tv%SpV_avg(i,j,k) |
| 612 | 0 | dz_bbl_rem(i) = dz_bbl_rem(i) - dz(i,k) |
| 613 | 0 | hb(i) = hb(i) + h(i,j,k) |
| 614 | 0 | k_bot(i) = k |
| 615 | 0 | do_any = .true. |
| 616 | else | |
| 617 | 0 | if (dz(i,k) > 0.0) then |
| 618 | 0 | frac_in = dz_bbl_rem(i) / dz(i,k) |
| 619 | 0 | if (frac_in >= 0.5) k_bot(i) = k ! update bbl top index if >= 50% of layer |
| 620 | else | |
| 621 | 0 | frac_in = 0.0 |
| 622 | endif | |
| 623 | 0 | if (use_EOS) then |
| 624 | ! Store the properties of this layer to determine the average | |
| 625 | ! specific volume of the portion that is within the BBL. | |
| 626 | 0 | T_bbl(i) = tv%T(i,j,k) ; S_bbl(i) = tv%S(i,j,k) |
| 627 | 0 | dp(i) = frac_in * (GV%g_Earth*GV%H_to_RZ * h(i,j,k)) |
| 628 | 0 | P_bbl(i) = pres_int(i,K) + (1.0-frac_in) * (GV%g_Earth*GV%H_to_RZ * h(i,j,k)) |
| 629 | else | |
| 630 | 0 | SpV_bbl(i) = tv%SpV_avg(i,j,k) |
| 631 | endif | |
| 632 | 0 | h_bbl_frac(i) = frac_in * h(i,j,k) |
| 633 | 0 | dz_bbl_rem(i) = 0.0 |
| 634 | 0 | do_i(i) = .false. |
| 635 | endif | |
| 636 | endif ; enddo | |
| 637 | 0 | if (.not.do_any) exit |
| 638 | enddo | |
| 639 | 0 | do i=is,ie ; if (do_i(i)) then |
| 640 | ! The nominal bottom boundary layer is thicker than the water column, but layer 1 is | |
| 641 | ! already included in the averages. These values are set so that the call to find | |
| 642 | ! the layer-average specific volume will behave sensibly. | |
| 643 | 0 | if (use_EOS) then |
| 644 | 0 | T_bbl(i) = tv%T(i,j,1) ; S_bbl(i) = tv%S(i,j,1) |
| 645 | 0 | dp(i) = 0.0 |
| 646 | 0 | P_bbl(i) = pres_int(i,1) |
| 647 | else | |
| 648 | 0 | SpV_bbl(i) = tv%SpV_avg(i,j,1) |
| 649 | endif | |
| 650 | 0 | h_bbl_frac(i) = 0.0 |
| 651 | endif ; enddo | |
| 652 | ||
| 653 | 0 | if (use_EOS) then |
| 654 | ! Find the average specific volume of the fractional layer atop the BBL. | |
| 655 | 0 | EOSdom(:) = EOS_domain(G%HI) |
| 656 | 0 | call average_specific_vol(T_bbl, S_bbl, P_bbl, dp, SpV_bbl, tv%eqn_of_state, EOSdom) |
| 657 | endif | |
| 658 | ||
| 659 | 0 | do i=is,ie |
| 660 | 0 | if (hb(i) + h_bbl_frac(i) < GV%H_subroundoff) h_bbl_frac(i) = GV%H_subroundoff |
| 661 | 0 | rho_bot(i) = G%mask2dT(i,j) * (hb(i) + h_bbl_frac(i)) / (SpV_h_bot(i) + h_bbl_frac(i)*SpV_bbl(i)) |
| 662 | 0 | h_bot(i) = hb(i) + h_bbl_frac(i) |
| 663 | enddo | |
| 664 | endif | |
| 665 | ||
| 666 | 720 | end subroutine find_rho_bottom |
| 667 | ||
| 668 | ||
| 669 | !> Converts thickness from geometric height units to thickness units, perhaps via an | |
| 670 | !! inversion of the integral of the density in pressure using variables stored in | |
| 671 | !! the thermo_var_ptrs type when in non-Boussinesq mode. | |
| 672 | 1 | subroutine dz_to_thickness_tv(dz, tv, h, G, GV, US, halo_size) |
| 673 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure | |
| 674 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure | |
| 675 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 676 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & | |
| 677 | intent(in) :: dz !< Geometric layer thicknesses in height units [Z ~> m] | |
| 678 | type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various | |
| 679 | !! thermodynamic variables | |
| 680 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & | |
| 681 | intent(inout) :: h !< Output thicknesses in thickness units [H ~> m or kg m-2]. | |
| 682 | !! This is essentially intent out, but declared as intent | |
| 683 | !! inout to preserve any initialized values in halo points. | |
| 684 | integer, optional, intent(in) :: halo_size !< Width of halo within which to | |
| 685 | !! calculate thicknesses | |
| 686 | ! Local variables | |
| 687 | integer :: i, j, k, is, ie, js, je, halo, nz | |
| 688 | ||
| 689 | 1 | halo = 0 ; if (present(halo_size)) halo = max(0,halo_size) |
| 690 | 1 | is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo ; nz = GV%ke |
| 691 | ||
| 692 | 1 | if (GV%Boussinesq) then |
| 693 | 544576 | do k=1,nz ; do j=js,je ; do i=is,ie |
| 694 | 544500 | h(i,j,k) = GV%Z_to_H * dz(i,j,k) |
| 695 | enddo ; enddo ; enddo | |
| 696 | else | |
| 697 | 0 | if (associated(tv%eqn_of_state)) then |
| 698 | 0 | if (associated(tv%p_surf)) then |
| 699 | 0 | call dz_to_thickness_EOS(dz, tv%T, tv%S, tv%eqn_of_state, h, G, GV, US, halo, tv%p_surf) |
| 700 | else | |
| 701 | 0 | call dz_to_thickness_EOS(dz, tv%T, tv%S, tv%eqn_of_state, h, G, GV, US, halo) |
| 702 | endif | |
| 703 | else | |
| 704 | 0 | do k=1,nz ; do j=js,je ; do i=is,ie |
| 705 | 0 | h(i,j,k) = (GV%RZ_to_H * GV%Rlay(k)) * dz(i,j,k) |
| 706 | enddo ; enddo ; enddo | |
| 707 | endif | |
| 708 | endif | |
| 709 | ||
| 710 | 1 | end subroutine dz_to_thickness_tv |
| 711 | ||
| 712 | !> Converts thickness from geometric height units to thickness units, working via an | |
| 713 | !! inversion of the integral of the density in pressure when in non-Boussinesq mode. | |
| 714 | 0 | subroutine dz_to_thickness_EOS(dz, Temp, Saln, EoS, h, G, GV, US, halo_size, p_surf) |
| 715 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure | |
| 716 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure | |
| 717 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 718 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & | |
| 719 | intent(in) :: dz !< Geometric layer thicknesses in height units [Z ~> m] | |
| 720 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & | |
| 721 | intent(in) :: Temp !< Input layer temperatures [C ~> degC] | |
| 722 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & | |
| 723 | intent(in) :: Saln !< Input layer salinities [S ~> ppt] | |
| 724 | type(EOS_type), intent(in) :: EoS !< Equation of state structure | |
| 725 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & | |
| 726 | intent(inout) :: h !< Output thicknesses in thickness units [H ~> m or kg m-2]. | |
| 727 | !! This is essentially intent out, but declared as intent | |
| 728 | !! inout to preserve any initialized values in halo points. | |
| 729 | integer, optional, intent(in) :: halo_size !< Width of halo within which to | |
| 730 | !! calculate thicknesses | |
| 731 | real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: p_surf !< Surface pressures [R L2 T-2 ~> Pa] | |
| 732 | ! Local variables | |
| 733 | real, dimension(SZI_(G),SZJ_(G)) :: & | |
| 734 | 0 | p_top, p_bot ! Pressure at the interfaces above and below a layer [R L2 T-2 ~> Pa] |
| 735 | 0 | real :: dp(SZI_(G),SZJ_(G)) ! Pressure change across a layer [R L2 T-2 ~> Pa] |
| 736 | 0 | real :: dz_geo(SZI_(G),SZJ_(G)) ! The change in geopotential height across a layer [L2 T-2 ~> m2 s-2] |
| 737 | 0 | real :: rho(SZI_(G)) ! The in situ density [R ~> kg m-3] |
| 738 | real :: dp_adj ! The amount by which to change the bottom pressure in an | |
| 739 | ! iteration [R L2 T-2 ~> Pa] | |
| 740 | real :: I_gEarth ! Unit conversion factors divided by the gravitational | |
| 741 | ! acceleration [H T2 R-1 L-2 ~> s2 m2 kg-1 or s2 m-1] | |
| 742 | 0 | logical :: do_more(SZI_(G),SZJ_(G)) ! If true, additional iterations would be beneficial. |
| 743 | logical :: do_any ! True if there are points in this layer that need more itertions. | |
| 744 | integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state | |
| 745 | integer :: i, j, k, is, ie, js, je, halo, nz | |
| 746 | integer :: itt, max_itt | |
| 747 | ||
| 748 | 0 | halo = 0 ; if (present(halo_size)) halo = max(0,halo_size) |
| 749 | 0 | is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo ; nz = GV%ke |
| 750 | 0 | max_itt = 10 |
| 751 | ||
| 752 | 0 | if (GV%Boussinesq) then |
| 753 | 0 | do k=1,nz ; do j=js,je ; do i=is,ie |
| 754 | 0 | h(i,j,k) = GV%Z_to_H * dz(i,j,k) |
| 755 | enddo ; enddo ; enddo | |
| 756 | else | |
| 757 | 0 | I_gEarth = GV%RZ_to_H / GV%g_Earth |
| 758 | ||
| 759 | 0 | if (present(p_surf)) then |
| 760 | 0 | do j=js,je ; do i=is,ie |
| 761 | 0 | p_bot(i,j) = 0.0 ; p_top(i,j) = p_surf(i,j) |
| 762 | enddo ; enddo | |
| 763 | else | |
| 764 | 0 | do j=js,je ; do i=is,ie |
| 765 | 0 | p_bot(i,j) = 0.0 ; p_top(i,j) = 0.0 |
| 766 | enddo ; enddo | |
| 767 | endif | |
| 768 | 0 | EOSdom(:) = EOS_domain(G%HI) |
| 769 | ||
| 770 | ! The iterative approach here is inherited from very old code that was in the | |
| 771 | ! MOM_state_initialization module. It does converge, but it is very inefficient and | |
| 772 | ! should be revised, although doing so would change answers in non-Boussinesq mode. | |
| 773 | 0 | do k=1,nz |
| 774 | 0 | do j=js,je |
| 775 | 0 | do i=is,ie ; p_top(i,j) = p_bot(i,j) ; enddo |
| 776 | call calculate_density(Temp(:,j,k), Saln(:,j,k), p_top(:,j), rho, & | |
| 777 | 0 | EoS, EOSdom) |
| 778 | ! The following two expressions are mathematically equivalent. | |
| 779 | 0 | if (GV%semi_Boussinesq) then |
| 780 | 0 | do i=is,ie |
| 781 | 0 | p_bot(i,j) = p_top(i,j) + (GV%g_Earth*GV%H_to_Z) * ((GV%Z_to_H*dz(i,j,k)) * rho(i)) |
| 782 | 0 | dp(i,j) = (GV%g_Earth*GV%H_to_Z) * ((GV%Z_to_H*dz(i,j,k)) * rho(i)) |
| 783 | enddo | |
| 784 | else | |
| 785 | 0 | do i=is,ie |
| 786 | 0 | p_bot(i,j) = p_top(i,j) + rho(i) * (GV%g_Earth * dz(i,j,k)) |
| 787 | 0 | dp(i,j) = rho(i) * (GV%g_Earth * dz(i,j,k)) |
| 788 | enddo | |
| 789 | endif | |
| 790 | enddo | |
| 791 | ||
| 792 | 0 | do_more(:,:) = .true. |
| 793 | 0 | do itt=1,max_itt |
| 794 | 0 | do_any = .false. |
| 795 | 0 | call int_specific_vol_dp(Temp(:,:,k), Saln(:,:,k), p_top, p_bot, 0.0, G%HI, EoS, US, dz_geo) |
| 796 | 0 | if (itt < max_itt) then ; do j=js,je |
| 797 | 0 | call calculate_density(Temp(:,j,k), Saln(:,j,k), p_bot(:,j), rho, EoS, EOSdom) |
| 798 | ! Use Newton's method to correct the bottom value. | |
| 799 | ! The hydrostatic equation is sufficiently linear that no bounds-checking is needed. | |
| 800 | 0 | if (GV%semi_Boussinesq) then |
| 801 | 0 | do i=is,ie |
| 802 | 0 | dp_adj = rho(i) * ((GV%g_Earth*GV%H_to_Z)*(GV%Z_to_H*dz(i,j,k)) - dz_geo(i,j)) |
| 803 | 0 | p_bot(i,j) = p_bot(i,j) + dp_adj |
| 804 | 0 | dp(i,j) = dp(i,j) + dp_adj |
| 805 | enddo | |
| 806 | 0 | do_any = .true. ! To avoid changing answers, always use the maximum number of itertions. |
| 807 | else | |
| 808 | 0 | do i=is,ie ; if (do_more(i,j)) then |
| 809 | 0 | dp_adj = rho(i) * (GV%g_Earth*dz(i,j,k) - dz_geo(i,j)) |
| 810 | 0 | p_bot(i,j) = p_bot(i,j) + dp_adj |
| 811 | 0 | dp(i,j) = dp(i,j) + dp_adj |
| 812 | ! Check for convergence to roundoff. | |
| 813 | 0 | do_more(i,j) = (abs(dp_adj) > 1.0e-15*dp(i,j)) |
| 814 | 0 | if (do_more(i,j)) do_any = .true. |
| 815 | endif ; enddo | |
| 816 | endif | |
| 817 | enddo ; endif | |
| 818 | 0 | if (.not.do_any) exit |
| 819 | enddo | |
| 820 | ||
| 821 | 0 | if (GV%semi_Boussinesq) then |
| 822 | 0 | do j=js,je ; do i=is,ie |
| 823 | 0 | h(i,j,k) = (p_bot(i,j) - p_top(i,j)) * I_gEarth |
| 824 | enddo ; enddo | |
| 825 | else | |
| 826 | 0 | do j=js,je ; do i=is,ie |
| 827 | 0 | h(i,j,k) = dp(i,j) * I_gEarth |
| 828 | enddo ; enddo | |
| 829 | endif | |
| 830 | enddo | |
| 831 | endif | |
| 832 | ||
| 833 | 0 | end subroutine dz_to_thickness_EOS |
| 834 | ||
| 835 | !> Converts thickness from geometric height units to thickness units, perhaps using | |
| 836 | !! a simple conversion factor that may be problematic in non-Boussinesq mode. | |
| 837 | 0 | subroutine dz_to_thickness_simple(dz, h, G, GV, US, halo_size, layer_mode) |
| 838 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure | |
| 839 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure | |
| 840 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 841 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & | |
| 842 | intent(in) :: dz !< Geometric layer thicknesses in height units [Z ~> m] | |
| 843 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & | |
| 844 | intent(inout) :: h !< Output thicknesses in thickness units [H ~> m or kg m-2]. | |
| 845 | !! This is essentially intent out, but declared as intent | |
| 846 | !! inout to preserve any initialized values in halo points. | |
| 847 | integer, optional, intent(in) :: halo_size !< Width of halo within which to | |
| 848 | !! calculate thicknesses | |
| 849 | logical, optional, intent(in) :: layer_mode !< If present and true, do the conversion that | |
| 850 | !! is appropriate in pure isopycnal layer mode with | |
| 851 | !! no state variables or equation of state. Otherwise | |
| 852 | !! use a simple constant rescaling factor and avoid the | |
| 853 | !! use of GV%Rlay. | |
| 854 | ! Local variables | |
| 855 | logical :: layered ! If true and the model is non-Boussinesq, do calculations appropriate for use | |
| 856 | ! in pure isopycnal layered mode with no state variables or equation of state. | |
| 857 | integer :: i, j, k, is, ie, js, je, halo, nz | |
| 858 | ||
| 859 | 0 | halo = 0 ; if (present(halo_size)) halo = max(0,halo_size) |
| 860 | 0 | layered = .false. ; if (present(layer_mode)) layered = layer_mode |
| 861 | 0 | is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo ; nz = GV%ke |
| 862 | ||
| 863 | 0 | if (GV%Boussinesq) then |
| 864 | 0 | do k=1,nz ; do j=js,je ; do i=is,ie |
| 865 | 0 | h(i,j,k) = GV%Z_to_H * dz(i,j,k) |
| 866 | enddo ; enddo ; enddo | |
| 867 | 0 | elseif (layered) then |
| 868 | 0 | do k=1,nz ; do j=js,je ; do i=is,ie |
| 869 | 0 | h(i,j,k) = (GV%RZ_to_H * GV%Rlay(k)) * dz(i,j,k) |
| 870 | enddo ; enddo ; enddo | |
| 871 | else | |
| 872 | 0 | do k=1,nz ; do j=js,je ; do i=is,ie |
| 873 | 0 | h(i,j,k) = (US%Z_to_m * GV%m_to_H) * dz(i,j,k) |
| 874 | enddo ; enddo ; enddo | |
| 875 | endif | |
| 876 | ||
| 877 | 0 | end subroutine dz_to_thickness_simple |
| 878 | ||
| 879 | !> Converts layer thicknesses in thickness units to the vertical distance between edges in height | |
| 880 | !! units, perhaps by multiplication by the precomputed layer-mean specific volume stored in an | |
| 881 | !! array in the thermo_var_ptrs type when in non-Boussinesq mode. | |
| 882 | 144 | subroutine thickness_to_dz_3d(h, tv, dz, G, GV, US, halo_size, do_offload) |
| 883 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure | |
| 884 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure | |
| 885 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 886 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & | |
| 887 | intent(in) :: h !< Input thicknesses in thickness units [H ~> m or kg m-2]. | |
| 888 | type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various | |
| 889 | !! thermodynamic variables | |
| 890 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & | |
| 891 | intent(inout) :: dz !< Geometric layer thicknesses in height units [Z ~> m] | |
| 892 | !! This is essentially intent out, but declared as intent | |
| 893 | !! inout to preserve any initialized values in halo points. | |
| 894 | integer, optional, intent(in) :: halo_size !< Width of halo within which to | |
| 895 | !! calculate thicknesses | |
| 896 | logical, optional, intent(in) :: do_offload !< If .true., only uses data calculates dz | |
| 897 | !! on GPU (default .false.) | |
| 898 | ! Local variables | |
| 899 | character(len=128) :: mesg ! A string for error messages | |
| 900 | integer :: i, j, k, is, ie, js, je, halo, nz | |
| 901 | logical :: use_doconcurrent | |
| 902 | ||
| 903 | ! guard to allow turning off/on do concurrent | |
| 904 | 144 | use_doconcurrent = .false. |
| 905 | 144 | if (present(do_offload)) use_doconcurrent = do_offload |
| 906 | ||
| 907 | 144 | halo = 0 ; if (present(halo_size)) halo = max(0,halo_size) |
| 908 | 144 | is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo ; nz = GV%ke |
| 909 | ||
| 910 | 144 | if ((.not.GV%Boussinesq) .and. allocated(tv%SpV_avg)) then |
| 911 | 0 | if ((allocated(tv%SpV_avg)) .and. (tv%valid_SpV_halo < halo)) then |
| 912 | 0 | if (tv%valid_SpV_halo < 0) then |
| 913 | 0 | mesg = "invalid values of SpV_avg." |
| 914 | else | |
| 915 | write(mesg, '("insufficiently large SpV_avg halos of width ", i2, " but ", i2," is needed.")') & | |
| 916 | 0 | tv%valid_SpV_halo, halo |
| 917 | endif | |
| 918 | 0 | call MOM_error(FATAL, "thickness_to_dz called in fully non-Boussinesq mode with "//trim(mesg)) |
| 919 | endif | |
| 920 | 0 | if (use_doconcurrent) then |
| 921 | 0 | do concurrent (k=1:nz, j=js:je, i=is:ie) |
| 922 | 0 | dz(i,j,k) = GV%H_to_RZ * h(i,j,k) * tv%SpV_avg(i,j,k) |
| 923 | enddo | |
| 924 | else | |
| 925 | 0 | do k=1,nz ; do j=js,je ; do i=is,ie |
| 926 | 0 | dz(i,j,k) = GV%H_to_RZ * h(i,j,k) * tv%SpV_avg(i,j,k) |
| 927 | enddo ; enddo ; enddo | |
| 928 | endif | |
| 929 | else | |
| 930 | 144 | if (use_doconcurrent) then |
| 931 | 84 | do concurrent (k=1:nz, j=js:je, i=is:ie) |
| 932 | 48298908 | dz(i,j,k) = GV%H_to_Z * h(i,j,k) |
| 933 | enddo | |
| 934 | else | |
| 935 | 33662760 | do k=1,nz ; do j=js,je ; do i=is,ie |
| 936 | 33658200 | dz(i,j,k) = GV%H_to_Z * h(i,j,k) |
| 937 | enddo ; enddo ; enddo | |
| 938 | endif | |
| 939 | endif | |
| 940 | ||
| 941 | 144 | end subroutine thickness_to_dz_3d |
| 942 | ||
| 943 | ||
| 944 | !> Converts a vertical i- / k- slice of layer thicknesses in thickness units to the vertical | |
| 945 | !! distance between edges in height units, perhaps by multiplication by the precomputed layer-mean | |
| 946 | !! specific volume stored in an array in the thermo_var_ptrs type when in non-Boussinesq mode. | |
| 947 | 2904 | subroutine thickness_to_dz_jslice(h, tv, dz, j, G, GV, halo_size) |
| 948 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure | |
| 949 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure | |
| 950 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & | |
| 951 | intent(in) :: h !< Input thicknesses in thickness units [H ~> m or kg m-2]. | |
| 952 | type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various | |
| 953 | !! thermodynamic variables | |
| 954 | real, dimension(SZI_(G),SZK_(GV)), & | |
| 955 | intent(inout) :: dz !< Geometric layer thicknesses in height units [Z ~> m] | |
| 956 | !! This is essentially intent out, but declared as intent | |
| 957 | !! inout to preserve any initialized values in halo points. | |
| 958 | integer, intent(in) :: j !< The second (j-) index of the input thicknesses to work with | |
| 959 | integer, optional, intent(in) :: halo_size !< Width of halo within which to | |
| 960 | !! calculate thicknesses | |
| 961 | ! Local variables | |
| 962 | character(len=128) :: mesg ! A string for error messages | |
| 963 | integer :: i, k, is, ie, halo, nz | |
| 964 | ||
| 965 | 2904 | halo = 0 ; if (present(halo_size)) halo = max(0,halo_size) |
| 966 | 2904 | is = G%isc-halo ; ie = G%iec+halo ; nz = GV%ke |
| 967 | ||
| 968 | 2904 | if ((.not.GV%Boussinesq) .and. allocated(tv%SpV_avg)) then |
| 969 | 0 | if ((allocated(tv%SpV_avg)) .and. (tv%valid_SpV_halo < halo)) then |
| 970 | 0 | if (tv%valid_SpV_halo < 0) then |
| 971 | 0 | mesg = "invalid values of SpV_avg." |
| 972 | else | |
| 973 | write(mesg, '("insufficiently large SpV_avg halos of width ", i2, " but ", i2," is needed.")') & | |
| 974 | 0 | tv%valid_SpV_halo, halo |
| 975 | endif | |
| 976 | 0 | call MOM_error(FATAL, "thickness_to_dz called in fully non-Boussinesq mode with "//trim(mesg)) |
| 977 | endif | |
| 978 | ||
| 979 | 0 | do k=1,nz ; do i=is,ie |
| 980 | 0 | dz(i,k) = GV%H_to_RZ * h(i,j,k) * tv%SpV_avg(i,j,k) |
| 981 | enddo ; enddo | |
| 982 | else | |
| 983 | 26468304 | do k=1,nz ; do i=is,ie |
| 984 | 26465400 | dz(i,k) = GV%H_to_Z * h(i,j,k) |
| 985 | enddo ; enddo | |
| 986 | endif | |
| 987 | ||
| 988 | 2904 | end subroutine thickness_to_dz_jslice |
| 989 | ||
| 990 | ||
| 991 | !> Convert mixed layer depths in height units into the thickness of water in the mixed | |
| 992 | !! in thickness units. | |
| 993 | 12 | subroutine convert_MLD_to_ML_thickness(MLD_in, h, h_MLD, tv, G, GV, halo) |
| 994 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure | |
| 995 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure | |
| 996 | real, dimension(SZI_(G),SZJ_(G)), & | |
| 997 | intent(in) :: MLD_in !< Input mixed layer depth [Z ~> m]. | |
| 998 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & | |
| 999 | intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] | |
| 1000 | real, dimension(SZI_(G),SZJ_(G)), & | |
| 1001 | intent(out) :: h_MLD !< Thickness of water in the mixed layer [H ~> m or kg m-2] | |
| 1002 | type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any available | |
| 1003 | !! thermodynamic fields. | |
| 1004 | integer, optional, intent(in) :: halo !< Halo width over which to calculate frazil | |
| 1005 | ||
| 1006 | ! Local variables | |
| 1007 | 24 | real :: MLD_rem(SZI_(G)) ! The vertical extent of the MLD_in that has not yet been accounted for [Z ~> m] |
| 1008 | character(len=128) :: mesg ! A string for error messages | |
| 1009 | logical :: keep_going | |
| 1010 | integer :: i, j, k, is, ie, js, je, nz, halos | |
| 1011 | ||
| 1012 | 12 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke |
| 1013 | ||
| 1014 | 12 | halos = 0 ; if (present(halo)) halos = halo |
| 1015 | 12 | if (present(halo)) then |
| 1016 | 0 | is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo |
| 1017 | endif | |
| 1018 | ||
| 1019 | 12 | if (GV%Boussinesq .or. (.not.allocated(tv%SpV_avg))) then |
| 1020 | 87132 | do j=js,je ; do i=is,ie |
| 1021 | 87120 | h_MLD(i,j) = GV%Z_to_H * MLD_in(i,j) |
| 1022 | enddo ; enddo | |
| 1023 | else ! The fully non-Boussinesq conversion between height in MLD_in and thickness. | |
| 1024 | 0 | if ((allocated(tv%SpV_avg)) .and. (tv%valid_SpV_halo < halos)) then |
| 1025 | 0 | if (tv%valid_SpV_halo < 0) then |
| 1026 | 0 | mesg = "invalid values of SpV_avg." |
| 1027 | else | |
| 1028 | write(mesg, '("insufficiently large SpV_avg halos of width ", i2, " but ", i2," is needed.")') & | |
| 1029 | 0 | tv%valid_SpV_halo, halos |
| 1030 | endif | |
| 1031 | 0 | call MOM_error(FATAL, "convert_MLD_to_ML_thickness called in fully non-Boussinesq mode with "//trim(mesg)) |
| 1032 | endif | |
| 1033 | ||
| 1034 | 0 | do j=js,je |
| 1035 | 0 | do i=is,ie ; MLD_rem(i) = MLD_in(i,j) ; h_MLD(i,j) = 0.0 ; enddo |
| 1036 | 0 | do k=1,nz |
| 1037 | 0 | keep_going = .false. |
| 1038 | 0 | do i=is,ie ; if (MLD_rem(i) > 0.0) then |
| 1039 | 0 | if (MLD_rem(i) > GV%H_to_RZ * h(i,j,k) * tv%SpV_avg(i,j,k)) then |
| 1040 | 0 | h_MLD(i,j) = h_MLD(i,j) + h(i,j,k) |
| 1041 | 0 | MLD_rem(i) = MLD_rem(i) - GV%H_to_RZ * h(i,j,k) * tv%SpV_avg(i,j,k) |
| 1042 | 0 | keep_going = .true. |
| 1043 | else | |
| 1044 | 0 | h_MLD(i,j) = h_MLD(i,j) + GV%RZ_to_H * MLD_rem(i) / tv%SpV_avg(i,j,k) |
| 1045 | 0 | MLD_rem(i) = 0.0 |
| 1046 | endif | |
| 1047 | endif ; enddo | |
| 1048 | 0 | if (.not.keep_going) exit |
| 1049 | enddo | |
| 1050 | enddo | |
| 1051 | endif | |
| 1052 | ||
| 1053 | 12 | end subroutine convert_MLD_to_ML_thickness |
| 1054 | ||
| 1055 | end module MOM_interface_heights |