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 | #include <do_concurrent_compat.h> | |
| 6 | ||
| 7 | !> This module contains the subroutines that advect tracers along coordinate surfaces. | |
| 8 | module MOM_tracer_advect | |
| 9 | ||
| 10 | use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end | |
| 11 | use MOM_cpu_clock, only : CLOCK_MODULE, CLOCK_ROUTINE | |
| 12 | use MOM_diag_mediator, only : post_data, query_averaging_enabled, diag_ctrl | |
| 13 | use MOM_diag_mediator, only : register_diag_field, safe_alloc_ptr, time_type | |
| 14 | use MOM_domains, only : sum_across_PEs, max_across_PEs | |
| 15 | use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type, pass_var | |
| 16 | use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg, is_root_pe | |
| 17 | use MOM_file_parser, only : get_param, log_version, param_file_type | |
| 18 | use MOM_grid, only : ocean_grid_type | |
| 19 | use MOM_open_boundary, only : ocean_OBC_type, OBC_NONE, OBC_DIRECTION_E | |
| 20 | use MOM_open_boundary, only : OBC_DIRECTION_W, OBC_DIRECTION_N, OBC_DIRECTION_S | |
| 21 | use MOM_open_boundary, only : OBC_segment_type | |
| 22 | use MOM_tracer_registry, only : tracer_registry_type, tracer_type | |
| 23 | use MOM_unit_scaling, only : unit_scale_type | |
| 24 | use MOM_verticalGrid, only : verticalGrid_type | |
| 25 | use MOM_tracer_advect_schemes, only : ADVECT_PLM, ADVECT_PPMH3, ADVECT_PPM | |
| 26 | use MOM_tracer_advect_schemes, only : set_tracer_advect_scheme, TracerAdvectionSchemeDoc | |
| 27 | implicit none ; private | |
| 28 | ||
| 29 | #include <MOM_memory.h> | |
| 30 | ||
| 31 | public advect_tracer | |
| 32 | public tracer_advect_init | |
| 33 | public tracer_advect_end | |
| 34 | ||
| 35 | !> Control structure for this module | |
| 36 | type, public :: tracer_advect_CS ; private | |
| 37 | real :: dt !< The baroclinic dynamics time step [T ~> s]. | |
| 38 | type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the | |
| 39 | !< timing of diagnostic output. | |
| 40 | logical :: debug !< If true, write verbose checksums for debugging purposes. | |
| 41 | logical :: useHuynhStencilBug = .false. !< If true, use the incorrect stencil width. | |
| 42 | !! This is provided for compatibility with legacy simuations. | |
| 43 | type(group_pass_type) :: pass_uhr_vhr_t_hprev !< A structure used for group passes | |
| 44 | integer :: default_advect_scheme = -1 !< Determines which reconstruction to use | |
| 45 | end type tracer_advect_CS | |
| 46 | ||
| 47 | !>@{ CPU time clocks | |
| 48 | integer :: id_clock_advect | |
| 49 | integer :: id_clock_pass | |
| 50 | integer :: id_clock_sync | |
| 51 | !>@} | |
| 52 | ||
| 53 | contains | |
| 54 | ||
| 55 | !> This routine time steps the tracer concentration using a | |
| 56 | !! monotonic, conservative, weakly diffusive scheme. | |
| 57 | 24 | subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, x_first_in, & |
| 58 | 0 | vol_prev, max_iter_in, update_vol_prev, uhr_out, vhr_out) |
| 59 | type(ocean_grid_type), intent(inout) :: G !< ocean grid structure | |
| 60 | type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure | |
| 61 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & | |
| 62 | intent(in) :: h_end !< Layer thickness after advection [H ~> m or kg m-2] | |
| 63 | real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & | |
| 64 | intent(in) :: uhtr !< Accumulated volume or mass flux through the | |
| 65 | !! zonal faces [H L2 ~> m3 or kg] | |
| 66 | real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & | |
| 67 | intent(in) :: vhtr !< Accumulated volume or mass flux through the | |
| 68 | !! meridional faces [H L2 ~> m3 or kg] | |
| 69 | type(ocean_OBC_type), pointer :: OBC !< specifies whether, where, and what OBCs are used | |
| 70 | real, intent(in) :: dt !< time increment [T ~> s] | |
| 71 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 72 | type(tracer_advect_CS), pointer :: CS !< control structure for module | |
| 73 | type(tracer_registry_type), pointer :: Reg !< pointer to tracer registry | |
| 74 | logical, optional, intent(in) :: x_first_in !< If present, indicate whether to update | |
| 75 | !! first in the x- or y-direction. | |
| 76 | ! The remaining optional arguments are only used in offline tracer mode. | |
| 77 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & | |
| 78 | optional, intent(inout) :: vol_prev !< Cell volume before advection [H L2 ~> m3 or kg]. | |
| 79 | !! If update_vol_prev is true, the returned value is | |
| 80 | !! the cell volume after the transport that was done | |
| 81 | !! by this call, and if all the transport could be | |
| 82 | !! accommodated it should be close to h_end*G%areaT. | |
| 83 | integer, optional, intent(in) :: max_iter_in !< The maximum number of iterations | |
| 84 | logical, optional, intent(in) :: update_vol_prev !< If present and true, update vol_prev to | |
| 85 | !! return its value after the tracer have been updated. | |
| 86 | real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & | |
| 87 | optional, intent(out) :: uhr_out !< Remaining accumulated volume or mass fluxes | |
| 88 | !! through the zonal faces [H L2 ~> m3 or kg] | |
| 89 | real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & | |
| 90 | optional, intent(out) :: vhr_out !< Remaining accumulated volume or mass fluxes | |
| 91 | !! through the meridional faces [H L2 ~> m3 or kg] | |
| 92 | ||
| 93 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: & | |
| 94 | 24 | hprev ! cell volume at the end of previous tracer change [H L2 ~> m3 or kg] |
| 95 | real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: & | |
| 96 | 24 | uhr ! The remaining zonal thickness flux [H L2 ~> m3 or kg] |
| 97 | real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: & | |
| 98 | 24 | vhr ! The remaining meridional thickness fluxes [H L2 ~> m3 or kg] |
| 99 | 24 | real :: uh_neglect(SZIB_(G),SZJ_(G)) ! uh_neglect and vh_neglect are the |
| 100 | 24 | real :: vh_neglect(SZI_(G),SZJB_(G)) ! magnitude of remaining transports that |
| 101 | ! can be simply discarded [H L2 ~> m3 or kg]. | |
| 102 | ||
| 103 | real :: landvolfill ! An arbitrary? nonzero cell volume [H L2 ~> m3 or kg]. | |
| 104 | logical :: use_PPM_stencil ! If true, use the correct PPM stencil width. | |
| 105 | real :: Idt ! 1/dt [T-1 ~> s-1]. | |
| 106 | 24 | logical :: domore_u(SZJ_(G),SZK_(GV)) ! domore_u and domore_v indicate whether there is more |
| 107 | 24 | logical :: domore_v(SZJB_(G),SZK_(GV)) ! advection to be done in the corresponding row or column. |
| 108 | logical :: x_first ! If true, advect in the x-direction first. | |
| 109 | integer :: max_iter ! maximum number of iterations in each layer | |
| 110 | 12 | integer :: domore_k(SZK_(GV)) |
| 111 | integer :: stencil ! stencil of the advection scheme | |
| 112 | integer :: nsten_halo ! number of stencils that fit in the halos | |
| 113 | integer :: i, j, k, m, is, ie, js, je, isd, ied, jsd, jed, nz, itt, ntr, do_any | |
| 114 | integer :: isv, iev, jsv, jev ! The valid range of the indices. | |
| 115 | integer :: IsdB, IedB, JsdB, JedB | |
| 116 | integer :: stencil_local ! Stencil for the local adection scheme | |
| 117 | 12 | integer :: local_advect_scheme(Reg%ntr) ! contains the list of the advection for each tracer |
| 118 | integer :: domore_k_tmp | |
| 119 | ||
| 120 | 12 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke |
| 121 | 12 | isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed |
| 122 | 12 | IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB |
| 123 | 12 | landvolfill = 1.0e-20 ! This is arbitrary, but must be positive. |
| 124 | 12 | stencil = 2 ! The scheme's stencil; 2 for PLM |
| 125 | ||
| 126 | 12 | ntr = Reg%ntr |
| 127 | 12 | Idt = 1.0 / dt |
| 128 | ||
| 129 | 12 | if (.not. associated(CS)) call MOM_error(FATAL, "MOM_tracer_advect: "// & |
| 130 | 0 | "tracer_advect_init must be called before advect_tracer.") |
| 131 | 12 | if (.not. associated(Reg)) call MOM_error(FATAL, "MOM_tracer_advect: "// & |
| 132 | 0 | "register_tracer must be called before advect_tracer.") |
| 133 | 12 | if (Reg%ntr==0) return |
| 134 | ||
| 135 | !$omp target enter data map(to: OBC, Reg, Reg%Tr(:)) map(alloc: domore_u, domore_v, uhr, vhr, uh_neglect, & | |
| 136 | !$omp vh_neglect, hprev, local_advect_scheme) | |
| 137 | ||
| 138 | 12 | do concurrent (k=1:nz, j=jsd:jed) |
| 139 | 62028 | domore_u(j,k) = .false. |
| 140 | enddo | |
| 141 | 12 | do concurrent (k=1:nz, j=jsdB:jedB) |
| 142 | 62940 | domore_v(j,k) = .false. |
| 143 | enddo | |
| 144 | ||
| 145 | 12 | call cpu_clock_begin(id_clock_advect) |
| 146 | 12 | x_first = (MOD(G%first_direction,2) == 0) |
| 147 | ||
| 148 | ! Choose the maximum stencil from all the local advection scheme | |
| 149 | 12 | do concurrent (m = 1:ntr) |
| 150 | ||
| 151 | 36 | local_advect_scheme(m) = Reg%Tr(m)%advect_scheme |
| 152 | 36 | if (local_advect_scheme(m) < 0) local_advect_scheme(m) = CS%default_advect_scheme |
| 153 | ||
| 154 | 36 | if (local_advect_scheme(m) == ADVECT_PLM) then |
| 155 | 0 | stencil_local = 2 |
| 156 | 36 | elseif (local_advect_scheme(m) == ADVECT_PPM) then |
| 157 | 0 | stencil_local = 3 |
| 158 | 36 | elseif (local_advect_scheme(m) == ADVECT_PPMH3) then |
| 159 | 36 | if (CS%useHuynhStencilBug) then |
| 160 | 0 | stencil_local = 2 |
| 161 | else | |
| 162 | 36 | stencil_local = 3 |
| 163 | endif | |
| 164 | endif | |
| 165 | 48 | stencil = max(stencil, stencil_local) |
| 166 | enddo | |
| 167 | ||
| 168 | !$omp target update from(local_advect_scheme) | |
| 169 | ||
| 170 | 12 | if (min(is-isd,ied-ie,js-jsd,jed-je) < stencil) then |
| 171 | call MOM_error(FATAL, "MOM_tracer_advect: "//& | |
| 172 | 0 | "stencil is wider than the halo.") |
| 173 | endif | |
| 174 | ||
| 175 | 12 | max_iter = 2*INT(CEILING(dt/CS%dt)) + 1 |
| 176 | ||
| 177 | 12 | if (present(max_iter_in)) max_iter = max_iter_in |
| 178 | 12 | if (present(x_first_in)) x_first = x_first_in |
| 179 | 12 | call cpu_clock_begin(id_clock_pass) |
| 180 | 12 | call create_group_pass(CS%pass_uhr_vhr_t_hprev, uhr, vhr, G%Domain) |
| 181 | 12 | call create_group_pass(CS%pass_uhr_vhr_t_hprev, hprev, G%Domain) |
| 182 | 48 | do m=1,ntr |
| 183 | 48 | call create_group_pass(CS%pass_uhr_vhr_t_hprev, Reg%Tr(m)%t, G%Domain) |
| 184 | enddo | |
| 185 | 12 | call cpu_clock_end(id_clock_pass) |
| 186 | ||
| 187 | ! This initializes the halos of uhr and vhr because pass_vector might do | |
| 188 | ! calculations on them, even though they are never used. | |
| 189 | 12 | do concurrent (k=1:nz) |
| 190 | 900 | do concurrent (j=jsd:jed, I=IsdB:IedB) |
| 191 | 8011800 | uhr(I,j,k) = 0.0 |
| 192 | enddo | |
| 193 | 900 | do concurrent (J=jsdB:jedB, i=Isd:Ied) |
| 194 | 8064900 | vhr(i,J,k) = 0.0 |
| 195 | enddo | |
| 196 | 900 | do concurrent (j=jsd:jed, i=Isd:Ied) |
| 197 | 7949700 | hprev(i,j,k) = 0.0 |
| 198 | enddo | |
| 199 | 900 | domore_k(k)=1 |
| 200 | ! Put the remaining (total) thickness fluxes into uhr and vhr. | |
| 201 | 900 | do concurrent (j=js:je, I=is-1:ie) |
| 202 | 6643800 | uhr(I,j,k) = uhtr(I,j,k) |
| 203 | enddo | |
| 204 | 900 | do concurrent (J=js-1:je, i=is:ie) |
| 205 | 6696900 | vhr(i,J,k) = vhtr(i,J,k) |
| 206 | enddo | |
| 207 | 1812 | if (.not. present(vol_prev)) then |
| 208 | ! This loop reconstructs the thickness field the last time that the | |
| 209 | ! tracers were updated, probably just after the diabatic forcing. A useful | |
| 210 | ! diagnostic could be to compare this reconstruction with that older value. | |
| 211 | 900 | do concurrent (j=js:je, i=is:ie) |
| 212 | hprev(i,j,k) = max(0.0, G%areaT(i,j)*h_end(i,j,k) + & | |
| 213 | 6480000 | ((uhr(I,j,k) - uhr(I-1,j,k)) + (vhr(i,J,k) - vhr(i,J-1,k)))) |
| 214 | ! In the case that the layer is now dramatically thinner than it was previously, | |
| 215 | ! add a bit of mass to avoid truncation errors. This will lead to | |
| 216 | ! non-conservation of tracers | |
| 217 | hprev(i,j,k) = hprev(i,j,k) + & | |
| 218 | 6588900 | max(0.0, 1.0e-13*hprev(i,j,k) - G%areaT(i,j)*h_end(i,j,k)) |
| 219 | enddo | |
| 220 | else | |
| 221 | 0 | do concurrent (j=js:je, i=is:ie) |
| 222 | 0 | hprev(i,j,k) = vol_prev(i,j,k) |
| 223 | enddo | |
| 224 | endif | |
| 225 | enddo | |
| 226 | ||
| 227 | 12 | do concurrent (j=jsd:jed, I=isd:ied-1) |
| 228 | 105168 | uh_neglect(I,j) = GV%H_subroundoff * MIN(G%areaT(i,j), G%areaT(i+1,j)) |
| 229 | enddo | |
| 230 | 1548 | do concurrent (J=jsd:jed-1, i=isd:ied) |
| 231 | 104460 | vh_neglect(i,J) = GV%H_subroundoff * MIN(G%areaT(i,j), G%areaT(i,j+1)) |
| 232 | enddo | |
| 233 | ||
| 234 | ! update GPU copy of Tr(:)%t | |
| 235 | ! only update t because other members are zeroed | |
| 236 | !$ do m=1,ntr | |
| 237 | !$omp target enter data map(to: Reg%Tr(m)%t) & | |
| 238 | !$omp map(alloc: Reg%Tr(m)%ad_x, Reg%Tr(m)%ad_y, Reg%Tr(m)%ad2d_x, Reg%Tr(m)%ad2d_y, & | |
| 239 | !$omp Reg%Tr(m)%advection_xy) | |
| 240 | !$ enddo | |
| 241 | ||
| 242 | ! initialize diagnostic fluxes and tendencies | |
| 243 | 12 | do concurrent (m=1:ntr) |
| 244 | 36 | if (associated(Reg%Tr(m)%ad_x)) then |
| 245 | 0 | do concurrent (k=1:nz, j=jsd:jed, I=IsdB:IedB) |
| 246 | 0 | Reg%Tr(m)%ad_x(i,j,k) = 0.0 |
| 247 | enddo | |
| 248 | endif | |
| 249 | 36 | if (associated(Reg%Tr(m)%ad_y)) then |
| 250 | 0 | do concurrent (k=1:nz, j=JsdB:jedB, i=isd:ied) |
| 251 | 0 | Reg%Tr(m)%ad_y(i,j,k) = 0.0 |
| 252 | enddo | |
| 253 | endif | |
| 254 | 36 | if (associated(Reg%Tr(m)%advection_xy)) then |
| 255 | 0 | do concurrent (k=1:nz, j=jsd:jed, i=isd:ied) |
| 256 | 0 | Reg%Tr(m)%advection_xy(i,j,k) = 0.0 |
| 257 | enddo | |
| 258 | endif | |
| 259 | 36 | if (associated(Reg%Tr(m)%ad2d_x)) then |
| 260 | 0 | do concurrent (j=jsd:jed, I=IsdB:IedB) |
| 261 | 0 | Reg%Tr(m)%ad2d_x(i,j) = 0.0 |
| 262 | enddo | |
| 263 | endif | |
| 264 | 48 | if (associated(Reg%Tr(m)%ad2d_y)) then |
| 265 | 0 | do concurrent (J=JsdB:JedB, i=isd:ied) |
| 266 | 0 | Reg%Tr(m)%ad2d_y(i,j) = 0.0 |
| 267 | enddo | |
| 268 | endif | |
| 269 | enddo | |
| 270 | ||
| 271 | 12 | isv = is ; iev = ie ; jsv = js ; jev = je |
| 272 | 12 | nsten_halo = min(is - isd, ied - ie, js - jsd, jed - je) / stencil |
| 273 | ||
| 274 | 19 | do itt=1,max_iter |
| 275 | ||
| 276 | 19 | if (isv > is-stencil) then |
| 277 | 19 | call do_group_pass(CS%pass_uhr_vhr_t_hprev, G%Domain, clock=id_clock_pass, omp_offload=.true.) |
| 278 | ||
| 279 | 19 | isv = is - nsten_halo * stencil ; jsv = js - nsten_halo * stencil |
| 280 | 19 | iev = ie + nsten_halo * stencil ; jev = je + nsten_halo * stencil |
| 281 | ! Reevaluate domore_u & domore_v unless the valid range is the same size as | |
| 282 | ! before. Also, do this if there is Strang splitting. | |
| 283 | 19 | if ((nsten_halo > 1) .or. (itt==1)) then |
| 284 | 912 | do concurrent (k=1:nz, domore_k(k) > 0) |
| 285 | 60300 | do concurrent (j=jsv:jev, .not.domore_u(j,k)) |
| 286 | 1220400 | do i=isv+stencil-1,iev-stencil ; if (uhr(I,j,k) /= 0.0) then |
| 287 | 54000 | domore_u(j,k) = .true. ; exit |
| 288 | endif ; enddo ! i-loop | |
| 289 | enddo | |
| 290 | 55800 | do concurrent (J=jsv+stencil-1:jev-stencil, .not.domore_v(J,k)) |
| 291 | 740700 | do i=isv+stencil,iev-stencil ; if (vhr(i,J,k) /= 0.0) then |
| 292 | 53100 | domore_v(J,k) = .true. ; exit |
| 293 | endif ; enddo ! i-loop | |
| 294 | enddo | |
| 295 | ||
| 296 | ! At this point, domore_k is global. Change it so that it indicates | |
| 297 | ! whether any work is needed on a layer on this processor. | |
| 298 | 900 | domore_k_tmp = 0 |
| 299 | 60300 | do concurrent (j=jsv:jev, domore_u(j,k)) DO_LOCALITY(reduce(max:domore_k_tmp)) |
| 300 | 60300 | domore_k_tmp = 1 |
| 301 | enddo | |
| 302 | 55800 | do concurrent (J=jsv+stencil-1:jev-stencil, domore_v(J,k)) DO_LOCALITY(reduce(max:domore_k_tmp)) |
| 303 | 55800 | domore_k_tmp = 1 |
| 304 | enddo | |
| 305 | 912 | domore_k(k) = domore_k_tmp |
| 306 | enddo ! k-loop | |
| 307 | endif | |
| 308 | endif | |
| 309 | ||
| 310 | ! Set the range of valid points after this iteration. | |
| 311 | 19 | isv = isv + stencil ; iev = iev - stencil |
| 312 | 19 | jsv = jsv + stencil ; jev = jev - stencil |
| 313 | ||
| 314 | ! To ensure positive definiteness of the thickness at each iteration, the | |
| 315 | ! mass fluxes out of each layer are checked each step, and limited to keep | |
| 316 | ! the thicknesses positive. This means that several iterations may be required | |
| 317 | ! for all the transport to happen. The sum over domore_k keeps the processors | |
| 318 | ! synchronized. This may not be very efficient, but it should be reliable. | |
| 319 | ||
| 320 | 19 | if (x_first) then |
| 321 | ||
| 322 | 1444 | do k=1,nz ; if (domore_k(k) > 0) then |
| 323 | ! First, advect zonally. | |
| 324 | call advect_x(Reg%Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & | |
| 325 | isv, iev, jsv-stencil, jev+stencil, k, G, GV, US, & | |
| 326 | 921 | local_advect_scheme) |
| 327 | endif ; enddo | |
| 328 | ||
| 329 | 1444 | do k=1,nz ; if (domore_k(k) > 0) then |
| 330 | ! Next, advect meridionally. | |
| 331 | call advect_y(Reg%Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & | |
| 332 | 921 | isv, iev, jsv, jev, k, G, GV, US, local_advect_scheme) |
| 333 | ||
| 334 | ! Update domore_k(k) for the next iteration | |
| 335 | 921 | domore_k_tmp = 0 |
| 336 | 61707 | do concurrent (j=jsv-stencil:jev+stencil, domore_u(j,k)) DO_LOCALITY(reduce(max:domore_k_tmp)) |
| 337 | 61707 | domore_k_tmp = 1 |
| 338 | enddo | |
| 339 | 57102 | do concurrent (J=jsv-1:jev, domore_v(J,k)) DO_LOCALITY(reduce(max:domore_k_tmp)) |
| 340 | 57102 | domore_k_tmp = 1 |
| 341 | enddo | |
| 342 | 921 | domore_k(k) = domore_k_tmp |
| 343 | endif ; enddo | |
| 344 | else | |
| 345 | 0 | do k=1,nz ; if (domore_k(k) > 0) then |
| 346 | ! First, advect meridionally. | |
| 347 | call advect_y(Reg%Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & | |
| 348 | isv-stencil, iev+stencil, jsv, jev, k, G, GV, US, & | |
| 349 | 0 | local_advect_scheme) |
| 350 | endif ; enddo | |
| 351 | ||
| 352 | 0 | do k=1,nz ; if (domore_k(k) > 0) then |
| 353 | ! Next, advect zonally. | |
| 354 | call advect_x(Reg%Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & | |
| 355 | 0 | isv, iev, jsv, jev, k, G, GV, US, local_advect_scheme) |
| 356 | ||
| 357 | ! Update domore_k(k) for the next iteration | |
| 358 | 0 | domore_k_tmp = 0 |
| 359 | 0 | do concurrent (j=jsv:jev, domore_u(j,k)) DO_LOCALITY(reduce(max:domore_k_tmp)) |
| 360 | 0 | domore_k_tmp = 1 |
| 361 | enddo | |
| 362 | 0 | do concurrent (J=jsv-1:jev, domore_v(J,k)) DO_LOCALITY(reduce(max:domore_k_tmp)) |
| 363 | 0 | domore_k_tmp = 1 |
| 364 | enddo | |
| 365 | 0 | domore_k(k) = domore_k_tmp |
| 366 | endif ; enddo | |
| 367 | ||
| 368 | endif ! x_first | |
| 369 | ||
| 370 | ! If the advection just isn't finishing after max_iter, move on. | |
| 371 | 19 | if (itt >= max_iter) then |
| 372 | 0 | exit |
| 373 | endif | |
| 374 | ||
| 375 | ! Exit if there are no layers that need more iterations. | |
| 376 | 19 | if (isv > is-stencil) then |
| 377 | 19 | do_any = 0 |
| 378 | 19 | call cpu_clock_begin(id_clock_sync) |
| 379 | 19 | call sum_across_PEs(domore_k(:), nz) |
| 380 | 19 | call cpu_clock_end(id_clock_sync) |
| 381 | 1444 | do k=1,nz ; do_any = do_any + domore_k(k) ; enddo |
| 382 | 19 | if (do_any == 0) then |
| 383 | 12 | exit |
| 384 | endif | |
| 385 | ||
| 386 | endif | |
| 387 | ||
| 388 | enddo ! Iterations loop | |
| 389 | ||
| 390 | !$ do m = 1, ntr | |
| 391 | !$omp target exit data map(from: Reg%Tr(m)%t, Reg%Tr(m)%ad_x, Reg%Tr(m)%ad_y, Reg%Tr(m)%ad2d_x, & | |
| 392 | !$omp Reg%Tr(m)%ad2d_y, Reg%Tr(m)%advection_xy) | |
| 393 | !$ enddo | |
| 394 | ||
| 395 | 12 | if (present(uhr_out)) then |
| 396 | 0 | do concurrent (k=1:nz, j=jsd:jed, i=isdB:iedB) |
| 397 | 0 | uhr_out(i,j,k) = uhr(i,j,k) |
| 398 | enddo | |
| 399 | endif | |
| 400 | 12 | if (present(vhr_out)) then |
| 401 | 0 | do concurrent (k=1:nz, j=jsdB:jedB, i=isd:ied) |
| 402 | 0 | vhr_out(i,j,k) = vhr(i,j,k) |
| 403 | enddo | |
| 404 | endif | |
| 405 | 12 | if (present(vol_prev) .and. present(update_vol_prev)) then |
| 406 | 0 | if (update_vol_prev) then |
| 407 | 0 | do concurrent (k=1:nz, j=jsd:jed, i=isd:ied) |
| 408 | 0 | vol_prev(i,j,k) = hprev(i,j,k) |
| 409 | enddo | |
| 410 | endif | |
| 411 | endif | |
| 412 | ||
| 413 | !$omp target exit data map(release: hprev, uhr, vhr, uh_neglect, vh_neglect, domore_u, & | |
| 414 | !$omp domore_v, local_advect_scheme, OBC, Reg, Reg%Tr(:)) | |
| 415 | ||
| 416 | 12 | call cpu_clock_end(id_clock_advect) |
| 417 | ||
| 418 | 12 | end subroutine advect_tracer |
| 419 | ||
| 420 | ||
| 421 | !> This subroutine does 1-d flux-form advection in the zonal direction using | |
| 422 | !! a monotonic piecewise linear scheme. | |
| 423 | 921 | subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & |
| 424 | 921 | is, ie, js, je, k, G, GV, US, advect_schemes) |
| 425 | type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure | |
| 426 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure | |
| 427 | integer, intent(in) :: ntr !< The number of tracers | |
| 428 | type(tracer_type), dimension(ntr), intent(inout) :: Tr !< The array of registered tracers to work on | |
| 429 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: hprev !< cell volume at the end of previous | |
| 430 | !! tracer change [H L2 ~> m3 or kg] | |
| 431 | real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: uhr !< accumulated volume/mass flux through | |
| 432 | !! the zonal face [H L2 ~> m3 or kg] | |
| 433 | real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: uh_neglect !< A tiny zonal mass flux that can | |
| 434 | !! be neglected [H L2 ~> m3 or kg] | |
| 435 | type(ocean_OBC_type), pointer :: OBC !< specifies whether, where, and what OBCs are used | |
| 436 | logical, dimension(SZJ_(G),SZK_(GV)), intent(inout) :: domore_u !< If true, there is more advection to be | |
| 437 | !! done in this u-row | |
| 438 | real, intent(in) :: Idt !< The inverse of dt [T-1 ~> s-1] | |
| 439 | integer, intent(in) :: is !< The starting tracer i-index to work on | |
| 440 | integer, intent(in) :: ie !< The ending tracer i-index to work on | |
| 441 | integer, intent(in) :: js !< The starting tracer j-index to work on | |
| 442 | integer, intent(in) :: je !< The ending tracer j-index to work on | |
| 443 | integer, intent(in) :: k !< The k-level to work on | |
| 444 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 445 | integer, dimension(ntr), intent(in) :: advect_schemes !< list of advection schemes to use | |
| 446 | ||
| 447 | real, dimension(SZI_(G),ntr) :: & | |
| 448 | 1842 | slope_x ! The concentration slope per grid point [conc]. |
| 449 | real, dimension(SZIB_(G),SZJ_(G),ntr) :: & | |
| 450 | 1842 | flux_x ! The tracer flux across a boundary [H L2 conc ~> m3 conc or kg conc]. |
| 451 | real, dimension(SZI_(G),ntr) :: & | |
| 452 | 1842 | T_tmp ! The copy of the tracer concentration at constant i,k [conc]. |
| 453 | ||
| 454 | real :: hup, hlos ! hup is the upwind volume, hlos is the | |
| 455 | ! part of that volume that might be lost | |
| 456 | ! due to advection out the other side of | |
| 457 | ! the grid box, both in [H L2 ~> m3 or kg]. | |
| 458 | 1842 | real :: uhh(SZIB_(G)) ! The zonal flux that occurs during the |
| 459 | ! current iteration [H L2 ~> m3 or kg]. | |
| 460 | real, dimension(SZIB_(G)) :: & | |
| 461 | 1842 | hlst, & ! Work variable [H L2 ~> m3 or kg]. |
| 462 | 1842 | Ihnew, & ! Work variable [H-1 L-2 ~> m-3 or kg-1]. |
| 463 | 1842 | CFL ! The absolute value of the advective upwind-cell CFL number [nondim]. |
| 464 | real :: min_h ! The minimum thickness that can be realized during | |
| 465 | ! any of the passes [H ~> m or kg m-2]. | |
| 466 | real :: tiny_h ! The smallest numerically invertible thickness [H ~> m or kg m-2]. | |
| 467 | real :: h_neglect ! A thickness that is so small it is usually lost | |
| 468 | ! in roundoff and can be neglected [H ~> m or kg m-2]. | |
| 469 | real :: aR, aL ! Reconstructed tracer concentrations at the right and left edges [conc] | |
| 470 | real :: dMx ! Difference between the maximum of the surrounding cell concentrations and | |
| 471 | ! the value in the cell whose reconstruction is being found [conc] | |
| 472 | real :: dMn ! Difference between the tracer concentration in the cell whose reconstruction | |
| 473 | ! is being found and the minimum of the surrounding values [conc] | |
| 474 | real :: Tp, Tc, Tm ! Tracer concentrations around the upstream cell [conc] | |
| 475 | real :: dA ! Difference between the reconstruction tracer edge values [conc] | |
| 476 | real :: mA ! Average of the reconstruction tracer edge values [conc] | |
| 477 | real :: a6 ! Curvature of the reconstruction tracer values [conc] | |
| 478 | 1842 | logical :: do_i(SZI_(G),SZJ_(G)) ! If true, work on given points. |
| 479 | logical :: usePLMslope, domore_u_jk | |
| 480 | integer :: i, j, m, n, i_up, stencil, ntr_id | |
| 481 | type(OBC_segment_type), pointer :: segment=>NULL() | |
| 482 | 1842 | logical, dimension(SZJ_(G),SZK_(GV)) :: domore_u_initial |
| 483 | ||
| 484 | ! keep a local copy of the initial values of domore_u, which is to be used when computing ad2d_x | |
| 485 | ! diagnostic at the end of this subroutine. | |
| 486 | 4767096 | domore_u_initial = domore_u |
| 487 | ||
| 488 | 921 | usePLMslope = .false. |
| 489 | ! stencil for calculating slope values | |
| 490 | 921 | stencil = 1 |
| 491 | 3684 | do m = 1,ntr |
| 492 | 2763 | if ((advect_schemes(m) == ADVECT_PLM) .or. (advect_schemes(m) == ADVECT_PPM)) & |
| 493 | 0 | usePLMslope = .true. |
| 494 | 3684 | if (advect_schemes(m) == ADVECT_PPM) stencil = 2 |
| 495 | enddo | |
| 496 | ||
| 497 | 921 | min_h = 0.1*GV%Angstrom_H |
| 498 | 921 | tiny_h = tiny(min_h) |
| 499 | 921 | h_neglect = GV%H_subroundoff |
| 500 | ||
| 501 | ! do I=is-1,ie ; CFL(I) = 0.0 ; enddo | |
| 502 | ||
| 503 | !$omp target enter data & | |
| 504 | !$omp map(alloc: slope_x, T_tmp, uhh, CFL, hlst, Ihnew, do_i, flux_x) | |
| 505 | !$omp target teams loop private(slope_x, T_tmp, uhh, CFL, hlst, Ihnew, Tp, dMx, dMn, m, i, n, & | |
| 506 | !$omp ntr_id, hup, hlos, i_up, Tc, Tm, aL, aR, dA, mA, a6, domore_u_jk) | |
| 507 | 61707 | do j=js,je ; if (domore_u(j,k)) then |
| 508 | 54098 | domore_u_jk = .false. |
| 509 | ||
| 510 | ! Calculate the i-direction profiles (slopes) of each tracer that is being advected. | |
| 511 | 54098 | if (usePLMslope) then |
| 512 | 0 | do concurrent (m=1:ntr, i=is-stencil:ie+stencil) |
| 513 | !if (ABS(Tr(m)%t(i+1,j,k)-Tr(m)%t(i,j,k)) < & | |
| 514 | ! ABS(Tr(m)%t(i,j,k)-Tr(m)%t(i-1,j,k))) then | |
| 515 | ! maxslope = 4.0*(Tr(m)%t(i+1,j,k)-Tr(m)%t(i,j,k)) | |
| 516 | !else | |
| 517 | ! maxslope = 4.0*(Tr(m)%t(i,j,k)-Tr(m)%t(i-1,j,k)) | |
| 518 | !endif | |
| 519 | !if ((Tr(m)%t(i+1,j,k)-Tr(m)%t(i,j,k)) * (Tr(m)%t(i,j,k)-Tr(m)%t(i-1,j,k)) < 0.0) then | |
| 520 | ! slope_x(i,m) = 0.0 | |
| 521 | !elseif (ABS(Tr(m)%t(i+1,j,k)-Tr(m)%t(i-1,j,k))<ABS(maxslope)) then | |
| 522 | ! slope_x(i,m) = G%mask2dCu(I,j)*G%mask2dCu(I-1,j) * & | |
| 523 | ! 0.5*(Tr(m)%t(i+1,j,k)-Tr(m)%t(i-1,j,k)) | |
| 524 | !else | |
| 525 | ! slope_x(i,m) = G%mask2dCu(I,j)*G%mask2dCu(I-1,j) * 0.5*maxslope | |
| 526 | !endif | |
| 527 | 0 | Tp = Tr(m)%t(i+1,j,k) ; Tc = Tr(m)%t(i,j,k) ; Tm = Tr(m)%t(i-1,j,k) |
| 528 | 0 | dMx = max( Tp, Tc, Tm ) - Tc |
| 529 | 0 | dMn= Tc - min( Tp, Tc, Tm ) |
| 530 | slope_x(i,m) = G%mask2dCu(I,j)*G%mask2dCu(I-1,j) * & | |
| 531 | 0 | sign( min(0.5*abs(Tp-Tm), 2.0*dMx, 2.0*dMn), Tp-Tm ) |
| 532 | enddo | |
| 533 | endif ! usePLMslope | |
| 534 | ||
| 535 | ! make a copy of the tracers in case values need to be overridden for OBCs | |
| 536 | 54098 | do concurrent (m = 1:ntr, i=G%isd:G%ied) |
| 537 | 27752274 | T_tmp(i,m) = Tr(m)%t(i,j,k) |
| 538 | enddo | |
| 539 | ! loop through open boundaries and recalculate flux terms | |
| 540 | 54098 | if (associated(OBC)) then ; if (OBC%OBC_pe) then |
| 541 | 0 | do n=1,OBC%number_of_segments |
| 542 | ! segment=>OBC%segment(n) | |
| 543 | 0 | if (.not. associated(OBC%segment(n)%tr_Reg)) cycle |
| 544 | 0 | if (OBC%segment(n)%is_E_or_W) then |
| 545 | 0 | if (j>=OBC%segment(n)%HI%jsd .and. j<=OBC%segment(n)%HI%jed) then |
| 546 | 0 | I = OBC%segment(n)%HI%IsdB |
| 547 | 0 | do concurrent (m = 1:OBC%segment(n)%tr_Reg%ntseg) ! replace tracers with OBC values |
| 548 | 0 | ntr_id = OBC%segment(n)%tr_reg%Tr(m)%ntr_index |
| 549 | 0 | if (allocated(OBC%segment(n)%tr_Reg%Tr(m)%tres)) then |
| 550 | 0 | if (OBC%segment(n)%direction == OBC_DIRECTION_W) then |
| 551 | 0 | T_tmp(i,ntr_id) = OBC%segment(n)%tr_Reg%Tr(m)%tres(i,j,k) |
| 552 | else | |
| 553 | 0 | T_tmp(i+1,ntr_id) = OBC%segment(n)%tr_Reg%Tr(m)%tres(i,j,k) |
| 554 | endif | |
| 555 | else | |
| 556 | 0 | if (OBC%segment(n)%direction == OBC_DIRECTION_W) then |
| 557 | 0 | T_tmp(i,ntr_id) = OBC%segment(n)%tr_Reg%Tr(m)%OBC_inflow_conc |
| 558 | else | |
| 559 | 0 | T_tmp(i+1,ntr_id) = OBC%segment(n)%tr_Reg%Tr(m)%OBC_inflow_conc |
| 560 | endif | |
| 561 | endif | |
| 562 | enddo | |
| 563 | ! Apply update tracer values for slope calculation | |
| 564 | 0 | do concurrent (m = 1:ntr, i=OBC%segment(n)%HI%IsdB-1:OBC%segment(n)%HI%IsdB+1) |
| 565 | 0 | Tp = T_tmp(i+1,m) ; Tc = T_tmp(i,m) ; Tm = T_tmp(i-1,m) |
| 566 | 0 | dMx = max( Tp, Tc, Tm ) - Tc |
| 567 | 0 | dMn= Tc - min( Tp, Tc, Tm ) |
| 568 | slope_x(i,m) = G%mask2dCu(I,j)*G%mask2dCu(I-1,j) * & | |
| 569 | 0 | sign( min(0.5*abs(Tp-Tm), 2.0*dMx, 2.0*dMn), Tp-Tm ) |
| 570 | enddo | |
| 571 | ||
| 572 | endif | |
| 573 | endif | |
| 574 | enddo | |
| 575 | endif ; endif | |
| 576 | ||
| 577 | ||
| 578 | ! Calculate the i-direction fluxes of each tracer, using as much | |
| 579 | ! the minimum of the remaining mass flux (uhr) and the half the mass | |
| 580 | ! in the cell plus whatever part of its half of the mass flux that | |
| 581 | ! the flux through the other side does not require. | |
| 582 | 54098 | do concurrent (I=is-1:ie) DO_LOCALITY(reduce(.or.:domore_u_jk)) |
| 583 | if ((uhr(I,j,k) == 0.0) .or. & | |
| 584 | 6545858 | ((uhr(I,j,k) < 0.0) .and. (hprev(i+1,j,k) <= tiny_h)) .or. & |
| 585 | 6599956 | ((uhr(I,j,k) > 0.0) .and. (hprev(i,j,k) <= tiny_h)) ) then |
| 586 | 2092551 | uhh(I) = 0.0 |
| 587 | 2092551 | CFL(I) = 0.0 |
| 588 | 4453307 | elseif (uhr(I,j,k) < 0.0) then |
| 589 | 2634018 | hup = hprev(i+1,j,k) - G%areaT(i+1,j)*min_h |
| 590 | 2634018 | hlos = MAX(0.0, uhr(I+1,j,k)) |
| 591 | 2634018 | if ((((hup - hlos) + uhr(I,j,k)) < 0.0) .and. & |
| 592 | ((0.5*hup + uhr(I,j,k)) < 0.0)) then | |
| 593 | 0 | uhh(I) = MIN(-0.5*hup, -hup+hlos, 0.0) |
| 594 | 0 | domore_u_jk = .true. |
| 595 | else | |
| 596 | 2634018 | uhh(I) = uhr(I,j,k) |
| 597 | endif | |
| 598 | 2634018 | CFL(I) = - uhh(I) / (hprev(i+1,j,k)) ! CFL is positive |
| 599 | else | |
| 600 | 1819289 | hup = hprev(i,j,k) - G%areaT(i,j)*min_h |
| 601 | 1819289 | hlos = MAX(0.0, -uhr(I-1,j,k)) |
| 602 | 1819289 | if ((((hup - hlos) - uhr(I,j,k)) < 0.0) .and. & |
| 603 | ((0.5*hup - uhr(I,j,k)) < 0.0)) then | |
| 604 | 107 | uhh(I) = MAX(0.5*hup, hup-hlos, 0.0) |
| 605 | 107 | domore_u_jk = .true. |
| 606 | else | |
| 607 | 1819182 | uhh(I) = uhr(I,j,k) |
| 608 | endif | |
| 609 | 1819289 | CFL(I) = uhh(I) / (hprev(i,j,k)) ! CFL is positive |
| 610 | endif | |
| 611 | enddo | |
| 612 | ||
| 613 | 54098 | domore_u(j,k) = domore_u_jk |
| 614 | ||
| 615 | 54098 | do concurrent (m=1:ntr) |
| 616 | ||
| 617 | 216392 | if ((advect_schemes(m) == ADVECT_PPM) .or. (advect_schemes(m) == ADVECT_PPMH3)) then |
| 618 | 162294 | do concurrent (I=is-1:ie) |
| 619 | ! centre cell depending on upstream direction | |
| 620 | 19637574 | if (uhh(I) >= 0.0) then |
| 621 | 11735520 | i_up = i |
| 622 | else | |
| 623 | 7902054 | i_up = i+1 |
| 624 | endif | |
| 625 | ||
| 626 | ! Implementation of PPM-H3 | |
| 627 | 19637574 | Tp = T_tmp(i_up+1,m) ; Tc = T_tmp(i_up,m) ; Tm = T_tmp(i_up-1,m) |
| 628 | ||
| 629 | 19637574 | if (advect_schemes(m) == ADVECT_PPMH3) then |
| 630 | 19637574 | aL = ( 5.*Tc + ( 2.*Tm - Tp ) )/6. ! H3 estimate |
| 631 | 19637574 | aL = max( min(Tc,Tm), aL) ; aL = min( max(Tc,Tm), aL) ! Bound |
| 632 | 19637574 | aR = ( 5.*Tc + ( 2.*Tp - Tm ) )/6. ! H3 estimate |
| 633 | 19637574 | aR = max( min(Tc,Tp), aR) ; aR = min( max(Tc,Tp), aR) ! Bound |
| 634 | else | |
| 635 | 0 | aL = 0.5 * ((Tm + Tc) + (slope_x(i_up-1,m) - slope_x(i_up,m)) / 3.) |
| 636 | 0 | aR = 0.5 * ((Tc + Tp) + (slope_x(i_up,m) - slope_x(i_up+1,m)) / 3.) |
| 637 | endif | |
| 638 | ||
| 639 | 19637574 | dA = aR - aL ; mA = 0.5*( aR + aL ) |
| 640 | 19637574 | if (G%mask2dCu(I_up,j)*G%mask2dCu(I_up-1,j)*(Tp-Tc)*(Tc-Tm) <= 0.) then |
| 641 | 12890671 | aL = Tc ; aR = Tc ! PCM for local extrema and boundary cells |
| 642 | 6746903 | elseif ( dA*(Tc-mA) > (dA*dA)/6. ) then |
| 643 | 652080 | aL = (3.*Tc) - 2.*aR |
| 644 | 6094823 | elseif ( dA*(Tc-mA) < - (dA*dA)/6. ) then |
| 645 | 675197 | aR = (3.*Tc) - 2.*aL |
| 646 | endif | |
| 647 | ||
| 648 | 19637574 | a6 = 6.*Tc - 3. * (aR + aL) ! Curvature |
| 649 | ||
| 650 | 39437442 | if (uhh(I) >= 0.0) then |
| 651 | flux_x(I,j,m) = uhh(I)*( aR - 0.5 * CFL(I) * ( & | |
| 652 | 11735520 | ( aR - aL ) - a6 * ( 1. - 2./3. * CFL(I) ) ) ) |
| 653 | else | |
| 654 | flux_x(I,j,m) = uhh(I)*( aL + 0.5 * CFL(I) * ( & | |
| 655 | 7902054 | ( aR - aL ) + a6 * ( 1. - 2./3. * CFL(I) ) ) ) |
| 656 | endif | |
| 657 | enddo | |
| 658 | else ! PLM | |
| 659 | 0 | do concurrent (I=is-1:ie) |
| 660 | 0 | if (uhh(I) >= 0.0) then |
| 661 | ! Indirect implementation of PLM | |
| 662 | !aL = Tr(m)%t(i,j,k) - 0.5 * slope_x(i,m) | |
| 663 | !aR = Tr(m)%t(i,j,k) + 0.5 * slope_x(i,m) | |
| 664 | !flux_x(I,j,m) = uhh(I)*( aR - 0.5 * (aR-aL) * CFL(I) ) | |
| 665 | ! Alternative implementation of PLM | |
| 666 | 0 | Tc = T_tmp(i,m) |
| 667 | 0 | flux_x(I,j,m) = uhh(I)*( Tc + 0.5 * slope_x(i,m) * ( 1. - CFL(I) ) ) |
| 668 | else | |
| 669 | ! Indirect implementation of PLM | |
| 670 | !aL = Tr(m)%t(i+1,j,k) - 0.5 * slope_x(i+1,m) | |
| 671 | !aR = Tr(m)%t(i+1,j,k) + 0.5 * slope_x(i+1,m) | |
| 672 | !flux_x(I,j,m) = uhh(I)*( aL + 0.5 * (aR-aL) * CFL(I) ) | |
| 673 | ! Alternative implementation of PLM | |
| 674 | 0 | Tc = T_tmp(i+1,m) |
| 675 | 0 | flux_x(I,j,m) = uhh(I)*( Tc - 0.5 * slope_x(i+1,m) * ( 1. - CFL(I) ) ) |
| 676 | endif | |
| 677 | enddo | |
| 678 | endif ! usePPM | |
| 679 | enddo | |
| 680 | ||
| 681 | 54098 | if (associated(OBC)) then ; if (OBC%OBC_pe) then |
| 682 | 0 | if (OBC%specified_u_BCs_exist_globally .or. OBC%open_u_BCs_exist_globally) then |
| 683 | 0 | do n=1,OBC%number_of_segments |
| 684 | ! segment=>OBC%segment(n) | |
| 685 | 0 | if (.not. associated(OBC%segment(n)%tr_Reg)) cycle |
| 686 | 0 | if (OBC%segment(n)%is_E_or_W) then |
| 687 | 0 | if (j>=OBC%segment(n)%HI%jsd .and. j<=OBC%segment(n)%HI%jed) then |
| 688 | 0 | I = OBC%segment(n)%HI%IsdB |
| 689 | ! Tracer fluxes are set to prescribed values only for inflows from masked areas. | |
| 690 | ! Now changing to simply fixed inflows. | |
| 691 | 0 | if ((uhr(I,j,k) > 0.0) .and. (OBC%segment(n)%direction == OBC_DIRECTION_W) .or. & |
| 692 | (uhr(I,j,k) < 0.0) .and. (OBC%segment(n)%direction == OBC_DIRECTION_E)) then | |
| 693 | 0 | uhh(I) = uhr(I,j,k) |
| 694 | ! should the reservoir evolve for this case Kate ?? - Nope | |
| 695 | 0 | do concurrent (m=1:OBC%segment(n)%tr_Reg%ntseg) |
| 696 | 0 | ntr_id = OBC%segment(n)%tr_reg%Tr(m)%ntr_index |
| 697 | 0 | if (allocated(OBC%segment(n)%tr_Reg%Tr(m)%tres)) then |
| 698 | 0 | flux_x(I,j,ntr_id) = uhh(I)*OBC%segment(n)%tr_Reg%Tr(m)%tres(I,j,k) |
| 699 | 0 | else ; flux_x(I,j,ntr_id) = uhh(I)*OBC%segment(n)%tr_Reg%Tr(m)%OBC_inflow_conc ; endif |
| 700 | enddo | |
| 701 | endif | |
| 702 | endif | |
| 703 | endif | |
| 704 | enddo | |
| 705 | endif | |
| 706 | ||
| 707 | 0 | if (OBC%open_u_BCs_exist_globally) then |
| 708 | 0 | do n=1,OBC%number_of_segments |
| 709 | ! segment=>OBC%segment(n) | |
| 710 | 0 | I = OBC%segment(n)%HI%IsdB |
| 711 | 0 | if (OBC%segment(n)%is_E_or_W .and. (j >= OBC%segment(n)%HI%jsd .and. j<= OBC%segment(n)%HI%jed)) then |
| 712 | 0 | if (OBC%segment(n)%specified) cycle |
| 713 | 0 | if (.not. associated(OBC%segment(n)%tr_Reg)) cycle |
| 714 | ||
| 715 | ! Tracer fluxes are set to prescribed values only for inflows from masked areas. | |
| 716 | 0 | if ((uhr(I,j,k) > 0.0) .and. (G%mask2dT(i,j) < 0.5) .or. & |
| 717 | (uhr(I,j,k) < 0.0) .and. (G%mask2dT(i+1,j) < 0.5)) then | |
| 718 | 0 | uhh(I) = uhr(I,j,k) |
| 719 | 0 | do concurrent (m=1:OBC%segment(n)%tr_Reg%ntseg) |
| 720 | 0 | ntr_id = OBC%segment(n)%tr_reg%Tr(m)%ntr_index |
| 721 | 0 | if (allocated(OBC%segment(n)%tr_Reg%Tr(m)%tres)) then |
| 722 | 0 | flux_x(I,j,ntr_id) = uhh(I)*OBC%segment(n)%tr_Reg%Tr(m)%tres(I,j,k) |
| 723 | 0 | else; flux_x(I,j,ntr_id) = uhh(I)*OBC%segment(n)%tr_Reg%Tr(m)%OBC_inflow_conc; endif |
| 724 | enddo | |
| 725 | endif | |
| 726 | endif | |
| 727 | enddo | |
| 728 | endif | |
| 729 | endif ; endif | |
| 730 | ||
| 731 | ! Calculate new tracer concentration in each cell after accounting | |
| 732 | ! for the i-direction fluxes. | |
| 733 | 54098 | do concurrent (I=is-1:ie) |
| 734 | 6545858 | uhr(I,j,k) = uhr(I,j,k) - uhh(I) |
| 735 | 6599956 | if (abs(uhr(I,j,k)) < uh_neglect(I,j)) uhr(I,j,k) = 0.0 |
| 736 | enddo | |
| 737 | 54098 | do concurrent (i=is:ie) |
| 738 | 6545858 | if ((uhh(I) /= 0.0) .or. (uhh(I-1) /= 0.0)) then |
| 739 | 4512814 | do_i(i,j) = .true. |
| 740 | 4512814 | hlst(i) = hprev(i,j,k) |
| 741 | 4512814 | hprev(i,j,k) = hprev(i,j,k) - (uhh(I) - uhh(I-1)) |
| 742 | 4512814 | if (hprev(i,j,k) <= 0.0) then ; do_i(i,j) = .false. |
| 743 | 4512814 | elseif (hprev(i,j,k) < h_neglect*G%areaT(i,j)) then |
| 744 | 0 | hlst(i) = hlst(i) + (h_neglect*G%areaT(i,j) - hprev(i,j,k)) |
| 745 | 0 | Ihnew(i) = 1.0 / (h_neglect*G%areaT(i,j)) |
| 746 | 4512814 | else ; Ihnew(i) = 1.0 / hprev(i,j,k) ; endif |
| 747 | else | |
| 748 | 1978946 | do_i(i,j) = .false. |
| 749 | endif | |
| 750 | enddo | |
| 751 | ||
| 752 | ! Update do_i so that nothing changes outside of the OBC (problem for interior OBCs only) | |
| 753 | 54098 | if (associated(OBC)) then |
| 754 | if (.not.OBC%exterior_OBC_bug .and. OBC%OBC_pe & | |
| 755 | 0 | .and. (OBC%specified_u_BCs_exist_globally .or. OBC%open_u_BCs_exist_globally)) then |
| 756 | ! OBC_DIRECTION_E / OBC_DIRECTION_W on the west / east edge | |
| 757 | 0 | do concurrent (i=is:ie, OBC%segnum_u(I-1,j) > 0 .or. OBC%segnum_u(I,j) < 0) |
| 758 | 0 | do_i(i,j) = .false. |
| 759 | enddo | |
| 760 | endif | |
| 761 | endif | |
| 762 | ||
| 763 | ! update tracer concentration from i-flux and save some diagnostics | |
| 764 | 54098 | do concurrent (m=1:ntr) |
| 765 | ||
| 766 | ! update tracer | |
| 767 | 162294 | do concurrent (i=is:ie) |
| 768 | 19637574 | if (do_i(i,j)) then |
| 769 | 13538442 | if (Ihnew(i) > 0.0) then |
| 770 | Tr(m)%t(i,j,k) = (Tr(m)%t(i,j,k) * hlst(i) - & | |
| 771 | 13538442 | (flux_x(I,j,m) - flux_x(I-1,j,m))) * Ihnew(i) |
| 772 | endif | |
| 773 | endif | |
| 774 | enddo | |
| 775 | ||
| 776 | ! diagnostics | |
| 777 | 162294 | if (associated(Tr(m)%ad_x)) then |
| 778 | 0 | do concurrent (I=is-1:ie, do_i(i,j) .or. do_i(i+1,j)) |
| 779 | 0 | Tr(m)%ad_x(I,j,k) = Tr(m)%ad_x(I,j,k) + flux_x(I,j,m)*Idt |
| 780 | enddo | |
| 781 | endif | |
| 782 | ||
| 783 | ! diagnose convergence of flux_x (do not use the Ihnew(i) part of the logic). | |
| 784 | ! division by areaT to get into W/m2 for heat and kg/(s*m2) for salt. | |
| 785 | 216392 | if (associated(Tr(m)%advection_xy)) then |
| 786 | 0 | do concurrent (i=is:ie, do_i(i,j)) |
| 787 | Tr(m)%advection_xy(i,j,k) = Tr(m)%advection_xy(i,j,k) - & | |
| 788 | (flux_x(I,j,m) - flux_x(I-1,j,m)) * & | |
| 789 | 0 | Idt * G%IareaT(i,j) |
| 790 | enddo | |
| 791 | endif | |
| 792 | ||
| 793 | enddo | |
| 794 | ||
| 795 | endif ; enddo ! End of j-loop. | |
| 796 | ||
| 797 | ! Do user controlled underflow of the tracer concentrations. | |
| 798 | 3684 | do concurrent (m=1:ntr, Tr(m)%conc_underflow > 0.0) |
| 799 | 3684 | do concurrent (j=js:je, i=is:ie) |
| 800 | 0 | if (abs(Tr(m)%t(i,j,k)) < Tr(m)%conc_underflow) Tr(m)%t(i,j,k) = 0.0 |
| 801 | enddo | |
| 802 | enddo | |
| 803 | ||
| 804 | ! compute ad2d_x diagnostic outside above j-loop so as to make the summation ordered when OMP is active. | |
| 805 | ||
| 806 | 3684 | do m=1,ntr ; if (associated(Tr(m)%ad2d_x)) then |
| 807 | 0 | do concurrent (j=js:je, domore_u_initial(j,k)) |
| 808 | 0 | do concurrent (I=is-1:ie, do_i(i,j) .or. do_i(i+1,j)) |
| 809 | 0 | Tr(m)%ad2d_x(I,j) = Tr(m)%ad2d_x(I,j) + flux_x(I,j,m)*Idt |
| 810 | enddo | |
| 811 | enddo | |
| 812 | endif ; enddo ! End of m-loop. | |
| 813 | ||
| 814 | !$omp target exit data & | |
| 815 | !$omp map(release: slope_x, T_tmp, uhh, CFL, hlst, Ihnew, do_i, flux_x) | |
| 816 | 921 | end subroutine advect_x |
| 817 | ||
| 818 | !> This subroutine does 1-d flux-form advection using a monotonic piecewise | |
| 819 | !! linear scheme. | |
| 820 | 921 | subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & |
| 821 | 921 | is, ie, js, je, k, G, GV, US, advect_schemes) |
| 822 | type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure | |
| 823 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure | |
| 824 | integer, intent(in) :: ntr !< The number of tracers | |
| 825 | type(tracer_type), dimension(ntr), intent(inout) :: Tr !< The array of registered tracers to work on | |
| 826 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: hprev !< cell volume at the end of previous | |
| 827 | !! tracer change [H L2 ~> m3 or kg] | |
| 828 | real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: vhr !< accumulated volume/mass flux through | |
| 829 | !! the meridional face [H L2 ~> m3 or kg] | |
| 830 | real, dimension(SZI_(G),SZJB_(G)), intent(inout) :: vh_neglect !< A tiny meridional mass flux that can | |
| 831 | !! be neglected [H L2 ~> m3 or kg] | |
| 832 | type(ocean_OBC_type), pointer :: OBC !< specifies whether, where, and what OBCs are used | |
| 833 | logical, dimension(SZJB_(G),SZK_(GV)), intent(inout) :: domore_v !< If true, there is more advection to be | |
| 834 | !! done in this v-row | |
| 835 | real, intent(in) :: Idt !< The inverse of dt [T-1 ~> s-1] | |
| 836 | integer, intent(in) :: is !< The starting tracer i-index to work on | |
| 837 | integer, intent(in) :: ie !< The ending tracer i-index to work on | |
| 838 | integer, intent(in) :: js !< The starting tracer j-index to work on | |
| 839 | integer, intent(in) :: je !< The ending tracer j-index to work on | |
| 840 | integer, intent(in) :: k !< The k-level to work on | |
| 841 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 842 | integer, dimension(ntr), intent(in) :: advect_schemes !< list of advection schemes to use | |
| 843 | ||
| 844 | real, dimension(SZI_(G),ntr,SZJ_(G)) :: & | |
| 845 | 1842 | slope_y ! The concentration slope per grid point [conc]. |
| 846 | real, dimension(SZI_(G),ntr,SZJB_(G)) :: & | |
| 847 | 1842 | flux_y ! The tracer flux across a boundary [H L2 conc ~> m3 conc or kg conc]. |
| 848 | real, dimension(SZI_(G),ntr,SZJB_(G)) :: & | |
| 849 | 1842 | T_tmp ! The copy of the tracer concentration at constant i,k [conc]. |
| 850 | 1842 | real :: vhh(SZI_(G),SZJB_(G)) ! The meridional flux that occurs during the |
| 851 | ! current iteration [H L2 ~> m3 or kg]. | |
| 852 | real :: hup, hlos ! hup is the upwind volume, hlos is the | |
| 853 | ! part of that volume that might be lost | |
| 854 | ! due to advection out the other side of | |
| 855 | ! the grid box, both in [H L2 ~> m3 or kg]. | |
| 856 | real, dimension(SZIB_(G)) :: & | |
| 857 | 1842 | hlst, & ! Work variable [H L2 ~> m3 or kg]. |
| 858 | 1842 | Ihnew, & ! Work variable [H-1 L-2 ~> m-3 or kg-1]. |
| 859 | 1842 | CFL ! The absolute value of the advective upwind-cell CFL number [nondim]. |
| 860 | real :: min_h ! The minimum thickness that can be realized during | |
| 861 | ! any of the passes [H ~> m or kg m-2]. | |
| 862 | real :: tiny_h ! The smallest numerically invertible thickness [H ~> m or kg m-2]. | |
| 863 | real :: h_neglect ! A thickness that is so small it is usually lost | |
| 864 | ! in roundoff and can be neglected [H ~> m or kg m-2]. | |
| 865 | real :: aR, aL ! Reconstructed tracer concentrations at the right and left edges [conc] | |
| 866 | real :: dMx ! Difference between the maximum of the surrounding cell concentrations and | |
| 867 | ! the value in the cell whose reconstruction is being found [conc] | |
| 868 | real :: dMn ! Difference between the tracer average in the cell whose reconstruction | |
| 869 | ! is being found and the minimum of the surrounding values [conc] | |
| 870 | real :: Tp, Tc, Tm ! Tracer concentrations around the upstream cell [conc] | |
| 871 | real :: dA ! Difference between the reconstruction tracer edge values [conc] | |
| 872 | real :: mA ! Average of the reconstruction tracer edge values [conc] | |
| 873 | real :: a6 ! Curvature of the reconstruction tracer values [conc] | |
| 874 | 1842 | logical :: do_j_tr(SZJ_(G)) ! If true, calculate the tracer profiles. |
| 875 | 1842 | logical :: do_i(SZI_(G), SZJ_(G)) ! If true, work on given points. |
| 876 | logical :: usePLMslope | |
| 877 | integer :: i, j, j2, m, n, j_up, stencil, ntr_id | |
| 878 | type(OBC_segment_type), pointer :: segment=>NULL() | |
| 879 | 1842 | logical :: domore_v_initial(SZJB_(G)), domore_v_jk ! Initial state of domore_v |
| 880 | ||
| 881 | 921 | usePLMslope = .false. |
| 882 | ! stencil for calculating slope values | |
| 883 | 921 | stencil = 1 |
| 884 | 3684 | do m = 1,ntr |
| 885 | 2763 | if ((advect_schemes(m) == ADVECT_PLM) .or. (advect_schemes(m) == ADVECT_PPM)) & |
| 886 | 0 | usePLMslope = .true. |
| 887 | 3684 | if (advect_schemes(m) == ADVECT_PPM) stencil = 2 |
| 888 | enddo | |
| 889 | ||
| 890 | 921 | min_h = 0.1*GV%Angstrom_H |
| 891 | 921 | tiny_h = tiny(min_h) |
| 892 | 921 | h_neglect = GV%H_subroundoff |
| 893 | ||
| 894 | !$omp target enter data map(alloc: vhh, T_tmp, slope_y, flux_y, domore_v_initial, do_j_tr, do_i) | |
| 895 | ||
| 896 | ! We conditionally perform work on tracer points: calculating the PLM slope, | |
| 897 | ! and updating tracer concentration within a cell | |
| 898 | ! this depends on whether there is a flux which would affect this tracer point, | |
| 899 | ! as indicated by domore_v. In the case of PPM reconstruction, a flux requires | |
| 900 | ! slope calculations at the two tracer points on either side (as indicated by | |
| 901 | ! the stencil variable), so we account for this with the do_j_tr flag array | |
| 902 | ! | |
| 903 | ! Note: this does lead to unnecessary work in updating tracer concentrations, | |
| 904 | ! since that doesn't need a wider stencil with the PPM advection scheme, but | |
| 905 | ! this would require an additional loop, etc. | |
| 906 | 921 | do concurrent (j=SZJ_(G)) |
| 907 | 63549 | do_j_tr(j) = .false. |
| 908 | enddo | |
| 909 | 57102 | do concurrent (J=js-1:je, domore_v(J,k)) |
| 910 | 57102 | do concurrent (j2=1-stencil:stencil) |
| 911 | 159300 | do_j_tr(j+j2) = .true. |
| 912 | enddo | |
| 913 | enddo | |
| 914 | 921 | do concurrent (j=SZJB_(G)) |
| 915 | 64470 | domore_v_initial(j) = domore_v(j,k) |
| 916 | enddo | |
| 917 | ||
| 918 | ! Calculate the j-direction profiles (slopes) of each tracer that | |
| 919 | ! is being advected. | |
| 920 | 921 | if (usePLMslope) then |
| 921 | 0 | do concurrent (j=js-stencil:je+stencil, do_j_tr(j)) |
| 922 | 0 | do concurrent (m=1:ntr, i=is:ie) |
| 923 | !if (ABS(Tr(m)%t(i,j+1,k)-Tr(m)%t(i,j,k)) < & | |
| 924 | ! ABS(Tr(m)%t(i,j,k)-Tr(m)%t(i,j-1,k))) then | |
| 925 | ! maxslope = 4.0*(Tr(m)%t(i,j+1,k)-Tr(m)%t(i,j,k)) | |
| 926 | !else | |
| 927 | ! maxslope = 4.0*(Tr(m)%t(i,j,k)-Tr(m)%t(i,j-1,k)) | |
| 928 | !endif | |
| 929 | !if ((Tr(m)%t(i,j+1,k)-Tr(m)%t(i,j,k))*(Tr(m)%t(i,j,k)-Tr(m)%t(i,j-1,k)) < 0.0) then | |
| 930 | ! slope_y(i,m,j) = 0.0 | |
| 931 | !elseif (ABS(Tr(m)%t(i,j+1,k)-Tr(m)%t(i,j-1,k))<ABS(maxslope)) then | |
| 932 | ! slope_y(i,m,j) = G%mask2dCv(i,J) * G%mask2dCv(i,J-1) * & | |
| 933 | ! 0.5*(Tr(m)%t(i,j+1,k)-Tr(m)%t(i,j-1,k)) | |
| 934 | !else | |
| 935 | ! slope_y(i,m,j) = G%mask2dCv(i,J) * G%mask2dCv(i,J-1) * 0.5*maxslope | |
| 936 | !endif | |
| 937 | 0 | Tp = Tr(m)%t(i,j+1,k) ; Tc = Tr(m)%t(i,j,k) ; Tm = Tr(m)%t(i,j-1,k) |
| 938 | 0 | dMx = max( Tp, Tc, Tm ) - Tc |
| 939 | 0 | dMn = Tc - min( Tp, Tc, Tm ) |
| 940 | slope_y(i,m,j) = G%mask2dCv(i,J)*G%mask2dCv(i,J-1) * & | |
| 941 | 0 | sign( min(0.5*abs(Tp-Tm), 2.0*dMx, 2.0*dMn), Tp-Tm ) |
| 942 | enddo | |
| 943 | enddo ! End of i-, m-, & j- loops. | |
| 944 | endif ! usePLMslope | |
| 945 | ||
| 946 | ||
| 947 | ! make a copy of the tracers in case values need to be overridden for OBCs | |
| 948 | ||
| 949 | 354585 | do concurrent (j=G%jsd:G%jed, m=1:ntr, i=G%isd:G%ied) |
| 950 | 24521625 | T_tmp(i,m,j) = Tr(m)%t(i,j,k) |
| 951 | enddo | |
| 952 | ||
| 953 | ! loop through open boundaries and recalculate flux terms | |
| 954 | 921 | if (associated(OBC)) then ; if (OBC%OBC_pe) then |
| 955 | 0 | do n=1,OBC%number_of_segments |
| 956 | 0 | segment=>OBC%segment(n) |
| 957 | 0 | if (.not. associated(segment%tr_Reg)) cycle |
| 958 | 0 | do i=is,ie |
| 959 | 0 | if (segment%is_N_or_S) then |
| 960 | 0 | if (i>=segment%HI%isd .and. i<=segment%HI%ied) then |
| 961 | 0 | J = segment%HI%JsdB |
| 962 | 0 | do m = 1,segment%tr_Reg%ntseg ! replace tracers with OBC values |
| 963 | 0 | ntr_id = segment%tr_reg%Tr(m)%ntr_index |
| 964 | 0 | if (allocated(segment%tr_Reg%Tr(m)%tres)) then |
| 965 | 0 | if (segment%direction == OBC_DIRECTION_S) then |
| 966 | 0 | T_tmp(i,ntr_id,j) = segment%tr_Reg%Tr(m)%tres(i,j,k) |
| 967 | else | |
| 968 | 0 | T_tmp(i,ntr_id,j+1) = segment%tr_Reg%Tr(m)%tres(i,j,k) |
| 969 | endif | |
| 970 | else | |
| 971 | 0 | if (segment%direction == OBC_DIRECTION_S) then |
| 972 | 0 | T_tmp(i,ntr_id,j) = segment%tr_Reg%Tr(m)%OBC_inflow_conc |
| 973 | else | |
| 974 | 0 | T_tmp(i,ntr_id,j+1) = segment%tr_Reg%Tr(m)%OBC_inflow_conc |
| 975 | endif | |
| 976 | endif | |
| 977 | enddo | |
| 978 | 0 | do m = 1,ntr ! Apply update tracer values for slope calculation |
| 979 | 0 | do j=segment%HI%JsdB-1,segment%HI%JsdB+1 |
| 980 | 0 | Tp = T_tmp(i,m,j+1) ; Tc = T_tmp(i,m,j) ; Tm = T_tmp(i,m,j-1) |
| 981 | 0 | dMx = max( Tp, Tc, Tm ) - Tc |
| 982 | 0 | dMn= Tc - min( Tp, Tc, Tm ) |
| 983 | slope_y(i,m,j) = G%mask2dCv(i,J)*G%mask2dCv(i,J-1) * & | |
| 984 | 0 | sign( min(0.5*abs(Tp-Tm), 2.0*dMx, 2.0*dMn), Tp-Tm ) |
| 985 | enddo | |
| 986 | enddo | |
| 987 | endif | |
| 988 | endif ! is_N_S | |
| 989 | enddo ! i-loop | |
| 990 | enddo ! segment loop | |
| 991 | endif ; endif | |
| 992 | ||
| 993 | ! Calculate the j-direction fluxes of each tracer, using as much | |
| 994 | ! the minimum of the remaining mass flux (vhr) and the half the mass | |
| 995 | ! in the cell plus whatever part of its half of the mass flux that | |
| 996 | ! the flux through the other side does not require. | |
| 997 | !$omp target teams loop & | |
| 998 | !$omp private(CFL, Ihnew, hup, hlos, j_up, Tp, Tc, Tm, i, m, aL, aR, dA, mA, a6, n, ntr_id, & | |
| 999 | !$omp domore_v_jk) | |
| 1000 | 57102 | do J=js-1,je ; if (domore_v(J,k)) then |
| 1001 | 53100 | domore_v_jk = .false. |
| 1002 | ||
| 1003 | 53100 | do concurrent (i=is:ie) DO_LOCALITY(reduce(.or.:domore_v_jk)) |
| 1004 | if ((vhr(i,J,k) == 0.0) .or. & | |
| 1005 | 6372000 | ((vhr(i,J,k) < 0.0) .and. (hprev(i,j+1,k) <= tiny_h)) .or. & |
| 1006 | 6425100 | ((vhr(i,J,k) > 0.0) .and. (hprev(i,j,k) <= tiny_h)) ) then |
| 1007 | 1967400 | vhh(i,J) = 0.0 |
| 1008 | 1967400 | CFL(i) = 0.0 |
| 1009 | 4404600 | elseif (vhr(i,J,k) < 0.0) then |
| 1010 | 3101027 | hup = hprev(i,j+1,k) - G%areaT(i,j+1)*min_h |
| 1011 | 3101027 | hlos = MAX(0.0, vhr(i,J+1,k)) |
| 1012 | 3101027 | if ((((hup - hlos) + vhr(i,J,k)) < 0.0) .and. & |
| 1013 | ((0.5*hup + vhr(i,J,k)) < 0.0)) then | |
| 1014 | 0 | vhh(i,J) = MIN(-0.5*hup, -hup+hlos, 0.0) |
| 1015 | 0 | domore_v_jk = .true. |
| 1016 | else | |
| 1017 | 3101027 | vhh(i,J) = vhr(i,J,k) |
| 1018 | endif | |
| 1019 | 3101027 | CFL(i) = - vhh(i,J) / hprev(i,j+1,k) ! CFL is positive |
| 1020 | else | |
| 1021 | 1303573 | hup = hprev(i,j,k) - G%areaT(i,j)*min_h |
| 1022 | 1303573 | hlos = MAX(0.0, -vhr(i,J-1,k)) |
| 1023 | 1303573 | if ((((hup - hlos) - vhr(i,J,k)) < 0.0) .and. & |
| 1024 | ((0.5*hup - vhr(i,J,k)) < 0.0)) then | |
| 1025 | 0 | vhh(i,J) = MAX(0.5*hup, hup-hlos, 0.0) |
| 1026 | 0 | domore_v_jk = .true. |
| 1027 | else | |
| 1028 | 1303573 | vhh(i,J) = vhr(i,J,k) |
| 1029 | endif | |
| 1030 | 1303573 | CFL(i) = vhh(i,J) / hprev(i,j,k) ! CFL is positive |
| 1031 | endif | |
| 1032 | enddo | |
| 1033 | ||
| 1034 | 53100 | domore_v(j,k) = domore_v_jk |
| 1035 | ||
| 1036 | 53100 | do concurrent (m=1:ntr) |
| 1037 | ||
| 1038 | 212400 | if ((advect_schemes(m) == ADVECT_PPM) .or. (advect_schemes(m) == ADVECT_PPMH3)) then |
| 1039 | 159300 | do concurrent (i=is:ie) |
| 1040 | ! centre cell depending on upstream direction | |
| 1041 | 19116000 | if (vhh(i,J) >= 0.0) then |
| 1042 | 9812919 | j_up = j |
| 1043 | else | |
| 1044 | 9303081 | j_up = j + 1 |
| 1045 | endif | |
| 1046 | ||
| 1047 | ! Implementation of PPM-H3 | |
| 1048 | 19116000 | Tp = T_tmp(i,m,j_up+1) ; Tc = T_tmp(i,m,j_up) ; Tm = T_tmp(i,m,j_up-1) |
| 1049 | ||
| 1050 | 19116000 | if (advect_schemes(m) == ADVECT_PPMH3) then |
| 1051 | 19116000 | aL = ( 5.*Tc + ( 2.*Tm - Tp ) )/6. ! H3 estimate |
| 1052 | 19116000 | aL = max( min(Tc,Tm), aL) ; aL = min( max(Tc,Tm), aL) ! Bound |
| 1053 | 19116000 | aR = ( 5.*Tc + ( 2.*Tp - Tm ) )/6. ! H3 estimate |
| 1054 | 19116000 | aR = max( min(Tc,Tp), aR) ; aR = min( max(Tc,Tp), aR) ! Bound |
| 1055 | else | |
| 1056 | 0 | aL = 0.5 * ((Tm + Tc) + (slope_y(i,m,j_up-1) - slope_y(i,m,j_up)) / 3.) |
| 1057 | 0 | aR = 0.5 * ((Tc + Tp) + (slope_y(i,m,j_up) - slope_y(i,m,j_up+1)) / 3.) |
| 1058 | endif | |
| 1059 | ||
| 1060 | 19116000 | dA = aR - aL ; mA = 0.5*( aR + aL ) |
| 1061 | 19116000 | if (G%mask2dCv(i,J_up)*G%mask2dCv(i,J_up-1)*(Tp-Tc)*(Tc-Tm) <= 0.) then |
| 1062 | 11795025 | aL = Tc ; aR = Tc ! PCM for local extrema and boundary cells |
| 1063 | 7320975 | elseif ( dA*(Tc-mA) > (dA*dA)/6. ) then |
| 1064 | 547064 | aL = (3.*Tc) - 2.*aR |
| 1065 | 6773911 | elseif ( dA*(Tc-mA) < - (dA*dA)/6. ) then |
| 1066 | 574778 | aR = (3.*Tc) - 2.*aL |
| 1067 | endif | |
| 1068 | ||
| 1069 | 19116000 | a6 = 6.*Tc - 3. * (aR + aL) ! Curvature |
| 1070 | ||
| 1071 | 38391300 | if (vhh(i,J) >= 0.0) then |
| 1072 | flux_y(i,m,J) = vhh(i,J)*( aR - 0.5 * CFL(i) * ( & | |
| 1073 | 9812919 | ( aR - aL ) - a6 * ( 1. - 2./3. * CFL(I) ) ) ) |
| 1074 | else | |
| 1075 | flux_y(i,m,J) = vhh(i,J)*( aL + 0.5 * CFL(i) * ( & | |
| 1076 | 9303081 | ( aR - aL ) + a6 * ( 1. - 2./3. * CFL(I) ) ) ) |
| 1077 | endif | |
| 1078 | enddo | |
| 1079 | else ! PLM | |
| 1080 | 0 | do concurrent (i=is:ie) |
| 1081 | 0 | if (vhh(i,J) >= 0.0) then |
| 1082 | ! Indirect implementation of PLM | |
| 1083 | !aL = Tr(m)%t(i,j,k) - 0.5 * slope_y(i,m,j) | |
| 1084 | !aR = Tr(m)%t(i,j,k) + 0.5 * slope_y(i,m,j) | |
| 1085 | !flux_y(i,m,J) = vhh(i,J)*( aR - 0.5 * (aR-aL) * CFL(i) ) | |
| 1086 | ! Alternative implementation of PLM | |
| 1087 | 0 | Tc = T_tmp(i,m,j) |
| 1088 | 0 | flux_y(i,m,J) = vhh(i,J)*( Tc + 0.5 * slope_y(i,m,j) * ( 1. - CFL(i) ) ) |
| 1089 | else | |
| 1090 | ! Indirect implementation of PLM | |
| 1091 | !aL = Tr(m)%t(i,j+1,k) - 0.5 * slope_y(i,m,j+1) | |
| 1092 | !aR = Tr(m)%t(i,j+1,k) + 0.5 * slope_y(i,m,j+1) | |
| 1093 | !flux_y(i,m,J) = vhh(i,J)*( aL + 0.5 * (aR-aL) * CFL(i) ) | |
| 1094 | ! Alternative implementation of PLM | |
| 1095 | 0 | Tc = T_tmp(i,m,j+1) |
| 1096 | 0 | flux_y(i,m,J) = vhh(i,J)*( Tc - 0.5 * slope_y(i,m,j+1) * ( 1. - CFL(i) ) ) |
| 1097 | endif | |
| 1098 | enddo | |
| 1099 | endif ! usePPM | |
| 1100 | enddo | |
| 1101 | ||
| 1102 | 53100 | if (associated(OBC)) then ; if (OBC%OBC_pe) then |
| 1103 | 0 | if (OBC%specified_v_BCs_exist_globally .or. OBC%open_v_BCs_exist_globally) then |
| 1104 | 0 | do n=1,OBC%number_of_segments |
| 1105 | 0 | if (.not. OBC%segment(n)%specified) cycle |
| 1106 | 0 | if (.not. associated(OBC%segment(n)%tr_Reg)) cycle |
| 1107 | 0 | if (OBC%segment(n)%is_N_or_S) then |
| 1108 | 0 | if (J >= OBC%segment(n)%HI%JsdB .and. J<= OBC%segment(n)%HI%JedB) then |
| 1109 | 0 | do i=OBC%segment(n)%HI%isd,OBC%segment(n)%HI%ied |
| 1110 | ! Tracer fluxes are set to prescribed values only for inflows from masked areas. | |
| 1111 | ! Now changing to simply fixed inflows. | |
| 1112 | 0 | if ((vhr(i,J,k) > 0.0) .and. (OBC%segment(n)%direction == OBC_DIRECTION_S) .or. & |
| 1113 | 0 | (vhr(i,J,k) < 0.0) .and. (OBC%segment(n)%direction == OBC_DIRECTION_N)) then |
| 1114 | 0 | vhh(i,J) = vhr(i,J,k) |
| 1115 | 0 | do m=1,OBC%segment(n)%tr_Reg%ntseg |
| 1116 | 0 | ntr_id = OBC%segment(n)%tr_reg%Tr(m)%ntr_index |
| 1117 | 0 | if (allocated(OBC%segment(n)%tr_Reg%Tr(m)%tres)) then |
| 1118 | 0 | flux_y(i,ntr_id,J) = vhh(i,J)*OBC%segment(n)%tr_Reg%Tr(m)%tres(i,J,k) |
| 1119 | else | |
| 1120 | 0 | flux_y(i,ntr_id,J) = vhh(i,J)*OBC%segment(n)%tr_Reg%Tr(m)%OBC_inflow_conc |
| 1121 | endif | |
| 1122 | enddo | |
| 1123 | endif | |
| 1124 | enddo | |
| 1125 | endif | |
| 1126 | endif | |
| 1127 | enddo | |
| 1128 | endif | |
| 1129 | ||
| 1130 | 0 | if (OBC%open_v_BCs_exist_globally) then |
| 1131 | 0 | do n=1,OBC%number_of_segments |
| 1132 | 0 | if (OBC%segment(n)%specified) cycle |
| 1133 | 0 | if (.not. associated(OBC%segment(n)%tr_Reg)) cycle |
| 1134 | 0 | if (OBC%segment(n)%is_N_or_S .and. (J >= OBC%segment(n)%HI%JsdB .and. J<= OBC%segment(n)%HI%JedB)) then |
| 1135 | 0 | do i=OBC%segment(n)%HI%isd,OBC%segment(n)%HI%ied |
| 1136 | ! Tracer fluxes are set to prescribed values only for inflows from masked areas. | |
| 1137 | 0 | if ((vhr(i,J,k) > 0.0) .and. (G%mask2dT(i,j) < 0.5) .or. & |
| 1138 | 0 | (vhr(i,J,k) < 0.0) .and. (G%mask2dT(i,j+1) < 0.5)) then |
| 1139 | 0 | vhh(i,J) = vhr(i,J,k) |
| 1140 | 0 | do m=1,OBC%segment(n)%tr_Reg%ntseg |
| 1141 | 0 | ntr_id = OBC%segment(n)%tr_reg%Tr(m)%ntr_index |
| 1142 | 0 | if (allocated(OBC%segment(n)%tr_Reg%Tr(m)%tres)) then |
| 1143 | 0 | flux_y(i,ntr_id,J) = vhh(i,J)*OBC%segment(n)%tr_Reg%Tr(m)%tres(i,J,k) |
| 1144 | 0 | else ; flux_y(i,ntr_id,J) = vhh(i,J)*OBC%segment(n)%tr_Reg%Tr(m)%OBC_inflow_conc ; endif |
| 1145 | enddo | |
| 1146 | endif | |
| 1147 | enddo | |
| 1148 | endif | |
| 1149 | enddo | |
| 1150 | endif | |
| 1151 | endif ; endif | |
| 1152 | ||
| 1153 | else ! not domore_v. | |
| 1154 | 3081 | do concurrent (i=is:ie) |
| 1155 | 372801 | vhh(i,J) = 0.0 |
| 1156 | enddo | |
| 1157 | 3081 | do concurrent (m=1:ntr, i=is:ie) |
| 1158 | 1481961 | flux_y(i,m,J) = 0.0 |
| 1159 | enddo | |
| 1160 | endif ; enddo ! End of j-loop | |
| 1161 | ||
| 1162 | 111441 | do concurrent (J=js-1:je, i=is:ie) |
| 1163 | 6741720 | vhr(i,J,k) = vhr(i,J,k) - vhh(i,J) |
| 1164 | 6853161 | if (abs(vhr(i,J,k)) < vh_neglect(i,J)) vhr(i,J,k) = 0.0 |
| 1165 | enddo | |
| 1166 | ||
| 1167 | ! Calculate new tracer concentration in each cell after accounting | |
| 1168 | ! for the j-direction fluxes. | |
| 1169 | !$omp target teams loop private(hlst, Ihnew, i, m) | |
| 1170 | 56181 | do j=js,je ; if (do_j_tr(j)) then |
| 1171 | 54000 | do concurrent (i=is:ie) |
| 1172 | 6534000 | if ((vhh(i,J) /= 0.0) .or. (vhh(i,J-1) /= 0.0)) then |
| 1173 | 4512600 | do_i(i,j) = .true. |
| 1174 | 4512600 | hlst(i) = hprev(i,j,k) |
| 1175 | 4512600 | hprev(i,j,k) = max(hprev(i,j,k) - (vhh(i,J) - vhh(i,J-1)), 0.0) |
| 1176 | 4512600 | if (hprev(i,j,k) <= 0.0) then ; do_i(i,j) = .false. |
| 1177 | 4512600 | elseif (hprev(i,j,k) < h_neglect*G%areaT(i,j)) then |
| 1178 | 0 | hlst(i) = hlst(i) + (h_neglect*G%areaT(i,j) - hprev(i,j,k)) |
| 1179 | 0 | Ihnew(i) = 1.0 / (h_neglect*G%areaT(i,j)) |
| 1180 | 4512600 | else ; Ihnew(i) = 1.0 / hprev(i,j,k) ; endif |
| 1181 | 1967400 | else ; do_i(i,j) = .false. ; endif |
| 1182 | enddo | |
| 1183 | ||
| 1184 | ! Update do_i so that nothing changes outside of the OBC (problem for interior OBCs only) | |
| 1185 | 54000 | if (associated(OBC)) then |
| 1186 | if (.not.OBC%exterior_OBC_bug .and. OBC%OBC_pe & | |
| 1187 | 0 | .and. (OBC%specified_v_BCs_exist_globally .or. OBC%open_v_BCs_exist_globally)) then |
| 1188 | ! OBC_DIRECTION_N / OBC_DIRECTION_S on the south / north edge | |
| 1189 | 0 | do concurrent (i=is:ie, OBC%segnum_v(i,J-1) > 0 .or. OBC%segnum_v(i,J) < 0) |
| 1190 | 0 | do_i(i,j) = .false. |
| 1191 | enddo | |
| 1192 | endif | |
| 1193 | endif | |
| 1194 | ||
| 1195 | ! update tracer and save some diagnostics | |
| 1196 | 54000 | do concurrent (m=1:ntr) |
| 1197 | 19602000 | do concurrent (i=is:ie, do_i(i,j)) |
| 1198 | Tr(m)%t(i,j,k) = (Tr(m)%t(i,j,k) * hlst(i) - & | |
| 1199 | 19602000 | (flux_y(i,m,J) - flux_y(i,m,J-1))) * Ihnew(i) |
| 1200 | enddo | |
| 1201 | ||
| 1202 | ! diagnose convergence of flux_y and add to convergence of flux_x. | |
| 1203 | ! division by areaT to get into W/m2 for heat and kg/(s*m2) for salt. | |
| 1204 | 216000 | if (associated(Tr(m)%advection_xy)) then |
| 1205 | 0 | do concurrent (i=is:ie, do_i(i,j)) |
| 1206 | Tr(m)%advection_xy(i,j,k) = Tr(m)%advection_xy(i,j,k) - & | |
| 1207 | (flux_y(i,m,J) - flux_y(i,m,J-1))* Idt * & | |
| 1208 | 0 | G%IareaT(i,j) |
| 1209 | enddo | |
| 1210 | endif | |
| 1211 | ||
| 1212 | enddo | |
| 1213 | endif ; enddo ! End of j-loop. | |
| 1214 | ||
| 1215 | ! Do user controlled underflow of the tracer concentrations. | |
| 1216 | 3684 | do concurrent (m=1:ntr, Tr(m)%conc_underflow > 0.0) |
| 1217 | 3684 | do concurrent (j=js:je, i=is:ie, abs(Tr(m)%t(i,j,k)) < Tr(m)%conc_underflow) |
| 1218 | 0 | Tr(m)%t(i,j,k) = 0.0 |
| 1219 | enddo | |
| 1220 | enddo | |
| 1221 | ||
| 1222 | ! compute ad_y and ad2d_y diagnostic outside above j-loop so as to make the summation ordered when OMP is active. | |
| 1223 | 3684 | do m=1,ntr ; if (associated(Tr(m)%ad_y)) then |
| 1224 | 0 | do concurrent (J=js-1:je, domore_v_initial(J)) |
| 1225 | 0 | do concurrent (i=is:ie, do_i(i,j) .or. do_i(i,j+1)) |
| 1226 | 0 | Tr(m)%ad_y(i,J,k) = Tr(m)%ad_y(i,J,k) + flux_y(i,m,J)*Idt |
| 1227 | enddo | |
| 1228 | enddo | |
| 1229 | endif ; enddo ! End of m-loop. | |
| 1230 | ||
| 1231 | 3684 | do m=1,ntr ; if (associated(Tr(m)%ad2d_y)) then |
| 1232 | 0 | do concurrent (J=js-1:je, domore_v_initial(J)) |
| 1233 | 0 | do concurrent (i=is:ie, do_i(i,j) .or. do_i(i,j+1)) |
| 1234 | 0 | Tr(m)%ad2d_y(i,J) = Tr(m)%ad2d_y(i,J) + flux_y(i,m,J)*Idt |
| 1235 | enddo | |
| 1236 | enddo | |
| 1237 | endif ; enddo ! End of m-loop. | |
| 1238 | ||
| 1239 | !$omp target exit data map(release: vhh, T_tmp, slope_y, flux_y, domore_v_initial, do_j_tr, do_i) | |
| 1240 | ||
| 1241 | 921 | end subroutine advect_y |
| 1242 | ||
| 1243 | !> Initialize lateral tracer advection module | |
| 1244 | 1 | subroutine tracer_advect_init(Time, G, US, param_file, diag, CS) |
| 1245 | type(time_type), target, intent(in) :: Time !< current model time | |
| 1246 | type(ocean_grid_type), intent(in) :: G !< ocean grid structure | |
| 1247 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 1248 | type(param_file_type), intent(in) :: param_file !< open file to parse for model parameters | |
| 1249 | type(diag_ctrl), target, intent(inout) :: diag !< regulates diagnostic output | |
| 1250 | type(tracer_advect_CS), pointer :: CS !< module control structure | |
| 1251 | ||
| 1252 | ! This include declares and sets the variable "version". | |
| 1253 | # include "version_variable.h" | |
| 1254 | character(len=40) :: mdl = "MOM_tracer_advect" ! This module's name. | |
| 1255 | character(len=256) :: mesg ! Message for error messages. | |
| 1256 | ||
| 1257 | 1 | if (associated(CS)) then |
| 1258 | 0 | call MOM_error(WARNING, "tracer_advect_init called with associated control structure.") |
| 1259 | 0 | return |
| 1260 | endif | |
| 1261 | 1 | allocate(CS) |
| 1262 | ||
| 1263 | 1 | CS%diag => diag |
| 1264 | ||
| 1265 | ! Read all relevant parameters and write them to the model log. | |
| 1266 | 1 | call log_version(param_file, mdl, version, "") |
| 1267 | call get_param(param_file, mdl, "DT", CS%dt, fail_if_missing=.true., & | |
| 1268 | 1 | desc="The (baroclinic) dynamics time step.", units="s", scale=US%s_to_T) |
| 1269 | 1 | call get_param(param_file, mdl, "DEBUG", CS%debug, default=.false.) |
| 1270 | call get_param(param_file, mdl, "TRACER_ADVECTION_SCHEME", mesg, & | |
| 1271 | desc="The horizontal transport scheme for tracers:\n"//& | |
| 1272 | 1 | trim(TracerAdvectionSchemeDoc), default='PLM') |
| 1273 | ||
| 1274 | ! Get the integer value of the tracer scheme | |
| 1275 | 1 | call set_tracer_advect_scheme(CS%default_advect_scheme, mesg) |
| 1276 | ||
| 1277 | 1 | if (CS%default_advect_scheme == ADVECT_PPMH3) then |
| 1278 | call get_param(param_file, mdl, "USE_HUYNH_STENCIL_BUG", & | |
| 1279 | CS%useHuynhStencilBug, & | |
| 1280 | desc="If true, use a stencil width of 2 in PPM:H3 tracer advection. " & | |
| 1281 | // "This is incorrect and will produce regressions in certain " & | |
| 1282 | // "configurations, but may be required to reproduce results in " & | |
| 1283 | // "legacy simulations.", & | |
| 1284 | 1 | default=.false.) |
| 1285 | endif | |
| 1286 | ||
| 1287 | !$omp target enter data map(to: CS) | |
| 1288 | ||
| 1289 | 1 | id_clock_advect = cpu_clock_id('(Ocean advect tracer)', grain=CLOCK_MODULE) |
| 1290 | 1 | id_clock_pass = cpu_clock_id('(Ocean tracer halo updates)', grain=CLOCK_ROUTINE) |
| 1291 | 1 | id_clock_sync = cpu_clock_id('(Ocean tracer global synch)', grain=CLOCK_ROUTINE) |
| 1292 | ||
| 1293 | end subroutine tracer_advect_init | |
| 1294 | ||
| 1295 | !> Close the tracer advection module | |
| 1296 | 1 | subroutine tracer_advect_end(CS) |
| 1297 | type(tracer_advect_CS), pointer :: CS !< module control structure | |
| 1298 | ||
| 1299 | !$omp target exit data map(delete: CS) | |
| 1300 | 1 | if (associated(CS)) deallocate(CS) |
| 1301 | ||
| 1302 | 1 | end subroutine tracer_advect_end |
| 1303 | ||
| 1304 | ||
| 1305 | !> \namespace mom_tracer_advect | |
| 1306 | !! | |
| 1307 | !! This program contains the subroutines that advect tracers | |
| 1308 | !! horizontally (i.e. along layers). | |
| 1309 | !! | |
| 1310 | !! \section section_mom_advect_intro | |
| 1311 | !! | |
| 1312 | !! * advect_tracer advects tracer concentrations using a combination | |
| 1313 | !! of the modified flux advection scheme from Easter (Mon. Wea. Rev., | |
| 1314 | !! 1993) with tracer distributions given by the monotonic piecewise | |
| 1315 | !! parabolic method, as described in Carpenter et al. (MWR, 1990). | |
| 1316 | !! This scheme conserves the total amount of tracer while avoiding | |
| 1317 | !! spurious maxima and minima of the tracer concentration. | |
| 1318 | !! | |
| 1319 | !! * advect_tracer subroutine determines the volume of a layer in | |
| 1320 | !! a grid cell at the previous instance when the tracer concentration | |
| 1321 | !! was changed, so it is essential that the volume fluxes should be | |
| 1322 | !! correct. It is also important that the tracer advection occurs | |
| 1323 | !! before each calculation of the diabatic forcing. | |
| 1324 | !! | |
| 1325 | !! The advection scheme of some tracers can be set to be different | |
| 1326 | !! to that used by active tracers. | |
| 1327 | ||
| 1328 | ||
| 1329 | 0 | end module MOM_tracer_advect |