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 | !> This module contains routines that implement physical fluxes of tracers (e.g. due | |
| 6 | !! to surface fluxes or mixing). These are intended to be called from call_tracer_column_fns | |
| 7 | !! in the MOM_tracer_flow_control module. | |
| 8 | module MOM_tracer_diabatic | |
| 9 | ||
| 10 | use MOM_grid, only : ocean_grid_type | |
| 11 | use MOM_verticalGrid, only : verticalGrid_type | |
| 12 | use MOM_forcing_type, only : forcing | |
| 13 | use MOM_error_handler, only : MOM_error, FATAL, WARNING | |
| 14 | ||
| 15 | implicit none ; private | |
| 16 | ||
| 17 | #include <MOM_memory.h> | |
| 18 | public tracer_vertdiff, tracer_vertdiff_Eulerian | |
| 19 | public applyTracerBoundaryFluxesInOut | |
| 20 | ||
| 21 | contains | |
| 22 | ||
| 23 | !> This subroutine solves a tridiagonal equation for the final tracer concentrations after the | |
| 24 | !! dual-entrainments, and possibly sinking or surface and bottom sources, are applied. The sinking | |
| 25 | !! is implemented with an fully implicit upwind advection scheme. Alternate time units can be | |
| 26 | !! used for the timestep, surface and bottom fluxes and sink_rate provided they are all consistent. | |
| 27 | 24 | subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, & |
| 28 | 0 | sfc_flux, btm_flux, btm_reservoir, sink_rate, convert_flux_in) |
| 29 | type(ocean_grid_type), intent(in) :: G !< ocean grid structure | |
| 30 | type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure | |
| 31 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_old !< layer thickness before entrainment | |
| 32 | !! [H ~> m or kg m-2] | |
| 33 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: ea !< amount of fluid entrained from the layer | |
| 34 | !! above [H ~> m or kg m-2] | |
| 35 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: eb !< amount of fluid entrained from the layer | |
| 36 | !! below [H ~> m or kg m-2] | |
| 37 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: tr !< tracer concentration in concentration units [CU] | |
| 38 | real, intent(in) :: dt !< amount of time covered by this call [T ~> s] | |
| 39 | real, dimension(SZI_(G),SZJ_(G)), optional,intent(in) :: sfc_flux !< surface flux of the tracer in units of | |
| 40 | !! [CU R Z T-1 ~> CU kg m-2 s-1] or | |
| 41 | !! [CU H ~> CU m or CU kg m-2] if | |
| 42 | !! convert_flux_in is .false. | |
| 43 | real, dimension(SZI_(G),SZJ_(G)), optional,intent(in) :: btm_flux !< The (negative upward) bottom flux of the | |
| 44 | !! tracer in [CU R Z T-1 ~> CU kg m-2 s-1] or | |
| 45 | !! [CU H ~> CU m or CU kg m-2] if | |
| 46 | !! convert_flux_in is .false. | |
| 47 | real, dimension(SZI_(G),SZJ_(G)), optional,intent(inout) :: btm_reservoir !< amount of tracer in a bottom reservoir | |
| 48 | !! [CU R Z ~> CU kg m-2] | |
| 49 | real, optional,intent(in) :: sink_rate !< rate at which the tracer sinks | |
| 50 | !! [Z T-1 ~> m s-1] | |
| 51 | logical, optional,intent(in) :: convert_flux_in !< True if the specified sfc_flux needs | |
| 52 | !! to be integrated in time | |
| 53 | ||
| 54 | ! local variables | |
| 55 | real :: sink_dist !< The distance the tracer sinks in a time step [H ~> m or kg m-2]. | |
| 56 | real, dimension(SZI_(G),SZJ_(G)) :: & | |
| 57 | 24 | sfc_src, & !< The time-integrated surface source of the tracer [CU H ~> CU m or CU kg m-2]. |
| 58 | 24 | btm_src !< The time-integrated bottom source of the tracer [CU H ~> CU m or CU kg m-2]. |
| 59 | real, dimension(SZI_(G)) :: & | |
| 60 | 12 | b1, & !< b1 is used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1]. |
| 61 | 24 | d1 !! d1=1-c1 is used by the tridiagonal solver [nondim]. |
| 62 | 24 | real :: c1(SZI_(G),SZK_(GV)) !< c1 is used by the tridiagonal solver [nondim]. |
| 63 | 24 | real :: h_minus_dsink(SZI_(G),SZK_(GV)) !< The layer thickness minus the |
| 64 | !! difference in sinking rates across the layer [H ~> m or kg m-2]. | |
| 65 | !! By construction, 0 <= h_minus_dsink < h_work. | |
| 66 | 12 | real :: sink(SZI_(G),SZK_(GV)+1) !< The tracer's sinking distances at the |
| 67 | !! interfaces, limited to prevent characteristics from | |
| 68 | !! crossing within a single timestep [H ~> m or kg m-2]. | |
| 69 | real :: b_denom_1 !< The first term in the denominator of b1 [H ~> m or kg m-2]. | |
| 70 | real :: h_tr !< h_tr is h at tracer points with a h_neglect added to | |
| 71 | !! ensure positive definiteness [H ~> m or kg m-2]. | |
| 72 | real :: h_neglect !< A thickness that is so small it is usually lost | |
| 73 | !! in roundoff and can be neglected [H ~> m or kg m-2]. | |
| 74 | logical :: convert_flux = .true. | |
| 75 | ||
| 76 | integer :: i, j, k, is, ie, js, je, nz | |
| 77 | 12 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke |
| 78 | ||
| 79 | 12 | if (nz == 1) then |
| 80 | call MOM_error(WARNING, "MOM_tracer_diabatic.F90, tracer_vertdiff called "//& | |
| 81 | 0 | "with only one vertical level") |
| 82 | 0 | return |
| 83 | endif | |
| 84 | ||
| 85 | 12 | if (present(convert_flux_in)) convert_flux = convert_flux_in |
| 86 | 12 | h_neglect = GV%H_subroundoff |
| 87 | 12 | sink_dist = 0.0 |
| 88 | 12 | if (present(sink_rate)) sink_dist = (dt*sink_rate) * GV%Z_to_H |
| 89 | !$OMP parallel default(shared) private(sink,h_minus_dsink,b_denom_1,b1,d1,h_tr,c1) | |
| 90 | !$OMP do | |
| 91 | 87132 | do j=js,je ; do i=is,ie ; sfc_src(i,j) = 0.0 ; btm_src(i,j) = 0.0 ; enddo ; enddo |
| 92 | 12 | if (present(sfc_flux)) then |
| 93 | 0 | if (convert_flux) then |
| 94 | !$OMP do | |
| 95 | 0 | do j=js,je ; do i=is,ie |
| 96 | 0 | sfc_src(i,j) = (sfc_flux(i,j)*dt) * GV%RZ_to_H |
| 97 | enddo ; enddo | |
| 98 | else | |
| 99 | !$OMP do | |
| 100 | 0 | do j=js,je ; do i=is,ie |
| 101 | 0 | sfc_src(i,j) = sfc_flux(i,j) |
| 102 | enddo ; enddo | |
| 103 | endif | |
| 104 | endif | |
| 105 | 12 | if (present(btm_flux)) then |
| 106 | 0 | if (convert_flux) then |
| 107 | !$OMP do | |
| 108 | 0 | do j=js,je ; do i=is,ie |
| 109 | 0 | btm_src(i,j) = (btm_flux(i,j)*dt) * GV%RZ_to_H |
| 110 | enddo ; enddo | |
| 111 | else | |
| 112 | !$OMP do | |
| 113 | 0 | do j=js,je ; do i=is,ie |
| 114 | 0 | btm_src(i,j) = btm_flux(i,j) |
| 115 | enddo ; enddo | |
| 116 | endif | |
| 117 | endif | |
| 118 | ||
| 119 | 12 | if (present(sink_rate)) then |
| 120 | !$OMP do | |
| 121 | 0 | do j=js,je |
| 122 | ! Find the sinking rates at all interfaces, limiting them if necesary | |
| 123 | ! so that the characteristics do not cross within a timestep. | |
| 124 | ! If a non-constant sinking rate were used, that would be incorprated | |
| 125 | ! here. | |
| 126 | 0 | if (present(btm_reservoir)) then |
| 127 | 0 | do i=is,ie ; sink(i,nz+1) = sink_dist ; enddo |
| 128 | 0 | do k=2,nz ; do i=is,ie |
| 129 | 0 | sink(i,K) = sink_dist ; h_minus_dsink(i,k) = h_old(i,j,k) |
| 130 | enddo ; enddo | |
| 131 | else | |
| 132 | 0 | do i=is,ie ; sink(i,nz+1) = 0.0 ; enddo |
| 133 | ! Find the limited sinking distance at the interfaces. | |
| 134 | 0 | do k=nz,2,-1 ; do i=is,ie |
| 135 | 0 | if (sink(i,K+1) >= sink_dist) then |
| 136 | 0 | sink(i,K) = sink_dist |
| 137 | 0 | h_minus_dsink(i,k) = h_old(i,j,k) + (sink(i,K+1) - sink(i,K)) |
| 138 | 0 | elseif (sink(i,K+1) + h_old(i,j,k) < sink_dist) then |
| 139 | 0 | sink(i,K) = sink(i,K+1) + h_old(i,j,k) |
| 140 | 0 | h_minus_dsink(i,k) = 0.0 |
| 141 | else | |
| 142 | 0 | sink(i,K) = sink_dist |
| 143 | 0 | h_minus_dsink(i,k) = (h_old(i,j,k) + sink(i,K+1)) - sink(i,K) |
| 144 | endif | |
| 145 | enddo ; enddo | |
| 146 | endif | |
| 147 | 0 | do i=is,ie |
| 148 | 0 | sink(i,1) = 0.0 ; h_minus_dsink(i,1) = (h_old(i,j,1) + sink(i,2)) |
| 149 | enddo | |
| 150 | ||
| 151 | ! Now solve the tridiagonal equation for the tracer concentrations. | |
| 152 | 0 | do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then |
| 153 | 0 | b_denom_1 = h_minus_dsink(i,1) + ea(i,j,1) + h_neglect |
| 154 | 0 | b1(i) = 1.0 / (b_denom_1 + eb(i,j,1)) |
| 155 | 0 | d1(i) = b_denom_1 * b1(i) |
| 156 | 0 | h_tr = h_old(i,j,1) + h_neglect |
| 157 | 0 | tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + b1(i)*sfc_src(i,j) |
| 158 | endif ; enddo | |
| 159 | 0 | do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then |
| 160 | 0 | c1(i,k) = eb(i,j,k-1) * b1(i) |
| 161 | b_denom_1 = h_minus_dsink(i,k) + d1(i) * (ea(i,j,k) + sink(i,K)) + & | |
| 162 | 0 | h_neglect |
| 163 | 0 | b1(i) = 1.0 / (b_denom_1 + eb(i,j,k)) |
| 164 | 0 | d1(i) = b_denom_1 * b1(i) |
| 165 | 0 | h_tr = h_old(i,j,k) + h_neglect |
| 166 | tr(i,j,k) = b1(i) * (h_tr * tr(i,j,k) + & | |
| 167 | 0 | (ea(i,j,k) + sink(i,K)) * tr(i,j,k-1)) |
| 168 | endif ; enddo ; enddo | |
| 169 | 0 | do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then |
| 170 | 0 | c1(i,nz) = eb(i,j,nz-1) * b1(i) |
| 171 | b_denom_1 = h_minus_dsink(i,nz) + d1(i) * (ea(i,j,nz) + sink(i,nz)) + & | |
| 172 | 0 | h_neglect |
| 173 | 0 | b1(i) = 1.0 / (b_denom_1 + eb(i,j,nz)) |
| 174 | 0 | h_tr = h_old(i,j,nz) + h_neglect |
| 175 | tr(i,j,nz) = b1(i) * ((h_tr * tr(i,j,nz) + btm_src(i,j)) + & | |
| 176 | 0 | (ea(i,j,nz) + sink(i,nz)) * tr(i,j,nz-1)) |
| 177 | endif ; enddo | |
| 178 | 0 | if (present(btm_reservoir)) then ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then |
| 179 | 0 | btm_reservoir(i,j) = btm_reservoir(i,j) + (sink(i,nz+1)*tr(i,j,nz)) * GV%H_to_RZ |
| 180 | endif ; enddo ; endif | |
| 181 | ||
| 182 | 0 | do k=nz-1,1,-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then |
| 183 | 0 | tr(i,j,k) = tr(i,j,k) + c1(i,k+1)*tr(i,j,k+1) |
| 184 | endif ; enddo ; enddo | |
| 185 | enddo | |
| 186 | else | |
| 187 | !$OMP do | |
| 188 | 732 | do j=js,je |
| 189 | 87120 | do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then |
| 190 | 60168 | h_tr = h_old(i,j,1) + h_neglect |
| 191 | 60168 | b_denom_1 = h_tr + ea(i,j,1) |
| 192 | 60168 | b1(i) = 1.0 / (b_denom_1 + eb(i,j,1)) |
| 193 | 60168 | d1(i) = h_tr * b1(i) |
| 194 | 60168 | tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + b1(i)*sfc_src(i,j) |
| 195 | endif ; enddo | |
| 196 | 6360480 | do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then |
| 197 | 4392264 | c1(i,k) = eb(i,j,k-1) * b1(i) |
| 198 | 4392264 | h_tr = h_old(i,j,k) + h_neglect |
| 199 | 4392264 | b_denom_1 = h_tr + d1(i) * ea(i,j,k) |
| 200 | 4392264 | b1(i) = 1.0 / (b_denom_1 + eb(i,j,k)) |
| 201 | 4392264 | d1(i) = b_denom_1 * b1(i) |
| 202 | 4392264 | tr(i,j,k) = b1(i) * (h_tr * tr(i,j,k) + ea(i,j,k) * tr(i,j,k-1)) |
| 203 | endif ; enddo ; enddo | |
| 204 | 87120 | do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then |
| 205 | 60168 | c1(i,nz) = eb(i,j,nz-1) * b1(i) |
| 206 | 60168 | h_tr = h_old(i,j,nz) + h_neglect |
| 207 | 60168 | b_denom_1 = h_tr + d1(i)*ea(i,j,nz) |
| 208 | 60168 | b1(i) = 1.0 / ( b_denom_1 + eb(i,j,nz)) |
| 209 | tr(i,j,nz) = b1(i) * (( h_tr * tr(i,j,nz) + btm_src(i,j)) + & | |
| 210 | 60168 | ea(i,j,nz) * tr(i,j,nz-1)) |
| 211 | endif ; enddo | |
| 212 | 6447612 | do k=nz-1,1,-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then |
| 213 | 4452432 | tr(i,j,k) = tr(i,j,k) + c1(i,k+1)*tr(i,j,k+1) |
| 214 | endif ; enddo ; enddo | |
| 215 | enddo | |
| 216 | endif | |
| 217 | !$OMP end parallel | |
| 218 | ||
| 219 | 12 | end subroutine tracer_vertdiff |
| 220 | ||
| 221 | ||
| 222 | !> This subroutine solves a tridiagonal equation for the final tracer concentrations after | |
| 223 | !! Eulerian mixing, and possibly sinking or surface and bottom sources, are applied. The sinking | |
| 224 | !! is implemented with an fully implicit upwind advection scheme. Alternate time units can be | |
| 225 | !! used for the timestep, surface and bottom fluxes and sink_rate provided they are all consistent. | |
| 226 | 0 | subroutine tracer_vertdiff_Eulerian(h_old, ent, dt, tr, G, GV, & |
| 227 | 0 | sfc_flux, btm_flux, btm_reservoir, sink_rate, convert_flux_in) |
| 228 | type(ocean_grid_type), intent(in) :: G !< ocean grid structure | |
| 229 | type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure | |
| 230 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_old !< layer thickness before entrainment | |
| 231 | !! [H ~> m or kg m-2] | |
| 232 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: ent !< Amount of fluid mixed across interfaces | |
| 233 | !! [H ~> m or kg m-2] | |
| 234 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: tr !< tracer concentration in concentration units [CU] | |
| 235 | real, intent(in) :: dt !< amount of time covered by this call [T ~> s] | |
| 236 | real, dimension(SZI_(G),SZJ_(G)), optional,intent(in) :: sfc_flux !< surface flux of the tracer in units of | |
| 237 | !! [CU R Z T-1 ~> CU kg m-2 s-1] or | |
| 238 | !! [CU H ~> CU m or CU kg m-2] if | |
| 239 | !! convert_flux_in is .false. | |
| 240 | real, dimension(SZI_(G),SZJ_(G)), optional,intent(in) :: btm_flux !< The (negative upward) bottom flux of the | |
| 241 | !! tracer in [CU kg m-2 T-1 ~> CU kg m-2 s-1] or | |
| 242 | !! [CU H ~> CU m or CU kg m-2] if | |
| 243 | !! convert_flux_in is .false. | |
| 244 | real, dimension(SZI_(G),SZJ_(G)), optional,intent(inout) :: btm_reservoir !< amount of tracer in a bottom reservoir | |
| 245 | !! [CU R Z ~> CU kg m-2] | |
| 246 | real, optional,intent(in) :: sink_rate !< rate at which the tracer sinks | |
| 247 | !! [Z T-1 ~> m s-1] | |
| 248 | logical, optional,intent(in) :: convert_flux_in !< True if the specified sfc_flux needs | |
| 249 | !! to be integrated in time | |
| 250 | ||
| 251 | ! local variables | |
| 252 | real :: sink_dist !< The distance the tracer sinks in a time step [H ~> m or kg m-2]. | |
| 253 | real, dimension(SZI_(G),SZJ_(G)) :: & | |
| 254 | 0 | sfc_src, & !< The time-integrated surface source of the tracer [CU H ~> CU m or CU kg m-2]. |
| 255 | 0 | btm_src !< The time-integrated bottom source of the tracer [CU H ~> CU m or CU kg m-2]. |
| 256 | real, dimension(SZI_(G)) :: & | |
| 257 | 0 | b1, & !< b1 is used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1]. |
| 258 | 0 | d1 !! d1=1-c1 is used by the tridiagonal solver [nondim]. |
| 259 | 0 | real :: c1(SZI_(G),SZK_(GV)) !< c1 is used by the tridiagonal solver [nondim]. |
| 260 | 0 | real :: h_minus_dsink(SZI_(G),SZK_(GV)) !< The layer thickness minus the |
| 261 | !! difference in sinking rates across the layer [H ~> m or kg m-2]. | |
| 262 | !! By construction, 0 <= h_minus_dsink < h_work. | |
| 263 | 0 | real :: sink(SZI_(G),SZK_(GV)+1) !< The tracer's sinking distances at the |
| 264 | !! interfaces, limited to prevent characteristics from | |
| 265 | !! crossing within a single timestep [H ~> m or kg m-2]. | |
| 266 | real :: b_denom_1 !< The first term in the denominator of b1 [H ~> m or kg m-2]. | |
| 267 | real :: h_tr !< h_tr is h at tracer points with a h_neglect added to | |
| 268 | !! ensure positive definiteness [H ~> m or kg m-2]. | |
| 269 | real :: h_neglect !< A thickness that is so small it is usually lost | |
| 270 | !! in roundoff and can be neglected [H ~> m or kg m-2]. | |
| 271 | logical :: convert_flux | |
| 272 | ||
| 273 | integer :: i, j, k, is, ie, js, je, nz | |
| 274 | 0 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke |
| 275 | ||
| 276 | 0 | if (nz == 1) then |
| 277 | call MOM_error(WARNING, "MOM_tracer_diabatic.F90, tracer_vertdiff called "//& | |
| 278 | 0 | "with only one vertical level") |
| 279 | 0 | return |
| 280 | endif | |
| 281 | ||
| 282 | 0 | convert_flux = .true. |
| 283 | 0 | if (present(convert_flux_in)) convert_flux = convert_flux_in |
| 284 | 0 | h_neglect = GV%H_subroundoff |
| 285 | 0 | sink_dist = 0.0 |
| 286 | 0 | if (present(sink_rate)) sink_dist = (dt*sink_rate) * GV%Z_to_H |
| 287 | !$OMP parallel default(shared) private(sink,h_minus_dsink,b_denom_1,b1,d1,h_tr,c1) | |
| 288 | !$OMP do | |
| 289 | 0 | do j=js,je ; do i=is,ie ; sfc_src(i,j) = 0.0 ; btm_src(i,j) = 0.0 ; enddo ; enddo |
| 290 | 0 | if (present(sfc_flux)) then |
| 291 | 0 | if (convert_flux) then |
| 292 | !$OMP do | |
| 293 | 0 | do j=js,je ; do i=is,ie |
| 294 | 0 | sfc_src(i,j) = (sfc_flux(i,j)*dt) * GV%RZ_to_H |
| 295 | enddo ; enddo | |
| 296 | else | |
| 297 | !$OMP do | |
| 298 | 0 | do j=js,je ; do i=is,ie |
| 299 | 0 | sfc_src(i,j) = sfc_flux(i,j) |
| 300 | enddo ; enddo | |
| 301 | endif | |
| 302 | endif | |
| 303 | 0 | if (present(btm_flux)) then |
| 304 | 0 | if (convert_flux) then |
| 305 | !$OMP do | |
| 306 | 0 | do j=js,je ; do i=is,ie |
| 307 | 0 | btm_src(i,j) = (btm_flux(i,j)*dt) * GV%kg_m2_to_H |
| 308 | enddo ; enddo | |
| 309 | else | |
| 310 | !$OMP do | |
| 311 | 0 | do j=js,je ; do i=is,ie |
| 312 | 0 | btm_src(i,j) = btm_flux(i,j) |
| 313 | enddo ; enddo | |
| 314 | endif | |
| 315 | endif | |
| 316 | ||
| 317 | 0 | if (present(sink_rate)) then |
| 318 | !$OMP do | |
| 319 | 0 | do j=js,je |
| 320 | ! Find the sinking rates at all interfaces, limiting them if necesary | |
| 321 | ! so that the characteristics do not cross within a timestep. | |
| 322 | ! If a non-constant sinking rate were used, that would be incorprated | |
| 323 | ! here. | |
| 324 | 0 | if (present(btm_reservoir)) then |
| 325 | 0 | do i=is,ie ; sink(i,nz+1) = sink_dist ; enddo |
| 326 | 0 | do k=2,nz ; do i=is,ie |
| 327 | 0 | sink(i,K) = sink_dist ; h_minus_dsink(i,k) = h_old(i,j,k) |
| 328 | enddo ; enddo | |
| 329 | else | |
| 330 | 0 | do i=is,ie ; sink(i,nz+1) = 0.0 ; enddo |
| 331 | ! Find the limited sinking distance at the interfaces. | |
| 332 | 0 | do k=nz,2,-1 ; do i=is,ie |
| 333 | 0 | if (sink(i,K+1) >= sink_dist) then |
| 334 | 0 | sink(i,K) = sink_dist |
| 335 | 0 | h_minus_dsink(i,k) = h_old(i,j,k) + (sink(i,K+1) - sink(i,K)) |
| 336 | 0 | elseif (sink(i,K+1) + h_old(i,j,k) < sink_dist) then |
| 337 | 0 | sink(i,K) = sink(i,K+1) + h_old(i,j,k) |
| 338 | 0 | h_minus_dsink(i,k) = 0.0 |
| 339 | else | |
| 340 | 0 | sink(i,K) = sink_dist |
| 341 | 0 | h_minus_dsink(i,k) = (h_old(i,j,k) + sink(i,K+1)) - sink(i,K) |
| 342 | endif | |
| 343 | enddo ; enddo | |
| 344 | endif | |
| 345 | 0 | do i=is,ie |
| 346 | 0 | sink(i,1) = 0.0 ; h_minus_dsink(i,1) = (h_old(i,j,1) + sink(i,2)) |
| 347 | enddo | |
| 348 | ||
| 349 | ! Now solve the tridiagonal equation for the tracer concentrations. | |
| 350 | 0 | do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then |
| 351 | 0 | b_denom_1 = h_minus_dsink(i,1) + ent(i,j,1) + h_neglect |
| 352 | 0 | b1(i) = 1.0 / (b_denom_1 + ent(i,j,2)) |
| 353 | 0 | d1(i) = b_denom_1 * b1(i) |
| 354 | 0 | h_tr = h_old(i,j,1) + h_neglect |
| 355 | 0 | tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + b1(i)*sfc_src(i,j) |
| 356 | endif ; enddo | |
| 357 | 0 | do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then |
| 358 | 0 | c1(i,k) = ent(i,j,K) * b1(i) |
| 359 | b_denom_1 = h_minus_dsink(i,k) + d1(i) * (ent(i,j,K) + sink(i,K)) + & | |
| 360 | 0 | h_neglect |
| 361 | 0 | b1(i) = 1.0 / (b_denom_1 + ent(i,j,K+1)) |
| 362 | 0 | d1(i) = b_denom_1 * b1(i) |
| 363 | 0 | h_tr = h_old(i,j,k) + h_neglect |
| 364 | tr(i,j,k) = b1(i) * (h_tr * tr(i,j,k) + & | |
| 365 | 0 | (ent(i,j,K) + sink(i,K)) * tr(i,j,k-1)) |
| 366 | endif ; enddo ; enddo | |
| 367 | 0 | do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then |
| 368 | 0 | c1(i,nz) = ent(i,j,nz) * b1(i) |
| 369 | b_denom_1 = h_minus_dsink(i,nz) + d1(i) * (ent(i,j,nz) + sink(i,nz)) + & | |
| 370 | 0 | h_neglect |
| 371 | 0 | b1(i) = 1.0 / (b_denom_1 + ent(i,j,nz+1)) |
| 372 | 0 | h_tr = h_old(i,j,nz) + h_neglect |
| 373 | tr(i,j,nz) = b1(i) * ((h_tr * tr(i,j,nz) + btm_src(i,j)) + & | |
| 374 | 0 | (ent(i,j,nz) + sink(i,nz)) * tr(i,j,nz-1)) |
| 375 | endif ; enddo | |
| 376 | 0 | if (present(btm_reservoir)) then ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then |
| 377 | 0 | btm_reservoir(i,j) = btm_reservoir(i,j) + (sink(i,nz+1)*tr(i,j,nz)) * GV%H_to_RZ |
| 378 | endif ; enddo ; endif | |
| 379 | ||
| 380 | 0 | do k=nz-1,1,-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then |
| 381 | 0 | tr(i,j,k) = tr(i,j,k) + c1(i,k+1)*tr(i,j,k+1) |
| 382 | endif ; enddo ; enddo | |
| 383 | enddo | |
| 384 | else | |
| 385 | !$OMP do | |
| 386 | 0 | do j=js,je |
| 387 | 0 | do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then |
| 388 | 0 | h_tr = h_old(i,j,1) + h_neglect |
| 389 | 0 | b_denom_1 = h_tr + ent(i,j,1) |
| 390 | 0 | b1(i) = 1.0 / (b_denom_1 + ent(i,j,2)) |
| 391 | 0 | d1(i) = h_tr * b1(i) |
| 392 | 0 | tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + b1(i)*sfc_src(i,j) |
| 393 | endif ; enddo | |
| 394 | 0 | do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then |
| 395 | 0 | c1(i,k) = ent(i,j,K) * b1(i) |
| 396 | 0 | h_tr = h_old(i,j,k) + h_neglect |
| 397 | 0 | b_denom_1 = h_tr + d1(i) * ent(i,j,K) |
| 398 | 0 | b1(i) = 1.0 / (b_denom_1 + ent(i,j,K+1)) |
| 399 | 0 | d1(i) = b_denom_1 * b1(i) |
| 400 | 0 | tr(i,j,k) = b1(i) * (h_tr * tr(i,j,k) + ent(i,j,K) * tr(i,j,k-1)) |
| 401 | endif ; enddo ; enddo | |
| 402 | 0 | do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then |
| 403 | 0 | c1(i,nz) = ent(i,j,nz) * b1(i) |
| 404 | 0 | h_tr = h_old(i,j,nz) + h_neglect |
| 405 | 0 | b_denom_1 = h_tr + d1(i)*ent(i,j,nz) |
| 406 | 0 | b1(i) = 1.0 / ( b_denom_1 + ent(i,j,nz+1)) |
| 407 | tr(i,j,nz) = b1(i) * (( h_tr * tr(i,j,nz) + btm_src(i,j)) + & | |
| 408 | 0 | ent(i,j,nz) * tr(i,j,nz-1)) |
| 409 | endif ; enddo | |
| 410 | 0 | do k=nz-1,1,-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then |
| 411 | 0 | tr(i,j,k) = tr(i,j,k) + c1(i,k+1)*tr(i,j,k+1) |
| 412 | endif ; enddo ; enddo | |
| 413 | enddo | |
| 414 | endif | |
| 415 | !$OMP end parallel | |
| 416 | ||
| 417 | 0 | end subroutine tracer_vertdiff_Eulerian |
| 418 | ||
| 419 | ||
| 420 | !> This routine is modeled after applyBoundaryFluxesInOut in MOM_diabatic_aux.F90 | |
| 421 | !! NOTE: Please note that in this routine sfc_flux gets set to zero to ensure that the surface | |
| 422 | !! flux of the tracer does not get applied again during a subsequent call to tracer_vertdif | |
| 423 | 24 | subroutine applyTracerBoundaryFluxesInOut(G, GV, Tr, dt, fluxes, h, evap_CFL_limit, minimum_forcing_depth, & |
| 424 | 0 | in_flux_optional, out_flux_optional, update_h_opt) |
| 425 | ||
| 426 | type(ocean_grid_type), intent(in ) :: G !< Grid structure | |
| 427 | type(verticalGrid_type), intent(in ) :: GV !< ocean vertical grid structure | |
| 428 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: Tr !< Tracer concentration on T-cell [conc] | |
| 429 | real, intent(in ) :: dt !< Time-step over which forcing is applied [T ~> s] | |
| 430 | type(forcing), intent(in ) :: fluxes !< Surface fluxes container | |
| 431 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2] | |
| 432 | real, intent(in ) :: evap_CFL_limit !< Limit on the fraction of the | |
| 433 | !! water that can be fluxed out of the top | |
| 434 | !! layer in a timestep [nondim] | |
| 435 | real, intent(in ) :: minimum_forcing_depth !< The smallest depth over | |
| 436 | !! which fluxes can be applied [H ~> m or kg m-2] | |
| 437 | real, dimension(SZI_(G),SZJ_(G)), optional, intent(in ) :: in_flux_optional !< The total time-integrated | |
| 438 | !! amount of tracer that enters with freshwater | |
| 439 | !! [conc H ~> conc m or conc kg m-2] | |
| 440 | real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: out_flux_optional !< The total time-integrated | |
| 441 | !! amount of tracer that leaves with freshwater | |
| 442 | !! [conc H ~> conc m or conc kg m-2] | |
| 443 | logical, optional, intent(in) :: update_h_opt !< Optional flag to determine whether | |
| 444 | !! h should be updated | |
| 445 | ||
| 446 | integer, parameter :: maxGroundings = 5 | |
| 447 | integer :: numberOfGroundings, iGround(maxGroundings), jGround(maxGroundings) | |
| 448 | real :: IforcingDepthScale ! The inverse of the scale over which to apply forcing [H-1 ~> m-1 or m2 kg-1] | |
| 449 | real :: dThickness ! The change in a layer's thickness [H ~> m or kg m-2] | |
| 450 | real :: dTracer ! The change in the integrated tracer content of a layer [conc H ~> conc m or conc kg m-2] | |
| 451 | real :: fractionOfForcing ! The fraction of the forcing to apply to a layer [nondim] | |
| 452 | real :: hOld ! The layer thickness before surface forcing is applied [H ~> m or kg m-2] | |
| 453 | real :: Ithickness ! The inverse of the new layer thickness [H-1 ~> m-1 or m2 kg-1] | |
| 454 | ||
| 455 | 12 | real :: h2d(SZI_(G),SZK_(GV)) ! A 2-d work copy of layer thicknesses [H ~> m or kg m-2] |
| 456 | 24 | real :: Tr2d(SZI_(G),SZK_(GV)) ! A 2-d work copy of tracer concentrations [conc] |
| 457 | 24 | real :: in_flux(SZI_(G),SZJ_(G)) ! The total time-integrated amount of tracer that |
| 458 | ! enters with freshwater [conc H ~> conc m or conc kg m-2] | |
| 459 | 24 | real :: out_flux(SZI_(G),SZJ_(G)) ! The total time-integrated amount of tracer that |
| 460 | ! leaves with freshwater [conc H ~> conc m or conc kg m-2] | |
| 461 | 24 | real :: netMassIn(SZI_(G)) ! The remaining mass entering ocean surface [H ~> m or kg m-2] |
| 462 | 24 | real :: netMassOut(SZI_(G)) ! The remaining mass leaving ocean surface [H ~> m or kg m-2] |
| 463 | 24 | real :: in_flux_1d(SZI_(G)) ! The remaining amount of tracer that enters with |
| 464 | ! the freshwater [conc H ~> conc m or conc kg m-2] | |
| 465 | 12 | real :: out_flux_1d(SZI_(G)) ! The remaining amount of tracer that leaves with |
| 466 | ! the freshwater [conc H ~> conc m or conc kg m-2] | |
| 467 | real :: hGrounding(maxGroundings) ! The remaining fresh water flux that was not able to be | |
| 468 | ! supplied from a column that grounded out [H ~> m or kg m-2] | |
| 469 | logical :: update_h | |
| 470 | integer :: i, j, is, ie, js, je, k, nz | |
| 471 | character(len=45) :: mesg | |
| 472 | ||
| 473 | 12 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke |
| 474 | ||
| 475 | ! If no freshwater fluxes, nothing needs to be done in this routine | |
| 476 | 12 | if ( (.not. associated(fluxes%netMassIn)) .or. (.not. associated(fluxes%netMassOut)) ) return |
| 477 | ||
| 478 | 210540 | in_flux(:,:) = 0.0 ; out_flux(:,:) = 0.0 |
| 479 | 12 | if (present(in_flux_optional)) then |
| 480 | 0 | do j=js,je ; do i=is,ie |
| 481 | 0 | in_flux(i,j) = in_flux_optional(i,j) |
| 482 | enddo ; enddo | |
| 483 | endif | |
| 484 | 12 | if (present(out_flux_optional)) then |
| 485 | 0 | do j=js,je ; do i=is,ie |
| 486 | 0 | out_flux(i,j) = out_flux_optional(i,j) |
| 487 | enddo ; enddo | |
| 488 | endif | |
| 489 | ||
| 490 | 12 | if (present(update_h_opt)) then |
| 491 | 0 | update_h = update_h_opt |
| 492 | else | |
| 493 | 12 | update_h = .true. |
| 494 | endif | |
| 495 | ||
| 496 | 12 | numberOfGroundings = 0 |
| 497 | ||
| 498 | !$OMP parallel do default(none) shared(is,ie,js,je,nz,h,Tr,G,GV,fluxes,dt, & | |
| 499 | !$OMP IforcingDepthScale,minimum_forcing_depth, & | |
| 500 | !$OMP numberOfGroundings,iGround,jGround,update_h, & | |
| 501 | !$OMP in_flux,out_flux,hGrounding,evap_CFL_limit) & | |
| 502 | !$OMP private(h2d,Tr2d,netMassIn,netMassOut, & | |
| 503 | !$OMP in_flux_1d,out_flux_1d,fractionOfForcing, & | |
| 504 | !$OMP dThickness,dTracer,hOld,Ithickness) | |
| 505 | ||
| 506 | ! Work in vertical slices for efficiency | |
| 507 | 732 | do j=js,je |
| 508 | ||
| 509 | ! Copy state into 2D-slice arrays | |
| 510 | 6534720 | do k=1,nz ; do i=is,ie |
| 511 | 6480000 | h2d(i,k) = h(i,j,k) |
| 512 | 6534000 | Tr2d(i,k) = Tr(i,j,k) |
| 513 | enddo ; enddo | |
| 514 | ||
| 515 | 87120 | do i = is,ie |
| 516 | 86400 | in_flux_1d(i) = in_flux(i,j) |
| 517 | 87120 | out_flux_1d(i) = out_flux(i,j) |
| 518 | enddo | |
| 519 | ! The surface forcing is contained in the fluxes type. | |
| 520 | ! We aggregate the thermodynamic forcing for a time step into the following: | |
| 521 | ! These should have been set and stored during a call to applyBoundaryFluxesInOut | |
| 522 | ! netMassIn = net mass entering at ocean surface over a timestep | |
| 523 | ! netMassOut = net mass leaving ocean surface [H ~> m or kg m-2] over a time step. | |
| 524 | ! netMassOut < 0 means mass leaves ocean. | |
| 525 | ||
| 526 | ! Note here that the aggregateFW flag has already been taken care of in the call to | |
| 527 | ! applyBoundaryFluxesInOut | |
| 528 | 87120 | do i=is,ie |
| 529 | 86400 | netMassOut(i) = fluxes%netMassOut(i,j) |
| 530 | 87120 | netMassIn(i) = fluxes%netMassIn(i,j) |
| 531 | enddo | |
| 532 | ||
| 533 | ! Apply the surface boundary fluxes in three steps: | |
| 534 | ! A/ update concentration from mass entering the ocean | |
| 535 | ! B/ update concentration from mass leaving ocean. | |
| 536 | 87120 | do i=is,ie |
| 537 | 86400 | if (G%mask2dT(i,j)>0.) then |
| 538 | ||
| 539 | ! A/ Update tracer due to incoming mass flux. | |
| 540 | 120336 | do k=1,1 |
| 541 | ||
| 542 | ! Change in state due to forcing | |
| 543 | 60168 | dThickness = netMassIn(i) ! Since we are adding mass, we can use all of it |
| 544 | 60168 | dTracer = 0. |
| 545 | ||
| 546 | ! Update the forcing by the part to be consumed within the present k-layer. | |
| 547 | ! If fractionOfForcing = 1, then updated netMassIn, netHeat, and netSalt vanish. | |
| 548 | 60168 | netMassIn(i) = netMassIn(i) - dThickness |
| 549 | 60168 | dTracer = dTracer + in_flux_1d(i) |
| 550 | 60168 | in_flux_1d(i) = 0.0 |
| 551 | ||
| 552 | ! Update state | |
| 553 | 60168 | hOld = h2d(i,k) ! Keep original thickness in hand |
| 554 | 60168 | h2d(i,k) = h2d(i,k) + dThickness ! New thickness |
| 555 | 120336 | if (h2d(i,k) > 0.0) then |
| 556 | 60168 | Ithickness = 1.0/h2d(i,k) ! Inverse new thickness |
| 557 | ! The "if"s below avoid changing T/S by roundoff unnecessarily | |
| 558 | 60168 | if (dThickness /= 0. .or. dTracer /= 0.) tr2d(i,k) = (hOld*tr2d(i,k)+ dTracer)*Ithickness |
| 559 | endif | |
| 560 | ||
| 561 | enddo ! k=1,1 | |
| 562 | ||
| 563 | ! B/ Update tracer from mass leaving ocean | |
| 564 | 4572768 | do k=1,nz |
| 565 | ||
| 566 | ! Place forcing into this layer if this layer has nontrivial thickness. | |
| 567 | ! For layers thin relative to 1/IforcingDepthScale, then distribute | |
| 568 | ! forcing into deeper layers. | |
| 569 | 4512600 | IforcingDepthScale = 1. / max(GV%H_subroundoff, minimum_forcing_depth - netMassOut(i) ) |
| 570 | ! fractionOfForcing = 1.0, unless h2d is less than IforcingDepthScale. | |
| 571 | 4512600 | fractionOfForcing = min(1.0, h2d(i,k)*IforcingDepthScale) |
| 572 | ||
| 573 | ! In the case with (-1)*netMassOut*fractionOfForcing greater than cfl*h, we | |
| 574 | ! limit the forcing applied to this cell, leaving the remaining forcing to | |
| 575 | ! be distributed downwards. | |
| 576 | 4512600 | if (-fractionOfForcing*netMassOut(i) > evap_CFL_limit*h2d(i,k)) then |
| 577 | 0 | fractionOfForcing = -evap_CFL_limit*h2d(i,k)/netMassOut(i) |
| 578 | endif | |
| 579 | ||
| 580 | ! Change in state due to forcing | |
| 581 | 4512600 | dThickness = max( fractionOfForcing*netMassOut(i), -h2d(i,k) ) |
| 582 | ! Note this is slightly different to how salt is currently treated | |
| 583 | 4512600 | dTracer = fractionOfForcing*out_flux_1d(i) |
| 584 | ||
| 585 | ! Update the forcing by the part to be consumed within the present k-layer. | |
| 586 | ! If fractionOfForcing = 1, then new netMassOut vanishes. | |
| 587 | 4512600 | netMassOut(i) = netMassOut(i) - dThickness |
| 588 | 4512600 | out_flux_1d(i) = out_flux_1d(i) - dTracer |
| 589 | ||
| 590 | ! Update state by the appropriate increment. | |
| 591 | 4512600 | hOld = h2d(i,k) ! Keep original thickness in hand |
| 592 | 4512600 | h2d(i,k) = h2d(i,k) + dThickness ! New thickness |
| 593 | 4572768 | if (h2d(i,k) > 0.) then |
| 594 | 4512600 | Ithickness = 1.0/h2d(i,k) ! Inverse of new thickness |
| 595 | 4512600 | Tr2d(i,k) = (hOld*Tr2d(i,k) + dTracer)*Ithickness |
| 596 | endif | |
| 597 | ||
| 598 | enddo ! k | |
| 599 | ||
| 600 | endif | |
| 601 | ||
| 602 | ! If anything remains after the k-loop, then we have grounded out, which is a problem. | |
| 603 | 87120 | if (abs(in_flux_1d(i))+abs(out_flux_1d(i)) /= 0.0) then |
| 604 | !$OMP critical | |
| 605 | 0 | numberOfGroundings = numberOfGroundings +1 |
| 606 | 0 | if (numberOfGroundings<=maxGroundings) then |
| 607 | 0 | iGround(numberOfGroundings) = i ! Record i,j location of event for |
| 608 | 0 | jGround(numberOfGroundings) = j ! warning message |
| 609 | 0 | hGrounding(numberOfGroundings) = abs(in_flux_1d(i))+abs(out_flux_1d(i)) |
| 610 | endif | |
| 611 | !$OMP end critical | |
| 612 | endif | |
| 613 | ||
| 614 | enddo ! i | |
| 615 | ||
| 616 | ! Step C/ copy updated tracer concentration from the 2d slice now back into model state. | |
| 617 | 6534720 | do k=1,nz ; do i=is,ie |
| 618 | 6534000 | Tr(i,j,k) = Tr2d(i,k) |
| 619 | enddo ; enddo | |
| 620 | ||
| 621 | 732 | if (update_h) then |
| 622 | 6534720 | do k=1,nz ; do i=is,ie |
| 623 | 6534000 | h(i,j,k) = h2d(i,k) |
| 624 | enddo ; enddo | |
| 625 | endif | |
| 626 | ||
| 627 | enddo ! j-loop finish | |
| 628 | ||
| 629 | 12 | if (numberOfGroundings>0) then |
| 630 | 0 | do i = 1, min(numberOfGroundings, maxGroundings) |
| 631 | 0 | write(mesg(1:45),'(3es15.3)') G%geoLonT( iGround(i), jGround(i) ), & |
| 632 | 0 | G%geoLatT( iGround(i), jGround(i)) , hGrounding(i) |
| 633 | call MOM_error(WARNING, "MOM_tracer_diabatic.F90, applyTracerBoundaryFluxesInOut(): "//& | |
| 634 | 0 | "Tracer created. x,y,dh= "//trim(mesg), all_print=.true.) |
| 635 | enddo | |
| 636 | ||
| 637 | 0 | if (numberOfGroundings - maxGroundings > 0) then |
| 638 | 0 | write(mesg, '(I0)') numberOfGroundings - maxGroundings |
| 639 | call MOM_error(WARNING, "MOM_tracer_vertical.F90, applyTracerBoundaryFluxesInOut(): "//& | |
| 640 | 0 | trim(mesg) // " groundings remaining", all_print=.true.) |
| 641 | endif | |
| 642 | endif | |
| 643 | ||
| 644 | 12 | end subroutine applyTracerBoundaryFluxesInOut |
| 645 | end module MOM_tracer_diabatic |