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