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 | #include "do_concurrent_compat.h" | |
| 5 | ||
| 6 | !> Main routine for lateral (along surface or neutral) diffusion of tracers | |
| 7 | module MOM_tracer_hor_diff | |
| 8 | ||
| 9 | use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end | |
| 10 | use MOM_cpu_clock, only : CLOCK_MODULE, CLOCK_ROUTINE | |
| 11 | use MOM_diag_mediator, only : post_data, diag_ctrl | |
| 12 | use MOM_diag_mediator, only : register_diag_field, safe_alloc_ptr, time_type | |
| 13 | use MOM_domains, only : sum_across_PEs, max_across_PEs | |
| 14 | use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type | |
| 15 | use MOM_domains, only : pass_vector | |
| 16 | use MOM_debugging, only : hchksum, uvchksum | |
| 17 | use MOM_diabatic_driver, only : diabatic_CS | |
| 18 | use MOM_EOS, only : calculate_density, EOS_type, EOS_domain | |
| 19 | use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg, is_root_pe | |
| 20 | use MOM_error_handler, only : MOM_set_verbosity, callTree_showQuery | |
| 21 | use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint | |
| 22 | use MOM_file_parser, only : get_param, log_version, param_file_type | |
| 23 | use MOM_grid, only : ocean_grid_type | |
| 24 | use MOM_lateral_mixing_coeffs, only : VarMix_CS | |
| 25 | use MOM_MEKE_types, only : MEKE_type | |
| 26 | use MOM_neutral_diffusion, only : neutral_diffusion_init, neutral_diffusion_end | |
| 27 | use MOM_neutral_diffusion, only : neutral_diffusion_CS | |
| 28 | use MOM_neutral_diffusion, only : neutral_diffusion_calc_coeffs, neutral_diffusion | |
| 29 | use MOM_hor_bnd_diffusion, only : hbd_CS, hor_bnd_diffusion_init | |
| 30 | use MOM_hor_bnd_diffusion, only : hor_bnd_diffusion, hor_bnd_diffusion_end | |
| 31 | use MOM_tracer_registry, only : tracer_registry_type, tracer_type, MOM_tracer_chksum | |
| 32 | use MOM_unit_scaling, only : unit_scale_type | |
| 33 | use MOM_variables, only : thermo_var_ptrs, vertvisc_type | |
| 34 | use MOM_verticalGrid, only : verticalGrid_type | |
| 35 | ||
| 36 | implicit none ; private | |
| 37 | ||
| 38 | #include <MOM_memory.h> | |
| 39 | ||
| 40 | public tracer_hordiff, tracer_hor_diff_init, tracer_hor_diff_end | |
| 41 | ||
| 42 | !> The control structure for along-layer and epineutral tracer diffusion | |
| 43 | type, public :: tracer_hor_diff_CS ; private | |
| 44 | real :: KhTr !< The along-isopycnal tracer diffusivity [L2 T-1 ~> m2 s-1]. | |
| 45 | real :: KhTr_Slope_Cff !< The non-dimensional coefficient in KhTr formula [nondim] | |
| 46 | real :: KhTr_min !< Minimum along-isopycnal tracer diffusivity [L2 T-1 ~> m2 s-1]. | |
| 47 | real :: KhTr_max !< Maximum along-isopycnal tracer diffusivity [L2 T-1 ~> m2 s-1]. | |
| 48 | real :: KhTr_passivity_coeff !< Passivity coefficient that scales Rd/dx (default = 0) | |
| 49 | !! where passivity is the ratio between along-isopycnal | |
| 50 | !! tracer mixing and thickness mixing [nondim] | |
| 51 | real :: KhTr_passivity_min !< Passivity minimum (default = 1/2) [nondim] | |
| 52 | real :: ML_KhTR_scale !< With Diffuse_ML_interior, the ratio of the | |
| 53 | !! truly horizontal diffusivity in the mixed | |
| 54 | !! layer to the epipycnal diffusivity [nondim]. | |
| 55 | real :: max_diff_CFL !< If positive, locally limit the along-isopycnal | |
| 56 | !! tracer diffusivity to keep the diffusive CFL | |
| 57 | !! locally at or below this value [nondim]. | |
| 58 | logical :: KhTr_use_vert_struct !< If true, uses the equivalent barotropic structure | |
| 59 | !! as the vertical structure of tracer diffusivity. | |
| 60 | logical :: full_depth_khtr_min !< If true, KHTR_MIN is enforced throughout the whole water column. | |
| 61 | !! Otherwise, KHTR_MIN is only enforced at the surface. This parameter | |
| 62 | !! is only available when KHTR_USE_EBT_STRUCT=True and KHTR_MIN>0. | |
| 63 | logical :: Diffuse_ML_interior !< If true, diffuse along isopycnals between | |
| 64 | !! the mixed layer and the interior. | |
| 65 | logical :: check_diffusive_CFL !< If true, automatically iterate the diffusion | |
| 66 | !! to ensure that the diffusive equivalent of | |
| 67 | !! the CFL limit is not violated. | |
| 68 | logical :: use_neutral_diffusion !< If true, use the neutral_diffusion module from within | |
| 69 | !! tracer_hor_diff. | |
| 70 | logical :: use_hor_bnd_diffusion !< If true, use the hor_bnd_diffusion module from within | |
| 71 | !! tracer_hor_diff. | |
| 72 | logical :: recalc_neutral_surf !< If true, recalculate the neutral surfaces if CFL has been | |
| 73 | !! exceeded | |
| 74 | logical :: limit_bug !< If true and the answer date is 20240330 or below, use a | |
| 75 | !! rotational symmetry breaking bug when limiting the tracer | |
| 76 | !! properties in tracer_epipycnal_ML_diff. | |
| 77 | integer :: answer_date !< The vintage of the order of arithmetic to use for the tracer | |
| 78 | !! diffusion. Values of 20240330 or below recover the answers | |
| 79 | !! from the original form of this code, while higher values use | |
| 80 | !! mathematically equivalent expressions that recover rotational symmetry | |
| 81 | !! when DIFFUSE_ML_TO_INTERIOR is true. | |
| 82 | type(neutral_diffusion_CS), pointer :: neutral_diffusion_CSp => NULL() !< Control structure for neutral diffusion. | |
| 83 | type(hbd_CS), pointer :: hor_bnd_diffusion_CSp => NULL() !< Control structure for | |
| 84 | !! horizontal boundary diffusion. | |
| 85 | type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to | |
| 86 | !! regulate the timing of diagnostic output. | |
| 87 | logical :: debug !< If true, write verbose checksums for debugging purposes. | |
| 88 | logical :: show_call_tree !< Display the call tree while running. Set by VERBOSITY level. | |
| 89 | logical :: first_call = .true. !< This is true until after the first call | |
| 90 | !>@{ Diagnostic IDs | |
| 91 | integer :: id_KhTr_u = -1 | |
| 92 | integer :: id_KhTr_v = -1 | |
| 93 | integer :: id_KhTr_h = -1 | |
| 94 | integer :: id_CFL = -1 | |
| 95 | integer :: id_khdt_x = -1 | |
| 96 | integer :: id_khdt_y = -1 | |
| 97 | !>@} | |
| 98 | ||
| 99 | type(group_pass_type) :: pass_t !< For group halo pass, used in both | |
| 100 | !! tracer_hordiff and tracer_epipycnal_ML_diff | |
| 101 | end type tracer_hor_diff_CS | |
| 102 | ||
| 103 | !> A type that can be used to create arrays of pointers to 2D arrays | |
| 104 | type p2d | |
| 105 | real, dimension(:,:), pointer :: p => NULL() !< A pointer to a 2D array of reals [various] | |
| 106 | end type p2d | |
| 107 | !> A type that can be used to create arrays of pointers to 2D integer arrays | |
| 108 | type p2di | |
| 109 | integer, dimension(:,:), pointer :: p => NULL() !< A pointer to a 2D array of integers | |
| 110 | end type p2di | |
| 111 | ||
| 112 | !>@{ CPU time clocks | |
| 113 | integer :: id_clock_diffuse, id_clock_epimix, id_clock_pass, id_clock_sync | |
| 114 | !>@} | |
| 115 | ||
| 116 | contains | |
| 117 | ||
| 118 | !> Compute along-coordinate diffusion of all tracers | |
| 119 | !! using the diffusivity in CS%KhTr, or using space-dependent diffusivity. | |
| 120 | !! Multiple iterations are used (if necessary) so that there is no limit | |
| 121 | !! on the acceptable time increment. | |
| 122 | 0 | subroutine tracer_hordiff(h, dt, MEKE, VarMix, visc, G, GV, US, CS, Reg, tv, do_online_flag, read_khdt_x, read_khdt_y) |
| 123 | type(ocean_grid_type), intent(inout) :: G !< Grid type | |
| 124 | type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure | |
| 125 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & | |
| 126 | intent(in) :: h !< Layer thickness [H ~> m or kg m-2] | |
| 127 | real, intent(in) :: dt !< time step [T ~> s] | |
| 128 | type(MEKE_type), intent(in) :: MEKE !< MEKE fields | |
| 129 | type(VarMix_CS), intent(in) :: VarMix !< Variable mixing type | |
| 130 | type(vertvisc_type), intent(in) :: visc !< Structure with vertical viscosities, | |
| 131 | !! boundary layer properties and related fields | |
| 132 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 133 | type(tracer_hor_diff_CS), pointer :: CS !< module control structure | |
| 134 | type(tracer_registry_type), pointer :: Reg !< registered tracers | |
| 135 | type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any available | |
| 136 | !! thermodynamic fields, including potential temp and | |
| 137 | !! salinity or mixed layer density. Absent fields have | |
| 138 | !! NULL ptrs, and these may (probably will) point to | |
| 139 | !! some of the same arrays as Tr does. tv is required | |
| 140 | !! for epipycnal mixing between mixed layer and the interior. | |
| 141 | ! Optional inputs for offline tracer transport | |
| 142 | logical, optional, intent(in) :: do_online_flag !< If present and true, do online | |
| 143 | !! tracer transport with stored velocities. | |
| 144 | ! The next two arguments do not appear to be used anywhere. | |
| 145 | real, dimension(SZIB_(G),SZJ_(G)), & | |
| 146 | optional, intent(in) :: read_khdt_x !< If present, these are the zonal diffusivities | |
| 147 | !! times a timestep from a previous run [L2 ~> m2] | |
| 148 | real, dimension(SZI_(G),SZJB_(G)), & | |
| 149 | optional, intent(in) :: read_khdt_y !< If present, these are the meridional diffusivities | |
| 150 | !! times a timestep from a previous run [L2 ~> m2] | |
| 151 | ||
| 152 | ||
| 153 | real, dimension(SZI_(G),SZJ_(G)) :: & | |
| 154 | 2 | Ihdxdy, & ! The inverse of the volume or mass of fluid in a layer in a |
| 155 | ! grid cell [H-1 L-2 ~> m-3 or kg-1]. | |
| 156 | 1 | CFL, & ! A diffusive CFL number for each cell [nondim]. |
| 157 | 2 | dTr ! The change in a tracer's concentration, in units of concentration [Conc]. |
| 158 | ||
| 159 | 2 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: Kh_h |
| 160 | ! The tracer diffusivity averaged to tracer points [L2 T-1 ~> m2 s-1]. | |
| 161 | real, dimension(SZIB_(G),SZJ_(G)) :: & | |
| 162 | 2 | khdt_x ! The value of Khtr*dt times the open face width divided by |
| 163 | ! the distance between adjacent tracer points [L2 ~> m2]. | |
| 164 | real, dimension(SZI_(G),SZJB_(G)) :: & | |
| 165 | 2 | khdt_y ! The value of Khtr*dt times the open face width divided by |
| 166 | ! the distance between adjacent tracer points [L2 ~> m2]. | |
| 167 | real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: & | |
| 168 | 2 | Coef_x, & ! The coefficients relating zonal tracer differences to time-integrated |
| 169 | ! fluxes, in [L2 ~> m2] for some schemes and [H L2 ~> m3 or kg] for others. | |
| 170 | 2 | Kh_u ! Tracer mixing coefficient at u-points [L2 T-1 ~> m2 s-1]. |
| 171 | real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: & | |
| 172 | 2 | Coef_y, & ! The coefficients relating meridional tracer differences to time-integrated |
| 173 | ! fluxes, in [L2 ~> m2] for some schemes and [H L2 ~> m3 or kg] for others. | |
| 174 | 1 | Kh_v ! Tracer mixing coefficient at u-points [L2 T-1 ~> m2 s-1]. |
| 175 | ||
| 176 | real :: khdt_max ! The local limiting value of khdt_x or khdt_y [L2 ~> m2]. | |
| 177 | real :: Coef_min ! The local limiting value of Coef_x or Coef_y, in [L2 ~> m2] for some | |
| 178 | ! schemes and [H L2 ~> m3 or kg] for others. | |
| 179 | real :: max_CFL ! The global maximum of the diffusive CFL number [nondim] | |
| 180 | logical :: use_VarMix, Resoln_scaled, do_online, use_Eady | |
| 181 | integer :: i, j, k, m, is, ie, js, je, nz, ntr, itt, num_itts | |
| 182 | real :: I_numitts ! The inverse of the number of iterations, num_itts [nondim] | |
| 183 | real :: scale ! The fraction of khdt_x or khdt_y that is applied in this | |
| 184 | ! layer for this iteration [nondim]. | |
| 185 | real :: Idt ! The inverse of the time step [T-1 ~> s-1]. | |
| 186 | real :: h_neglect ! A thickness that is so small it is usually lost | |
| 187 | ! in roundoff and can be neglected [H ~> m or kg m-2]. | |
| 188 | real :: Kh_loc ! The local value of Kh [L2 T-1 ~> m2 s-1]. | |
| 189 | real :: Res_Fn ! The local value of the resolution function [nondim]. | |
| 190 | real :: Rd_dx ! The local value of deformation radius over grid-spacing [nondim]. | |
| 191 | real :: normalize ! normalization used for diagnostic Kh_h [nondim]; diffusivity averaged to h-points. | |
| 192 | real :: MEKE_KhTr_fac | |
| 193 | ||
| 194 | 1 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke |
| 195 | ||
| 196 | 1 | do_online = .true. |
| 197 | 1 | if (present(do_online_flag)) do_online = do_online_flag |
| 198 | ||
| 199 | 1 | if (.not. associated(CS)) call MOM_error(FATAL, "MOM_tracer_hor_diff: "// & |
| 200 | 0 | "register_tracer must be called before tracer_hordiff.") |
| 201 | 1 | if (.not. associated(Reg)) call MOM_error(FATAL, "MOM_tracer_hor_diff: "// & |
| 202 | 0 | "register_tracer must be called before tracer_hordiff.") |
| 203 | 1 | if (Reg%ntr == 0 .or. (CS%KhTr <= 0.0 .and. .not. VarMix%use_variable_mixing)) return |
| 204 | ||
| 205 | 1 | if (CS%show_call_tree) call callTree_enter("tracer_hordiff(), MOM_tracer_hor_diff.F90") |
| 206 | ||
| 207 | 1 | call cpu_clock_begin(id_clock_diffuse) |
| 208 | ||
| 209 | !$omp target enter data map(to: Reg, Reg%Tr, CS) map(alloc: khdt_x, khdt_y, kh_u, kh_v) | |
| 210 | !$ do m = 1, Reg%ntr | |
| 211 | !$omp target enter data map(to: Reg%Tr(m)%t) | |
| 212 | !$omp target enter data map(to: Reg%Tr(m)%df_x) if(associated(Reg%Tr(m)%df_x)) | |
| 213 | !$omp target enter data map(to: Reg%Tr(m)%df_y) if(associated(Reg%Tr(m)%df_y)) | |
| 214 | !$omp target enter data map(to: Reg%Tr(m)%df2d_x) if(associated(Reg%Tr(m)%df2d_x)) | |
| 215 | !$omp target enter data map(to: Reg%Tr(m)%df2d_y) if(associated(Reg%Tr(m)%df2d_y)) | |
| 216 | !$ enddo | |
| 217 | ||
| 218 | 1 | ntr = Reg%ntr |
| 219 | 1 | Idt = 1.0 / dt |
| 220 | 1 | h_neglect = GV%H_subroundoff |
| 221 | ||
| 222 | 1 | if (CS%Diffuse_ML_interior .and. CS%first_call) then ; if (is_root_pe()) then |
| 223 | 4 | do m=1,ntr ; if (associated(Reg%Tr(m)%df_x) .or. associated(Reg%Tr(m)%df_y)) & |
| 224 | call MOM_error(WARNING, "tracer_hordiff: Tracer "//trim(Reg%Tr(m)%name)// & | |
| 225 | " has associated 3-d diffusive flux diagnostics. These are not "//& | |
| 226 | "valid when DIFFUSE_ML_TO_INTERIOR is defined. Use 2-d tracer "//& | |
| 227 | 1 | "diffusion diagnostics instead to get accurate total fluxes.") |
| 228 | enddo | |
| 229 | endif ; endif | |
| 230 | 1 | CS%first_call = .false. |
| 231 | ||
| 232 | 1 | if (CS%debug) call MOM_tracer_chksum("Before tracer diffusion ", Reg, G) |
| 233 | ||
| 234 | 1 | use_VarMix = .false. ; Resoln_scaled = .false. ; use_Eady = .false. |
| 235 | 1 | if (VarMix%use_variable_mixing) then |
| 236 | 1 | use_VarMix = VarMix%use_variable_mixing |
| 237 | 1 | Resoln_scaled = VarMix%Resoln_scaled_KhTr |
| 238 | 1 | use_Eady = CS%KhTr_Slope_Cff > 0. |
| 239 | 1 | CS%KhTr_use_vert_struct = allocated(VarMix%khtr_struct) |
| 240 | endif | |
| 241 | ||
| 242 | 1 | call cpu_clock_begin(id_clock_pass) |
| 243 | 4 | do m=1,ntr |
| 244 | 4 | call create_group_pass(CS%pass_t, Reg%Tr(m)%t(:,:,:), G%Domain) |
| 245 | enddo | |
| 246 | 1 | call cpu_clock_end(id_clock_pass) |
| 247 | ||
| 248 | 1 | if (CS%show_call_tree) call callTree_waypoint("Calculating diffusivity (tracer_hordiff)") |
| 249 | ||
| 250 | 1 | if (do_online) then |
| 251 | 1 | if (use_VarMix) then |
| 252 | 1 | MEKE_KhTr_fac = MEKE%KhTr_fac |
| 253 | !$omp target enter data map(to: VarMix, VarMix%SN_u, VarMix%L2u, VarMix%Res_fn_h, & | |
| 254 | !$omp VarMix%Rd_dx_h, MEKE, MEKE%Kh) | |
| 255 | 1 | do concurrent (j=js:je, I=is-1:ie) |
| 256 | 64980 | Kh_loc = CS%KhTr |
| 257 | 64980 | if (use_Eady) Kh_loc = Kh_loc + CS%KhTr_Slope_Cff*VarMix%L2u(I,j)*VarMix%SN_u(I,j) |
| 258 | 64980 | if (allocated(MEKE%Kh)) & |
| 259 | 0 | Kh_loc = Kh_loc + MEKE_KhTr_fac*sqrt(MEKE%Kh(i,j)*MEKE%Kh(i+1,j)) |
| 260 | 64980 | if (CS%KhTr_max > 0.) Kh_loc = min(Kh_loc, CS%KhTr_max) |
| 261 | 64980 | if (Resoln_scaled) & |
| 262 | 64980 | Kh_loc = Kh_loc * 0.5*(VarMix%Res_fn_h(i,j) + VarMix%Res_fn_h(i+1,j)) |
| 263 | 64980 | Kh_u(I,j,1) = max(Kh_loc, CS%KhTr_min) |
| 264 | 65342 | if (CS%KhTr_passivity_coeff>0.) then ! Apply passivity |
| 265 | 0 | Rd_dx=0.5*( VarMix%Rd_dx_h(i,j)+VarMix%Rd_dx_h(i+1,j) ) ! Rd/dx at u-points |
| 266 | 0 | Kh_loc = Kh_u(I,j,1)*max( CS%KhTr_passivity_min, CS%KhTr_passivity_coeff*Rd_dx ) |
| 267 | 0 | if (CS%KhTr_max > 0.) Kh_loc = min(Kh_loc, CS%KhTr_max) ! Re-apply max |
| 268 | 0 | Kh_u(I,j,1) = max(Kh_loc, CS%KhTr_min) ! Re-apply min |
| 269 | endif | |
| 270 | enddo | |
| 271 | 1 | do concurrent (J=js-1:je, i=is:ie) |
| 272 | 65160 | Kh_loc = CS%KhTr |
| 273 | 65160 | if (use_Eady) Kh_loc = Kh_loc + CS%KhTr_Slope_Cff*VarMix%L2v(i,J)*VarMix%SN_v(i,J) |
| 274 | 65160 | if (allocated(MEKE%Kh)) & |
| 275 | 0 | Kh_loc = Kh_loc + MEKE_KhTr_fac*sqrt(MEKE%Kh(i,j)*MEKE%Kh(i,j+1)) |
| 276 | 65160 | if (CS%KhTr_max > 0.) Kh_loc = min(Kh_loc, CS%KhTr_max) |
| 277 | 65160 | if (Resoln_scaled) & |
| 278 | 65160 | Kh_loc = Kh_loc * 0.5*(VarMix%Res_fn_h(i,j) + VarMix%Res_fn_h(i,j+1)) |
| 279 | 65160 | Kh_v(i,J,1) = max(Kh_loc, CS%KhTr_min) |
| 280 | 65521 | if (CS%KhTr_passivity_coeff>0.) then ! Apply passivity |
| 281 | 0 | Rd_dx = 0.5*( VarMix%Rd_dx_h(i,j)+VarMix%Rd_dx_h(i,j+1) ) ! Rd/dx at v-points |
| 282 | 0 | Kh_loc = Kh_v(i,J,1)*max( CS%KhTr_passivity_min, CS%KhTr_passivity_coeff*Rd_dx ) |
| 283 | 0 | if (CS%KhTr_max > 0.) Kh_loc = min(Kh_loc, CS%KhTr_max) ! Re-apply max |
| 284 | 0 | Kh_v(i,J,1) = max(Kh_loc, CS%KhTr_min) ! Re-apply min |
| 285 | endif | |
| 286 | enddo | |
| 287 | !$omp target exit data map(release: VarMix, VarMix%SN_u, VarMix%L2u, VarMix%SN_v, & | |
| 288 | !$omp VarMix%L2v, VarMix%Res_fn_h, VarMix%Rd_dx_h, MEKE, MEKE%Kh) | |
| 289 | ||
| 290 | 1 | do concurrent (j=js:je, I=is-1:ie) |
| 291 | 65342 | khdt_x(I,j) = dt*(Kh_u(I,j,1)*(G%dy_Cu(I,j)*G%IdxCu(I,j))) |
| 292 | enddo | |
| 293 | 1 | do concurrent (J=js-1:je, i=is:ie) |
| 294 | 65521 | khdt_y(i,J) = dt*(Kh_v(i,J,1)*(G%dx_Cv(i,J)*G%IdyCv(i,J))) |
| 295 | enddo | |
| 296 | 0 | elseif (Resoln_scaled) then |
| 297 | !$OMP parallel do default(shared) private(Res_fn) | |
| 298 | 0 | do j=js,je ; do I=is-1,ie |
| 299 | 0 | Res_fn = 0.5 * (VarMix%Res_fn_h(i,j) + VarMix%Res_fn_h(i+1,j)) |
| 300 | 0 | Kh_u(I,j,1) = max(CS%KhTr * Res_fn, CS%KhTr_min) |
| 301 | 0 | khdt_x(I,j) = dt*(CS%KhTr*(G%dy_Cu(I,j)*G%IdxCu(I,j))) * Res_fn |
| 302 | enddo ; enddo | |
| 303 | !$OMP parallel do default(shared) private(Res_fn) | |
| 304 | 0 | do J=js-1,je ; do i=is,ie |
| 305 | 0 | Res_fn = 0.5*(VarMix%Res_fn_h(i,j) + VarMix%Res_fn_h(i,j+1)) |
| 306 | 0 | Kh_v(i,J,1) = max(CS%KhTr * Res_fn, CS%KhTr_min) |
| 307 | 0 | khdt_y(i,J) = dt*(CS%KhTr*(G%dx_Cv(i,J)*G%IdyCv(i,J))) * Res_fn |
| 308 | enddo ; enddo | |
| 309 | !$omp target update to(khdt_x, khdt_y, Kh_u, Kh_v) | |
| 310 | else ! Use a simple constant diffusivity. | |
| 311 | 0 | if (CS%id_KhTr_u > 0) then |
| 312 | !$OMP parallel do default(shared) | |
| 313 | 0 | do j=js,je ; do I=is-1,ie |
| 314 | 0 | Kh_u(I,j,1) = CS%KhTr |
| 315 | 0 | khdt_x(I,j) = dt*(CS%KhTr*(G%dy_Cu(I,j)*G%IdxCu(I,j))) |
| 316 | enddo ; enddo | |
| 317 | else | |
| 318 | !$OMP parallel do default(shared) | |
| 319 | 0 | do j=js,je ; do I=is-1,ie |
| 320 | 0 | khdt_x(I,j) = dt*(CS%KhTr*(G%dy_Cu(I,j)*G%IdxCu(I,j))) |
| 321 | enddo ; enddo | |
| 322 | endif | |
| 323 | 0 | if (CS%id_KhTr_v > 0) then |
| 324 | !$OMP parallel do default(shared) | |
| 325 | 0 | do J=js-1,je ; do i=is,ie |
| 326 | 0 | Kh_v(i,J,1) = CS%KhTr |
| 327 | 0 | khdt_y(i,J) = dt*(CS%KhTr*(G%dx_Cv(i,J)*G%IdyCv(i,J))) |
| 328 | enddo ; enddo | |
| 329 | else | |
| 330 | !$OMP parallel do default(shared) | |
| 331 | 0 | do J=js-1,je ; do i=is,ie |
| 332 | 0 | khdt_y(i,J) = dt*(CS%KhTr*(G%dx_Cv(i,J)*G%IdyCv(i,J))) |
| 333 | enddo ; enddo | |
| 334 | endif | |
| 335 | !$omp target update to(khdt_x, khdt_y, Kh_u, Kh_v) | |
| 336 | endif ! VarMix | |
| 337 | ||
| 338 | 1 | if (CS%max_diff_CFL > 0.0) then |
| 339 | !$omp target update from(khdt_x, khdt_y, Kh_u, Kh_v) | |
| 340 | 0 | if ((CS%id_KhTr_u > 0) .or. (CS%id_KhTr_h > 0)) then |
| 341 | !$OMP parallel do default(shared) private(khdt_max) | |
| 342 | 0 | do j=js,je ; do I=is-1,ie |
| 343 | 0 | khdt_max = 0.125*CS%max_diff_CFL * min(G%areaT(i,j), G%areaT(i+1,j)) |
| 344 | 0 | if (khdt_x(I,j) > khdt_max) then |
| 345 | 0 | khdt_x(I,j) = khdt_max |
| 346 | 0 | if (dt*(G%dy_Cu(I,j)*G%IdxCu(I,j)) > 0.0) & |
| 347 | 0 | Kh_u(I,j,1) = khdt_x(I,j) / (dt*(G%dy_Cu(I,j)*G%IdxCu(I,j))) |
| 348 | endif | |
| 349 | enddo ; enddo | |
| 350 | else | |
| 351 | !$OMP parallel do default(shared) private(khdt_max) | |
| 352 | 0 | do j=js,je ; do I=is-1,ie |
| 353 | 0 | khdt_max = 0.125*CS%max_diff_CFL * min(G%areaT(i,j), G%areaT(i+1,j)) |
| 354 | 0 | khdt_x(I,j) = min(khdt_x(I,j), khdt_max) |
| 355 | enddo ; enddo | |
| 356 | endif | |
| 357 | 0 | if ((CS%id_KhTr_v > 0) .or. (CS%id_KhTr_h > 0)) then |
| 358 | !$OMP parallel do default(shared) private(khdt_max) | |
| 359 | 0 | do J=js-1,je ; do i=is,ie |
| 360 | 0 | khdt_max = 0.125*CS%max_diff_CFL * min(G%areaT(i,j), G%areaT(i,j+1)) |
| 361 | 0 | if (khdt_y(i,J) > khdt_max) then |
| 362 | 0 | khdt_y(i,J) = khdt_max |
| 363 | 0 | if (dt*(G%dx_Cv(i,J)*G%IdyCv(i,J)) > 0.0) & |
| 364 | 0 | Kh_v(i,J,1) = khdt_y(i,J) / (dt*(G%dx_Cv(i,J)*G%IdyCv(i,J))) |
| 365 | endif | |
| 366 | enddo ; enddo | |
| 367 | else | |
| 368 | !$OMP parallel do default(shared) private(khdt_max) | |
| 369 | 0 | do J=js-1,je ; do i=is,ie |
| 370 | 0 | khdt_max = 0.125*CS%max_diff_CFL * min(G%areaT(i,j), G%areaT(i,j+1)) |
| 371 | 0 | khdt_y(i,J) = min(khdt_y(i,J), khdt_max) |
| 372 | enddo ; enddo | |
| 373 | endif | |
| 374 | !$omp target update to(khdt_x, khdt_y, Kh_u, Kh_v) | |
| 375 | endif | |
| 376 | ||
| 377 | else ! .not. do_online | |
| 378 | !$OMP parallel do default(shared) | |
| 379 | 0 | do j=js,je ; do I=is-1,ie |
| 380 | 0 | khdt_x(I,j) = read_khdt_x(I,j) |
| 381 | enddo ; enddo | |
| 382 | !$OMP parallel do default(shared) | |
| 383 | 0 | do J=js-1,je ; do i=is,ie |
| 384 | 0 | khdt_y(i,J) = read_khdt_y(i,J) |
| 385 | enddo ; enddo | |
| 386 | 0 | call pass_vector(khdt_x, khdt_y, G%Domain) |
| 387 | !$omp target update to(khdt_x, khdt_y, Kh_u, Kh_v) | |
| 388 | endif ! do_online | |
| 389 | ||
| 390 | 1 | if (CS%check_diffusive_CFL) then |
| 391 | 0 | if (CS%show_call_tree) call callTree_waypoint("Checking diffusive CFL (tracer_hordiff)") |
| 392 | 0 | max_CFL = 0.0 |
| 393 | 0 | do j=js,je ; do i=is,ie |
| 394 | CFL(i,j) = 2.0*((khdt_x(I-1,j) + khdt_x(I,j)) + & | |
| 395 | 0 | (khdt_y(i,J-1) + khdt_y(i,J))) * G%IareaT(i,j) |
| 396 | 0 | if (max_CFL < CFL(i,j)) max_CFL = CFL(i,j) |
| 397 | enddo ; enddo | |
| 398 | 0 | call cpu_clock_begin(id_clock_sync) |
| 399 | 0 | call max_across_PEs(max_CFL) |
| 400 | 0 | call cpu_clock_end(id_clock_sync) |
| 401 | 0 | num_itts = max(1, ceiling(max_CFL - 4.0*EPSILON(max_CFL))) |
| 402 | 0 | I_numitts = 1.0 / (real(num_itts)) |
| 403 | 0 | if (CS%id_CFL > 0) call post_data(CS%id_CFL, CFL, CS%diag) |
| 404 | 1 | elseif (CS%max_diff_CFL > 0.0) then |
| 405 | 0 | num_itts = max(1, ceiling(CS%max_diff_CFL - 4.0*EPSILON(CS%max_diff_CFL))) |
| 406 | 0 | I_numitts = 1.0 / (real(num_itts)) |
| 407 | else | |
| 408 | 1 | num_itts = 1 ; I_numitts = 1.0 |
| 409 | endif | |
| 410 | ||
| 411 | 4 | do m=1,ntr |
| 412 | 3 | if (associated(Reg%Tr(m)%df_x)) then |
| 413 | 0 | do k=1,nz ; do j=js,je ; do I=is-1,ie |
| 414 | 0 | Reg%Tr(m)%df_x(I,j,k) = 0.0 |
| 415 | enddo ; enddo ; enddo | |
| 416 | !$omp target update to(Reg%Tr(m)%df_x) | |
| 417 | endif | |
| 418 | 3 | if (associated(Reg%Tr(m)%df_y)) then |
| 419 | 0 | do k=1,nz ; do J=js-1,je ; do i=is,ie |
| 420 | 0 | Reg%Tr(m)%df_y(i,J,k) = 0.0 |
| 421 | enddo ; enddo ; enddo | |
| 422 | !$omp target update to(Reg%Tr(m)%df_y) | |
| 423 | endif | |
| 424 | 3 | if (associated(Reg%Tr(m)%df2d_x)) then |
| 425 | 0 | do j=js,je ; do I=is-1,ie ; Reg%Tr(m)%df2d_x(I,j) = 0.0 ; enddo ; enddo |
| 426 | !$omp target update to(Reg%Tr(m)%df2d_x) | |
| 427 | endif | |
| 428 | 4 | if (associated(Reg%Tr(m)%df2d_y)) then |
| 429 | 0 | do J=js-1,je ; do i=is,ie ; Reg%Tr(m)%df2d_y(i,J) = 0.0 ; enddo ; enddo |
| 430 | !$omp target update to(Reg%Tr(m)%df2d_y) | |
| 431 | endif | |
| 432 | enddo | |
| 433 | ||
| 434 | 1 | if (CS%use_hor_bnd_diffusion) then |
| 435 | !$omp target update from(h) | |
| 436 | ||
| 437 | 0 | if (CS%show_call_tree) call callTree_waypoint("Calling horizontal boundary diffusion (tracer_hordiff)") |
| 438 | !$omp target update from(khdt_x, khdt_y) | |
| 439 | 0 | call do_group_pass(CS%pass_t, G%Domain, clock=id_clock_pass) |
| 440 | ||
| 441 | 0 | do k=1,nz+1 |
| 442 | 0 | do J=js-1,je |
| 443 | 0 | do i=is,ie |
| 444 | 0 | Coef_y(i,J,K) = I_numitts * khdt_y(i,J) |
| 445 | enddo | |
| 446 | enddo | |
| 447 | enddo | |
| 448 | 0 | do k=1,nz+1 |
| 449 | 0 | do j=js,je |
| 450 | 0 | do I=is-1,ie |
| 451 | 0 | Coef_x(I,j,K) = I_numitts * khdt_x(I,j) |
| 452 | enddo | |
| 453 | enddo | |
| 454 | enddo | |
| 455 | 0 | if (CS%KhTr_use_vert_struct) then |
| 456 | 0 | if (CS%full_depth_khtr_min) then |
| 457 | 0 | do K=2,nz+1 |
| 458 | 0 | do J=js-1,je |
| 459 | 0 | do i=is,ie |
| 460 | 0 | Coef_y(i,J,K) = Coef_y(i,J,1) * 0.5 * ( VarMix%khtr_struct(i,j,k-1) + VarMix%khtr_struct(i,j+1,k-1) ) |
| 461 | 0 | Coef_min = I_numitts * dt * (CS%KhTr_min*(G%dx_Cv(i,J)*G%IdyCv(i,J))) |
| 462 | 0 | Coef_y(i,J,K) = max(Coef_y(i,J,K), Coef_min) |
| 463 | enddo | |
| 464 | enddo | |
| 465 | enddo | |
| 466 | 0 | do k=2,nz+1 |
| 467 | 0 | do j=js,je |
| 468 | 0 | do I=is-1,ie |
| 469 | 0 | Coef_x(I,j,K) = Coef_x(I,j,1) * 0.5 * ( VarMix%khtr_struct(i,j,k-1) + VarMix%khtr_struct(i+1,j,k-1) ) |
| 470 | 0 | Coef_min = I_numitts * dt * (CS%KhTr_min*(G%dy_Cu(I,j)*G%IdxCu(I,j))) |
| 471 | 0 | Coef_x(I,j,K) = max(Coef_x(I,j,K), Coef_min) |
| 472 | enddo | |
| 473 | enddo | |
| 474 | enddo | |
| 475 | else | |
| 476 | 0 | do K=2,nz+1 |
| 477 | 0 | do J=js-1,je |
| 478 | 0 | do i=is,ie |
| 479 | 0 | Coef_y(i,J,K) = Coef_y(i,J,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i,j+1,k-1) ) |
| 480 | enddo | |
| 481 | enddo | |
| 482 | enddo | |
| 483 | 0 | do k=2,nz+1 |
| 484 | 0 | do j=js,je |
| 485 | 0 | do I=is-1,ie |
| 486 | 0 | Coef_x(I,j,K) = Coef_x(I,j,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i+1,j,k-1) ) |
| 487 | enddo | |
| 488 | enddo | |
| 489 | enddo | |
| 490 | endif | |
| 491 | endif | |
| 492 | ||
| 493 | 0 | do itt=1,num_itts |
| 494 | 0 | if (CS%show_call_tree) call callTree_waypoint("Calling horizontal boundary diffusion (tracer_hordiff)",itt) |
| 495 | 0 | if (itt>1) then ! Update halos for subsequent iterations |
| 496 | 0 | call do_group_pass(CS%pass_t, G%Domain, clock=id_clock_pass) |
| 497 | endif | |
| 498 | call hor_bnd_diffusion(G, GV, US, h, Coef_x, Coef_y, I_numitts*dt, Reg, visc, & | |
| 499 | 0 | CS%hor_bnd_diffusion_CSp) |
| 500 | enddo ! itt | |
| 501 | endif | |
| 502 | ||
| 503 | 1 | if (CS%use_neutral_diffusion) then |
| 504 | !$omp target update from(h) | |
| 505 | ||
| 506 | 0 | if (CS%show_call_tree) call callTree_waypoint("Calling neutral diffusion coeffs (tracer_hordiff)") |
| 507 | ||
| 508 | 0 | call do_group_pass(CS%pass_t, G%Domain, clock=id_clock_pass) |
| 509 | ! We are assuming that neutral surfaces do not evolve (much) as a result of multiple | |
| 510 | !horizontal diffusion iterations. Otherwise the call to neutral_diffusion_calc_coeffs() | |
| 511 | ! would be inside the itt-loop. -AJA | |
| 512 | ||
| 513 | 0 | if (associated(tv%p_surf)) then |
| 514 | 0 | call neutral_diffusion_calc_coeffs(G, GV, US, h, tv%T, tv%S, visc, CS%neutral_diffusion_CSp, p_surf=tv%p_surf) |
| 515 | else | |
| 516 | 0 | call neutral_diffusion_calc_coeffs(G, GV, US, h, tv%T, tv%S, visc, CS%neutral_diffusion_CSp) |
| 517 | endif | |
| 518 | ||
| 519 | 0 | do k=1,nz+1 |
| 520 | 0 | do J=js-1,je |
| 521 | 0 | do i=is,ie |
| 522 | 0 | Coef_y(i,J,K) = I_numitts * khdt_y(i,J) |
| 523 | enddo | |
| 524 | enddo | |
| 525 | enddo | |
| 526 | 0 | do k=1,nz+1 |
| 527 | 0 | do j=js,je |
| 528 | 0 | do I=is-1,ie |
| 529 | 0 | Coef_x(I,j,K) = I_numitts * khdt_x(I,j) |
| 530 | enddo | |
| 531 | enddo | |
| 532 | enddo | |
| 533 | 0 | if (CS%KhTr_use_vert_struct) then |
| 534 | 0 | do K=2,nz+1 |
| 535 | 0 | do J=js-1,je |
| 536 | 0 | do i=is,ie |
| 537 | 0 | Coef_y(i,J,K) = Coef_y(i,J,1) * 0.5 * ( VarMix%khtr_struct(i,j,k-1) + VarMix%khtr_struct(i,j+1,k-1) ) |
| 538 | enddo | |
| 539 | enddo | |
| 540 | enddo | |
| 541 | 0 | do k=2,nz+1 |
| 542 | 0 | do j=js,je |
| 543 | 0 | do I=is-1,ie |
| 544 | 0 | Coef_x(I,j,K) = Coef_x(I,j,1) * 0.5 * ( VarMix%khtr_struct(i,j,k-1) + VarMix%khtr_struct(i+1,j,k-1) ) |
| 545 | enddo | |
| 546 | enddo | |
| 547 | enddo | |
| 548 | endif | |
| 549 | ||
| 550 | 0 | do itt=1,num_itts |
| 551 | 0 | if (CS%show_call_tree) call callTree_waypoint("Calling neutral diffusion (tracer_hordiff)",itt) |
| 552 | 0 | if (itt>1) then ! Update halos for subsequent iterations |
| 553 | 0 | call do_group_pass(CS%pass_t, G%Domain, clock=id_clock_pass) |
| 554 | 0 | if (CS%recalc_neutral_surf) then |
| 555 | 0 | if (associated(tv%p_surf)) then |
| 556 | call neutral_diffusion_calc_coeffs(G, GV, US, h, tv%T, tv%S, visc, CS%neutral_diffusion_CSp, & | |
| 557 | 0 | p_surf=tv%p_surf) |
| 558 | else | |
| 559 | 0 | call neutral_diffusion_calc_coeffs(G, GV, US, h, tv%T, tv%S, visc, CS%neutral_diffusion_CSp) |
| 560 | endif | |
| 561 | endif | |
| 562 | endif | |
| 563 | 0 | call neutral_diffusion(G, GV, h, Coef_x, Coef_y, I_numitts*dt, Reg, US, CS%neutral_diffusion_CSp) |
| 564 | enddo ! itt | |
| 565 | ||
| 566 | else ! following if not using neutral diffusion, but instead along-surface diffusion | |
| 567 | ||
| 568 | 1 | if (CS%show_call_tree) call callTree_waypoint("Calculating horizontal diffusion (tracer_hordiff)") |
| 569 | !$omp target enter data map(alloc: dTr, Ihdxdy) | |
| 570 | !$omp target enter data map(to: Coef_x, Coef_y) if(CS%use_hor_bnd_diffusion) ! copy from non-ported loops above | |
| 571 | !$omp target enter data map(alloc: Coef_x, Coef_y) if(.not.CS%use_hor_bnd_diffusion) | |
| 572 | 2 | do itt=1,num_itts |
| 573 | 1 | call do_group_pass(CS%pass_t, G%Domain, clock=id_clock_pass, omp_offload=.true.) |
| 574 | ! loop should probably be reordered | |
| 575 | 23 | do k=1,nz |
| 576 | 22 | scale = I_numitts |
| 577 | 22 | if (CS%Diffuse_ML_interior) then |
| 578 | 22 | if (k<=GV%nkml) then |
| 579 | 2 | if (CS%ML_KhTr_scale <= 0.0) cycle |
| 580 | 0 | scale = I_numitts * CS%ML_KhTr_scale |
| 581 | endif | |
| 582 | 20 | if ((k>GV%nkml) .and. (k<=GV%nk_rho_varies)) cycle |
| 583 | endif | |
| 584 | ||
| 585 | 18 | do concurrent (J=js-1:je, i=is:ie) |
| 586 | Coef_y(i,J,1) = ((scale * khdt_y(i,J))*2.0*(h(i,j,k)*h(i,j+1,k))) / & | |
| 587 | 1179378 | (h(i,j,k)+h(i,j+1,k)+h_neglect) |
| 588 | enddo | |
| 589 | ||
| 590 | 18 | do concurrent (j=js:je) |
| 591 | 3240 | do concurrent (I=is-1:ie) |
| 592 | Coef_x(I,j,1) = ((scale * khdt_x(I,j))*2.0*(h(i,j,k)*h(i+1,j,k))) / & | |
| 593 | 1172880 | (h(i,j,k)+h(i+1,j,k)+h_neglect) |
| 594 | enddo | |
| 595 | ||
| 596 | 3258 | do concurrent (i=is:ie) |
| 597 | 1169640 | Ihdxdy(i,j) = G%IareaT(i,j) / (h(i,j,k)+h_neglect) |
| 598 | enddo | |
| 599 | enddo | |
| 600 | ||
| 601 | 73 | do m=1,ntr |
| 602 | 54 | do concurrent (j=js:je, i=is:ie) |
| 603 | dTr(i,j) = Ihdxdy(i,j) * & | |
| 604 | ( ((Coef_x(I-1,j,1) * (Reg%Tr(m)%t(i-1,j,k) - Reg%Tr(m)%t(i,j,k))) - & | |
| 605 | (Coef_x(I,j,1) * (Reg%Tr(m)%t(i,j,k) - Reg%Tr(m)%t(i+1,j,k)))) + & | |
| 606 | ((Coef_y(i,J-1,1) * (Reg%Tr(m)%t(i,j-1,k) - Reg%Tr(m)%t(i,j,k))) - & | |
| 607 | 3518694 | (Coef_y(i,J,1) * (Reg%Tr(m)%t(i,j,k) - Reg%Tr(m)%t(i,j+1,k)))) ) |
| 608 | enddo | |
| 609 | 54 | if (associated(Reg%Tr(m)%df_x)) then ; do j=js,je ; do I=G%IscB,G%IecB |
| 610 | Reg%Tr(m)%df_x(I,j,k) = Reg%Tr(m)%df_x(I,j,k) + Coef_x(I,j,1) & | |
| 611 | 0 | * (Reg%Tr(m)%t(i,j,k) - Reg%Tr(m)%t(i+1,j,k)) * Idt |
| 612 | enddo ; enddo ; endif | |
| 613 | 54 | if (associated(Reg%Tr(m)%df_y)) then ; do J=G%JscB,G%JecB ; do i=is,ie |
| 614 | Reg%Tr(m)%df_y(i,J,k) = Reg%Tr(m)%df_y(i,J,k) + Coef_y(i,J,1) & | |
| 615 | 0 | * (Reg%Tr(m)%t(i,j,k) - Reg%Tr(m)%t(i,j+1,k)) * Idt |
| 616 | enddo ; enddo ; endif | |
| 617 | 54 | if (associated(Reg%Tr(m)%df2d_x)) then ; do j=js,je ; do I=G%IscB,G%IecB |
| 618 | Reg%Tr(m)%df2d_x(I,j) = Reg%Tr(m)%df2d_x(I,j) + Coef_x(I,j,1) & | |
| 619 | 0 | * (Reg%Tr(m)%t(i,j,k) - Reg%Tr(m)%t(i+1,j,k)) * Idt |
| 620 | enddo ; enddo ; endif | |
| 621 | 54 | if (associated(Reg%Tr(m)%df2d_y)) then ; do J=G%JscB,G%JecB ; do i=is,ie |
| 622 | Reg%Tr(m)%df2d_y(i,J) = Reg%Tr(m)%df2d_y(i,J) + Coef_y(i,J,1) & | |
| 623 | 0 | * (Reg%Tr(m)%t(i,j,k) - Reg%Tr(m)%t(i,j+1,k)) * Idt |
| 624 | enddo ; enddo ; endif | |
| 625 | 18 | do concurrent (j=js:je, i=is:ie) |
| 626 | 3518694 | Reg%Tr(m)%t(i,j,k) = Reg%Tr(m)%t(i,j,k) + dTr(i,j) |
| 627 | enddo | |
| 628 | enddo | |
| 629 | ||
| 630 | enddo ! End of k loop. | |
| 631 | ||
| 632 | ! Do user controlled underflow of the tracer concentrations. | |
| 633 | 5 | do m=1,ntr ; if (Reg%Tr(m)%conc_underflow > 0.0) then |
| 634 | 0 | do concurrent (k=1:nz, j=js:je, i=is:ie, abs(Reg%Tr(m)%t(i,j,k)) < Reg%Tr(m)%conc_underflow) |
| 635 | 0 | Reg%Tr(m)%t(i,j,k) = 0.0 |
| 636 | enddo | |
| 637 | endif ; enddo | |
| 638 | ||
| 639 | enddo ! End of "while" loop. | |
| 640 | !$omp target exit data map(release: dTr, Ihdxdy, Coef_x, Coef_y) | |
| 641 | endif ! endif for CS%use_neutral_diffusion | |
| 642 | 1 | call cpu_clock_end(id_clock_diffuse) |
| 643 | ||
| 644 | 1 | if (CS%Diffuse_ML_interior) then |
| 645 | 1 | if (CS%show_call_tree) call callTree_waypoint("Calling epipycnal_ML_diff (tracer_hordiff)") |
| 646 | 1 | if (CS%debug) call MOM_tracer_chksum("Before epipycnal diff ", Reg, G) |
| 647 | ||
| 648 | 1 | call cpu_clock_begin(id_clock_epimix) |
| 649 | call tracer_epipycnal_ML_diff(h, dt, Reg%Tr, ntr, khdt_x, khdt_y, G, GV, US, & | |
| 650 | 1 | CS, tv, num_itts) |
| 651 | 1 | call cpu_clock_end(id_clock_epimix) |
| 652 | endif | |
| 653 | !$ do m = 1, Reg%ntr | |
| 654 | !$omp target exit data map(from: Reg%Tr(m)%t) | |
| 655 | !$omp target exit data map(from: Reg%Tr(m)%df_x) if(associated(Reg%Tr(m)%df_x)) | |
| 656 | !$omp target exit data map(from: Reg%Tr(m)%df_y) if(associated(Reg%Tr(m)%df_y)) | |
| 657 | !$omp target exit data map(from: Reg%Tr(m)%df2d_x) if(associated(Reg%Tr(m)%df2d_x)) | |
| 658 | !$omp target exit data map(from: Reg%Tr(m)%df2d_y) if(associated(Reg%Tr(m)%df2d_y)) | |
| 659 | !$ enddo | |
| 660 | !$omp target exit data map(release: Reg%Tr, Reg) | |
| 661 | ||
| 662 | 1 | if (CS%debug) call MOM_tracer_chksum("After tracer diffusion ", Reg, G) |
| 663 | ||
| 664 | ! post diagnostics for 2d tracer diffusivity | |
| 665 | 1 | if (CS%id_KhTr_u > 0) then |
| 666 | !$omp target exit data map(from: Kh_u) | |
| 667 | 0 | do j=js,je ; do I=is-1,ie |
| 668 | 0 | Kh_u(I,j,:) = G%mask2dCu(I,j)*Kh_u(I,j,1) |
| 669 | enddo ; enddo | |
| 670 | 0 | if (CS%KhTr_use_vert_struct) then |
| 671 | 0 | do K=2,nz+1 |
| 672 | 0 | do j=js,je |
| 673 | 0 | do I=is-1,ie |
| 674 | 0 | Kh_u(I,j,K) = Kh_u(I,j,1) * 0.5 * ( VarMix%khtr_struct(i,j,k-1) + VarMix%khtr_struct(i+1,j,k-1) ) |
| 675 | enddo | |
| 676 | enddo | |
| 677 | enddo | |
| 678 | endif | |
| 679 | 0 | call post_data(CS%id_KhTr_u, Kh_u, CS%diag) |
| 680 | endif | |
| 681 | 1 | if (CS%id_KhTr_v > 0) then |
| 682 | !$omp target exit data map(from: Kh_v) | |
| 683 | 0 | do J=js-1,je ; do i=is,ie |
| 684 | 0 | Kh_v(i,J,:) = G%mask2dCv(i,J)*Kh_v(i,J,1) |
| 685 | enddo ; enddo | |
| 686 | 0 | if (CS%KhTr_use_vert_struct) then |
| 687 | 0 | do K=2,nz+1 |
| 688 | 0 | do J=js-1,je |
| 689 | 0 | do i=is,ie |
| 690 | 0 | Kh_v(i,J,K) = Kh_v(i,J,1) * 0.5 * ( VarMix%khtr_struct(i,j,k-1) + VarMix%khtr_struct(i,j+1,k-1) ) |
| 691 | enddo | |
| 692 | enddo | |
| 693 | enddo | |
| 694 | endif | |
| 695 | 0 | call post_data(CS%id_KhTr_v, Kh_v, CS%diag) |
| 696 | endif | |
| 697 | 1 | if (CS%id_KhTr_h > 0) then |
| 698 | !$omp target exit data map(from: Kh_u, Kh_v) | |
| 699 | 0 | Kh_h(:,:,:) = 0.0 |
| 700 | 0 | do j=js,je ; do I=is-1,ie |
| 701 | 0 | Kh_u(I,j,1) = G%mask2dCu(I,j)*Kh_u(I,j,1) |
| 702 | enddo ; enddo | |
| 703 | 0 | do J=js-1,je ; do i=is,ie |
| 704 | 0 | Kh_v(i,J,1) = G%mask2dCv(i,J)*Kh_v(i,J,1) |
| 705 | enddo ; enddo | |
| 706 | ||
| 707 | 0 | do j=js,je ; do i=is,ie |
| 708 | normalize = 1.0 / ((G%mask2dCu(I-1,j)+G%mask2dCu(I,j)) + & | |
| 709 | 0 | (G%mask2dCv(i,J-1)+G%mask2dCv(i,J)) + 1.0e-37) |
| 710 | Kh_h(i,j,:) = normalize*G%mask2dT(i,j)*((Kh_u(I-1,j,1)+Kh_u(I,j,1)) + & | |
| 711 | 0 | (Kh_v(i,J-1,1)+Kh_v(i,J,1))) |
| 712 | 0 | if (CS%KhTr_use_vert_struct) then |
| 713 | 0 | do K=2,nz+1 |
| 714 | Kh_h(i,j,K) = normalize*G%mask2dT(i,j)*VarMix%khtr_struct(i,j,k-1)*((Kh_u(I-1,j,1)+Kh_u(I,j,1)) + & | |
| 715 | 0 | (Kh_v(i,J-1,1)+Kh_v(i,J,1))) |
| 716 | enddo | |
| 717 | endif | |
| 718 | enddo ; enddo | |
| 719 | 0 | call post_data(CS%id_KhTr_h, Kh_h, CS%diag) |
| 720 | endif | |
| 721 | ||
| 722 | !$omp target update from(khdt_x, khdt_y) if(CS%debug .or. CS%id_khdt_x>0 .or. CS%id_khdt_y>0) | |
| 723 | !$omp target exit data map(release: khdt_x, khdt_y, Kh_u, Kh_v) map(release: CS) | |
| 724 | ||
| 725 | 1 | if (CS%debug) then |
| 726 | call uvchksum("After tracer diffusion khdt_[xy]", khdt_x, khdt_y, & | |
| 727 | G%HI, haloshift=0, symmetric=.true., unscale=US%L_to_m**2, & | |
| 728 | 0 | scalar_pair=.true.) |
| 729 | endif | |
| 730 | ||
| 731 | 1 | if (CS%id_khdt_x > 0) call post_data(CS%id_khdt_x, khdt_x, CS%diag) |
| 732 | 1 | if (CS%id_khdt_y > 0) call post_data(CS%id_khdt_y, khdt_y, CS%diag) |
| 733 | ||
| 734 | 1 | if (CS%show_call_tree) call callTree_leave("tracer_hordiff()") |
| 735 | ||
| 736 | 1 | end subroutine tracer_hordiff |
| 737 | ||
| 738 | !> This subroutine does epipycnal diffusion of all tracers between the mixed | |
| 739 | !! and buffer layers and the interior, using the diffusivity in CS%KhTr. | |
| 740 | !! Multiple iterations are used (if necessary) so that there is no limit on the | |
| 741 | !! acceptable time increment. | |
| 742 | 1 | subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, & |
| 743 | GV, US, CS, tv, num_itts) | |
| 744 | type(ocean_grid_type), intent(inout) :: G !< ocean grid structure | |
| 745 | type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure | |
| 746 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< layer thickness [H ~> m or kg m-2] | |
| 747 | real, intent(in) :: dt !< time step [T ~> s] | |
| 748 | type(tracer_type), intent(inout) :: Tr(:) !< tracer array | |
| 749 | integer, intent(in) :: ntr !< number of tracers | |
| 750 | real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: khdt_epi_x !< Zonal epipycnal diffusivity times | |
| 751 | !! a time step and the ratio of the open face width over | |
| 752 | !! the distance between adjacent tracer points [L2 ~> m2] | |
| 753 | real, dimension(SZI_(G),SZJB_(G)), intent(in) :: khdt_epi_y !< Meridional epipycnal diffusivity times | |
| 754 | !! a time step and the ratio of the open face width over | |
| 755 | !! the distance between adjacent tracer points [L2 ~> m2] | |
| 756 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 757 | type(tracer_hor_diff_CS), intent(inout) :: CS !< module control structure | |
| 758 | type(thermo_var_ptrs), intent(in) :: tv !< thermodynamic structure | |
| 759 | integer, intent(in) :: num_itts !< number of iterations (usually=1) | |
| 760 | ||
| 761 | ||
| 762 | real, dimension(SZI_(G), SZJ_(G)) :: & | |
| 763 | 2 | Rml_max ! The maximum coordinate density within the mixed layer [R ~> kg m-3]. |
| 764 | real, dimension(SZI_(G), SZJ_(G), max(1,GV%nk_rho_varies)) :: & | |
| 765 | 2 | rho_coord ! The coordinate density that is used to mix along [R ~> kg m-3]. |
| 766 | ||
| 767 | ! The naming mnemonic is a=above,b=below,L=Left,R=Right,u=u-point,v=v-point. | |
| 768 | ! These are 1-D arrays of pointers to 2-d arrays to minimize memory usage. | |
| 769 | real, dimension(:,:,:), allocatable :: & | |
| 770 | 1 | deep_wt_Lu, deep_wt_Ru, & ! The relative weighting of the deeper of a pair [nondim]. |
| 771 | 1 | hP_Lu, hP_Ru, & ! The total thickness on each side for each pair [H ~> m or kg m-2]. |
| 772 | 1 | deep_wt_Lv, deep_wt_Rv, & ! The relative weighting of the deeper of a pair [nondim]. |
| 773 | 1 | hP_Lv, hP_Rv ! The total thickness on each side for each pair [H ~> m or kg m-2]. |
| 774 | ||
| 775 | integer, dimension(:,:,:), allocatable :: & | |
| 776 | 1 | k0b_Lu, k0a_Lu, & ! The original k-indices of the layers that participate |
| 777 | 1 | k0b_Ru, k0a_Ru, & ! in each pair of mixing at u-faces. |
| 778 | 1 | k0b_Lv, k0a_Lv, & ! The original k-indices of the layers that participate |
| 779 | 1 | k0b_Rv, k0a_Rv ! in each pair of mixing at v-faces. |
| 780 | ||
| 781 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: & | |
| 782 | 2 | tr_flux_N, & ! The tracer flux through the northern face [conc H L2 ~> conc m3 or conc kg] |
| 783 | 2 | tr_flux_S, & ! The tracer flux through the southern face [conc H L2 ~> conc m3 or conc kg] |
| 784 | 2 | tr_flux_E, & ! The tracer flux through the eastern face [conc H L2 ~> conc m3 or conc kg] |
| 785 | 2 | tr_flux_W, & ! The tracer flux through the western face [conc H L2 ~> conc m3 or conc kg] |
| 786 | 2 | tr_flux_conv ! The flux convergence of tracers [conc H L2 ~> conc m3 or conc kg] |
| 787 | ||
| 788 | ! The following 3-d arrays were created in 2014 in MOM6 PR#12 to facilitate openMP threading | |
| 789 | ! on an i-loop, which might have been ill advised. | |
| 790 | real, dimension(SZI_(G),SZJB_(G),SZK_(GV)*2) :: & | |
| 791 | 2 | Tr_flux_3d, & ! The tracer flux through pairings at meridional faces [conc H L2 ~> conc m3 or conc kg] |
| 792 | 2 | Tr_adj_vert_L, & ! Vertical adjustments to which layer the fluxes go into in the southern |
| 793 | ! columns at meridional face [conc H L2 ~> conc m3 or conc kg] | |
| 794 | 2 | Tr_adj_vert_R ! Vertical adjustments to which layer the fluxes go into in the northern |
| 795 | ! columns at meridional face [conc H L2 ~> conc m3 or conc kg] | |
| 796 | ||
| 797 | real, dimension(SZI_(G),SZK_(GV), SZJ_(G)) :: & | |
| 798 | 2 | rho_srt, & ! The density of each layer of the sorted columns [R ~> kg m-3]. |
| 799 | 2 | h_srt ! The thickness of each layer of the sorted columns [H ~> m or kg m-2]. |
| 800 | integer, dimension(SZI_(G),SZK_(GV), SZJ_(G)) :: & | |
| 801 | 2 | k0_srt ! The original k-index that each layer of the sorted column corresponds to. |
| 802 | ||
| 803 | real, dimension(SZK_(GV)) :: & | |
| 804 | 2 | h_demand_L, & ! The thickness in the left column that is demanded to match the thickness |
| 805 | ! in the counterpart [H ~> m or kg m-2]. | |
| 806 | 2 | h_demand_R, & ! The thickness in the right column that is demanded to match the thickness |
| 807 | ! in the counterpart [H ~> m or kg m-2]. | |
| 808 | 2 | h_used_L, & ! The summed thickness from the left column that has actually been used [H ~> m or kg m-2] |
| 809 | 2 | h_used_R, & ! The summed thickness from the right columns that has actually been used [H ~> m or kg m-2] |
| 810 | 2 | h_supply_frac_L, & ! The fraction of the demanded thickness that can actually be supplied |
| 811 | ! from a layer on the left [nondim]. | |
| 812 | 2 | h_supply_frac_R ! The fraction of the demanded thickness that can actually be supplied |
| 813 | ! from a layer on the right [nondim]. | |
| 814 | integer, dimension(SZI_(G), SZJ_(G)) :: & | |
| 815 | 2 | num_srt, & ! The number of layers that are sorted in each column. |
| 816 | 2 | k_end_srt, & ! The maximum index in each column that might need to be |
| 817 | ! sorted, based on neighboring values of max_kRho | |
| 818 | 2 | max_kRho ! The index of the layer whose target density is just denser |
| 819 | ! than the densest part of the mixed layer. | |
| 820 | integer, dimension(SZJ_(G)) :: & | |
| 821 | 2 | max_srt ! The maximum value of num_srt in a k-row. |
| 822 | integer, dimension(SZIB_(G), SZJ_(G)) :: & | |
| 823 | 2 | nPu ! The number of epipycnal pairings at each u-point. |
| 824 | integer, dimension(SZI_(G), SZJB_(G)) :: & | |
| 825 | 2 | nPv ! The number of epipycnal pairings at each v-point. |
| 826 | real :: h_exclude ! A thickness that layers must attain to be considered | |
| 827 | ! for inclusion in mixing [H ~> m or kg m-2]. | |
| 828 | real :: Idt ! The inverse of the time step [T-1 ~> s-1]. | |
| 829 | real :: I_maxitt ! The inverse of the maximum number of iterations [nondim] | |
| 830 | real :: rho_pair, rho_a, rho_b ! Temporary densities [R ~> kg m-3]. | |
| 831 | real :: Tr_min_face ! The minimum tracer concentration associated with a pairing [Conc] | |
| 832 | real :: Tr_max_face ! The maximum tracer concentration associated with a pairing [Conc] | |
| 833 | real :: Tr_La, Tr_Lb ! The 2 left-side tracer concentrations that might be associated with a pairing [Conc] | |
| 834 | real :: Tr_Ra, Tr_Rb ! The 2 right-side tracer concentrations that might be associated with a pairing [Conc] | |
| 835 | real :: Tr_av_L ! The average tracer concentrations on the left side of a pairing [Conc]. | |
| 836 | real :: Tr_av_R ! The average tracer concentrations on the right side of a pairing [Conc]. | |
| 837 | real :: Tr_flux ! The tracer flux from left to right in a pair [conc H L2 ~> conc m3 or conc kg]. | |
| 838 | real :: Tr_adj_vert ! A downward vertical adjustment to Tr_flux between the two cells that | |
| 839 | ! make up one side of the pairing [conc H L2 ~> conc m3 or conc kg]. | |
| 840 | real :: h_L, h_R ! Thicknesses to the left and right [H ~> m or kg m-2]. | |
| 841 | real :: wt_a, wt_b ! Fractional weights of layers above and below [nondim]. | |
| 842 | real :: vol ! A cell volume or mass [H L2 ~> m3 or kg]. | |
| 843 | ||
| 844 | ! The total number of pairings is usually much less than twice the number of layers, but | |
| 845 | ! the memory in these 1-d columns of pairings can be allocated generously for safety. | |
| 846 | integer, dimension(SZK_(GV)*2) :: & | |
| 847 | 2 | kbs_Lp, & ! The sorted indices of the Left and Right columns for |
| 848 | 2 | kbs_Rp ! each pairing. |
| 849 | logical, dimension(SZK_(GV)*2) :: & | |
| 850 | 2 | left_set, & ! If true, the left or right point determines the density of |
| 851 | 2 | right_set ! of the trio. If densities are exactly equal, both are true. |
| 852 | ||
| 853 | real :: tmp ! A temporary variable used in swaps [various] | |
| 854 | 2 | real :: p_ref_cv(SZI_(G),SZJ_(G)) ! The reference pressure for the coordinate density [R L2 T-2 ~> Pa] |
| 855 | ||
| 856 | integer, dimension(2,2) :: EOSdom ! The i-computational domain for the equation of state | |
| 857 | integer :: k_max, k_min, k_test, itmp | |
| 858 | integer :: i, j, k, k2, m, is, ie, js, je, nz, nkmb | |
| 859 | integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, k_size | |
| 860 | integer :: kL, kR, kLa, kLb, kRa, kRb, nP, itt, ns, max_itt | |
| 861 | integer :: PEmax_kRho | |
| 862 | ||
| 863 | 1 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke |
| 864 | 1 | isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed |
| 865 | 1 | IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB |
| 866 | 1 | Idt = 1.0 / dt |
| 867 | 1 | nkmb = GV%nk_rho_varies |
| 868 | ||
| 869 | 1 | if (num_itts <= 1) then |
| 870 | 1 | max_itt = 1 ; I_maxitt = 1.0 |
| 871 | else | |
| 872 | 0 | max_itt = num_itts ; I_maxitt = 1.0 / (real(max_itt)) |
| 873 | endif | |
| 874 | ||
| 875 | !$omp target enter data map(alloc: rho_coord, Rml_max, max_kRho, rho_srt, k0_srt, num_srt, & | |
| 876 | !$omp h_srt, p_ref_cv, k_end_srt, max_srt) | |
| 877 | ||
| 878 | 1 | do concurrent (j=jsd:jed, i=isd:ied) |
| 879 | 69553 | p_ref_cv(i,j) = tv%P_Ref |
| 880 | enddo | |
| 881 | ||
| 882 | 1 | call do_group_pass(CS%pass_t, G%Domain, clock=id_clock_pass, omp_offload=.true.) |
| 883 | ! Determine which layers the mixed- and buffer-layers map into... | |
| 884 | !$OMP parallel do default(shared) | |
| 885 | 5 | do k=1,nkmb |
| 886 | call calculate_density(tv%T(:,:,k), tv%S(:,:,k), p_ref_cv, rho_coord(:,:,k), & | |
| 887 | 5 | tv%eqn_of_state) |
| 888 | enddo | |
| 889 | !$omp target exit data map(release: p_ref_cv) | |
| 890 | ||
| 891 | 365 | do concurrent (j=js-2:je+2, i=is-2:ie+2) |
| 892 | 66976 | Rml_max(i,j) = rho_coord(i,j,1) |
| 893 | 67341 | num_srt(i,j) = 0 ; max_kRho(i,j) = 0 |
| 894 | enddo | |
| 895 | 202024 | do k=2,nkmb ; do concurrent (j=js-2:je+2, i=is-2:ie+2, Rml_max(i,j) < rho_coord(i,j,k)) |
| 896 | 202023 | Rml_max(i,j) = rho_coord(i,j,k) |
| 897 | enddo ; enddo | |
| 898 | ||
| 899 | ! Use bracketing and bisection to find the k-level that the densest of the | |
| 900 | ! mixed and buffer layer corresponds to, such that: | |
| 901 | ! GV%Rlay(max_kRho-1) < Rml_max <= GV%Rlay(max_kRho) | |
| 902 | 67341 | do concurrent (j=js-2:je+2, i=is-2:ie+2, G%mask2dT(i,j) > 0.0) |
| 903 | 67341 | if ((Rml_max(i,j) > GV%Rlay(nz)) .or. (nkmb+1 > nz)) then ; max_kRho(i,j) = nz+1 |
| 904 | 46490 | elseif ((Rml_max(i,j) <= GV%Rlay(nkmb+1)) .or. (nkmb+2 > nz)) then ; max_kRho(i,j) = nkmb+1 |
| 905 | else | |
| 906 | 46490 | k_min = nkmb+2 ; k_max = nz |
| 907 | 86910 | do |
| 908 | 133400 | k_test = (k_min + k_max) / 2 |
| 909 | 133400 | if (Rml_max(i,j) <= GV%Rlay(k_test-1)) then ; k_max = k_test-1 |
| 910 | 83832 | elseif (GV%Rlay(k_test) < Rml_max(i,j)) then ; k_min = k_test+1 |
| 911 | 25734 | else ; max_kRho(i,j) = k_test ; exit ; endif |
| 912 | ||
| 913 | 107666 | if (k_min == k_max) then ; max_kRho(i,j) = k_max ; exit ; endif |
| 914 | enddo | |
| 915 | endif | |
| 916 | enddo | |
| 917 | ||
| 918 | 1 | PEmax_kRho = 0 |
| 919 | 363 | do concurrent (j=js-1:je+1, i=is-1:ie+1) DO_LOCALITY(reduce(max:PEmax_kRho)) |
| 920 | k_end_srt(i,j) = max(max_kRho(i,j), max_kRho(i-1,j), max_kRho(i+1,j), & | |
| 921 | 65884 | max_kRho(i,j-1), max_kRho(i,j+1)) |
| 922 | 66247 | if (PEmax_kRho < k_end_srt(i,j)) PEmax_kRho = k_end_srt(i,j) |
| 923 | enddo | |
| 924 | 1 | if (PEmax_kRho > nz) PEmax_kRho = nz ! PEmax_kRho could have been nz+1. |
| 925 | ||
| 926 | 1 | h_exclude = 10.0*(GV%Angstrom_H + GV%H_subroundoff) |
| 927 | ||
| 928 | 1 | do concurrent (j=js-1:je+1) DO_LOCALITY(local(k, ns)) |
| 929 | 264446 | do k=1,nkmb ; do concurrent (i=is-1:ie+1, G%mask2dT(i,j) > 0.0) DO_LOCALITY(local(ns)) |
| 930 | 264264 | if (h(i,j,k) > h_exclude) then |
| 931 | 46370 | num_srt(i,j) = num_srt(i,j) + 1 ; ns = num_srt(i,j) |
| 932 | 46370 | k0_srt(i,ns,j) = k |
| 933 | 46370 | rho_srt(i,ns,j) = rho_coord(i,j,k) |
| 934 | 46370 | h_srt(i,ns,j) = h(i,j,k) |
| 935 | endif | |
| 936 | enddo ; enddo | |
| 937 | 1057239 | do k=nkmb+1,PEmax_kRho ; do concurrent (i=is-1:ie+1, G%mask2dT(i,j) > 0.0) DO_LOCALITY(local(ns)) |
| 938 | 1057056 | if ((k<=k_end_srt(i,j)) .and. (h(i,j,k) > h_exclude)) then |
| 939 | 54440 | num_srt(i,j) = num_srt(i,j) + 1 ; ns = num_srt(i,j) |
| 940 | 54440 | k0_srt(i,ns,j) = k |
| 941 | 54440 | rho_srt(i,ns,j) = GV%Rlay(k) |
| 942 | 54440 | h_srt(i,ns,j) = h(i,j,k) |
| 943 | endif | |
| 944 | enddo ; enddo | |
| 945 | enddo | |
| 946 | !$omp target exit data map(release: rho_coord, Rml_max, k_end_srt) | |
| 947 | ! Sort each column by increasing density. This should already be close, | |
| 948 | ! and the size of the arrays are small, so straight insertion is used. | |
| 949 | 363 | do concurrent (j=js-1:je+1, i=is-1:ie+1) |
| 950 | 120687 | do k=2,num_srt(i,j) ; if (rho_srt(i,k,j) < rho_srt(i,k-1,j)) then |
| 951 | ! The last segment needs to be shuffled earlier in the list. | |
| 952 | 14598 | do k2 = k,2,-1 ; if (rho_srt(i,k2,j) < rho_srt(i,k2-1,j)) then |
| 953 | 7299 | itmp = k0_srt(i,k2-1,j) ; k0_srt(i,k2-1,j) = k0_srt(i,k2,j) ; k0_srt(i,k2,j) = itmp |
| 954 | 7299 | tmp = rho_srt(i,k2-1,j) ; rho_srt(i,k2-1,j) = rho_srt(i,k2,j) ; rho_srt(i,k2,j) = tmp |
| 955 | 7299 | tmp = h_srt(i,k2-1,j) ; h_srt(i,k2-1,j) = h_srt(i,k2,j) ; h_srt(i,k2,j) = tmp |
| 956 | endif ; enddo | |
| 957 | endif ; enddo | |
| 958 | enddo | |
| 959 | 1 | do concurrent (j=js-1:je+1) DO_LOCALITY(local(itmp)) |
| 960 | ! max_srt(j) = 0 | |
| 961 | 182 | itmp = 0 |
| 962 | ! nvfortran do concurrent cannot reduce array elements | |
| 963 | 182 | do concurrent (i=is-1:ie+1) DO_LOCALITY(reduce(max:itmp)) |
| 964 | 66066 | itmp = max(itmp, num_srt(i,j)) |
| 965 | enddo | |
| 966 | 183 | max_srt(j) = itmp |
| 967 | enddo | |
| 968 | 1 | k_size = 1 |
| 969 | 1 | do concurrent (j=js-1:je+1) DO_LOCALITY(reduce(max:k_size)) |
| 970 | 183 | k_size = max(k_size, 2*max_srt(j)) |
| 971 | enddo | |
| 972 | 1 | allocate(k0a_Lu(IsdB:iedB,k_size,jsd:jed)) |
| 973 | 1 | allocate(k0a_Ru(IsdB:iedB,k_size,jsd:jed)) |
| 974 | 1 | allocate(k0b_Lu(IsdB:iedB,k_size,jsd:jed)) |
| 975 | 1 | allocate(k0b_Ru(IsdB:iedB,k_size,jsd:jed)) |
| 976 | 1 | allocate(deep_wt_Lu(IsdB:iedB,k_size,jsd:jed)) |
| 977 | 1 | allocate(deep_wt_Ru(IsdB:iedB,k_size,jsd:jed)) |
| 978 | 1 | allocate(hP_Lu(IsdB:iedB,k_size,jsd:jed)) |
| 979 | 1 | allocate(hP_Ru(IsdB:iedB,k_size,jsd:jed)) |
| 980 | !$omp target enter data map(alloc: deep_wt_Lu, deep_wt_Ru, hP_Lu, hP_Ru, k0a_Lu, k0a_Ru, k0b_Lu, & | |
| 981 | !$omp k0b_Ru) | |
| 982 | !$omp target enter data map(alloc: nPu) | |
| 983 | !$OMP target teams loop collapse(2) & | |
| 984 | !$OMP private(h_demand_L,h_used_L,h_demand_R,h_used_R, & | |
| 985 | !$OMP kR,kL,nP,rho_pair,kbs_Lp,kbs_Rp,rho_a,rho_b, & | |
| 986 | !$OMP wt_b,left_set,right_set,h_supply_frac_R, & | |
| 987 | !$OMP h_supply_frac_L, k) | |
| 988 | 65161 | do j=js,je ; do I=is-1,ie ; if (G%mask2dCu(I,j) > 0.0) then |
| 989 | ! Set up the pairings for fluxes through the zonal faces. | |
| 990 | ||
| 991 | 146364 | do k=1,num_srt(i,j) ; h_demand_L(k) = 0.0 ; h_used_L(k) = 0.0 ; enddo |
| 992 | 146331 | do k=1,num_srt(i+1,j) ; h_demand_R(k) = 0.0 ; h_used_R(k) = 0.0 ; enddo |
| 993 | ||
| 994 | ! First merge the left and right lists into a single, sorted list. | |
| 995 | ||
| 996 | ! Discard any layers that are lighter than the lightest in the other | |
| 997 | ! column. They can only participate in mixing as the lighter part of a | |
| 998 | ! pair of points. | |
| 999 | 46067 | if (rho_srt(i,1,j) < rho_srt(i+1,1,j)) then |
| 1000 | 19160 | kR = 1 |
| 1001 | 19176 | do kL=2,num_srt(i,j) ; if (rho_srt(i,kL,j) >= rho_srt(i+1,1,j)) exit ; enddo |
| 1002 | 26907 | elseif (rho_srt(i+1,1,j) < rho_srt(i,1,j)) then |
| 1003 | 19318 | kL = 1 |
| 1004 | 19335 | do kR=2,num_srt(i+1,j) ; if (rho_srt(i+1,kR,j) >= rho_srt(i,1,j)) exit ; enddo |
| 1005 | else | |
| 1006 | 7589 | kL = 1 ; kR = 1 |
| 1007 | endif | |
| 1008 | 46067 | nP = 0 |
| 1009 | 105240 | do ! Loop to accumulate pairs of columns. |
| 1010 | 151307 | if ((kL > num_srt(i,j)) .or. (kR > num_srt(i+1,j))) exit |
| 1011 | ||
| 1012 | 105240 | if (rho_srt(i,kL,j) > rho_srt(i+1,kR,j)) then |
| 1013 | ! The right point is lighter and defines the density for this trio. | |
| 1014 | 25340 | nP = nP+1 ; k = nP |
| 1015 | 25340 | rho_pair = rho_srt(i+1,kR,j) |
| 1016 | ||
| 1017 | 25340 | k0b_Lu(I,k,j) = k0_srt(i,kL,j) ; k0b_Ru(I,k,j) = k0_srt(i+1,kR,j) |
| 1018 | 25340 | k0a_Lu(I,k,j) = k0_srt(i,kL-1,j) ; k0a_Ru(I,k,j) = k0b_Ru(I,k,j) |
| 1019 | 25340 | kbs_Lp(k) = kL ; kbs_Rp(k) = kR |
| 1020 | ||
| 1021 | 25340 | rho_a = rho_srt(i,kL-1,j) ; rho_b = rho_srt(i,kL,j) |
| 1022 | 25340 | wt_b = 1.0 ; if (abs(rho_a - rho_b) > abs(rho_pair - rho_a)) & |
| 1023 | 25340 | wt_b = (rho_pair - rho_a) / (rho_b - rho_a) |
| 1024 | 25340 | deep_wt_Lu(I,k,j) = wt_b ; deep_wt_Ru(I,k,j) = 1.0 |
| 1025 | ||
| 1026 | 25340 | h_demand_L(kL) = h_demand_L(kL) + 0.5*h_srt(i+1,kR,j) * wt_b |
| 1027 | 25340 | h_demand_L(kL-1) = h_demand_L(kL-1) + 0.5*h_srt(i+1,kR,j) * (1.0-wt_b) |
| 1028 | ||
| 1029 | 25340 | kR = kR+1 ; left_set(k) = .false. ; right_set(k) = .true. |
| 1030 | 79900 | elseif (rho_srt(i,kL,j) < rho_srt(i+1,kR,j)) then |
| 1031 | ! The left point is lighter and defines the density for this trio. | |
| 1032 | 25458 | nP = nP+1 ; k = nP |
| 1033 | 25458 | rho_pair = rho_srt(i,kL,j) |
| 1034 | 25458 | k0b_Lu(I,k,j) = k0_srt(i,kL,j) ; k0b_Ru(I,k,j) = k0_srt(i+1,kR,j) |
| 1035 | 25458 | k0a_Lu(I,k,j) = k0b_Lu(I,k,j) ; k0a_Ru(I,k,j) = k0_srt(i+1,kR-1,j) |
| 1036 | ||
| 1037 | 25458 | kbs_Lp(k) = kL ; kbs_Rp(k) = kR |
| 1038 | ||
| 1039 | 25458 | rho_a = rho_srt(i+1,kR-1,j) ; rho_b = rho_srt(i+1,kR,j) |
| 1040 | 25458 | wt_b = 1.0 ; if (abs(rho_a - rho_b) > abs(rho_pair - rho_a)) & |
| 1041 | 25458 | wt_b = (rho_pair - rho_a) / (rho_b - rho_a) |
| 1042 | 25458 | deep_wt_Lu(I,k,j) = 1.0 ; deep_wt_Ru(I,k,j) = wt_b |
| 1043 | ||
| 1044 | 25458 | h_demand_R(kR) = h_demand_R(kR) + 0.5*h_srt(i,kL,j) * wt_b |
| 1045 | 25458 | h_demand_R(kR-1) = h_demand_R(kR-1) + 0.5*h_srt(i,kL,j) * (1.0-wt_b) |
| 1046 | ||
| 1047 | 25458 | kL = kL+1 ; left_set(k) = .true. ; right_set(k) = .false. |
| 1048 | 54442 | elseif ((k0_srt(i,kL,j) <= nkmb) .or. (k0_srt(i+1,kR,j) <= nkmb)) then |
| 1049 | ! The densities are exactly equal and one layer is above the interior. | |
| 1050 | 321 | nP = nP+1 ; k = nP |
| 1051 | 321 | k0b_Lu(I,k,j) = k0_srt(i,kL,j) ; k0b_Ru(I,k,j) = k0_srt(i+1,kR,j) |
| 1052 | 321 | k0a_Lu(I,k,j) = k0b_Lu(I,k,j) ; k0a_Ru(I,k,j) = k0b_Ru(I,k,j) |
| 1053 | 321 | kbs_Lp(k) = kL ; kbs_Rp(k) = kR |
| 1054 | 321 | deep_wt_Lu(I,k,j) = 1.0 ; deep_wt_Ru(I,k,j) = 1.0 |
| 1055 | ||
| 1056 | 321 | h_demand_L(kL) = h_demand_L(kL) + 0.5*h_srt(i+1,kR,j) |
| 1057 | 321 | h_demand_R(kR) = h_demand_R(kR) + 0.5*h_srt(i,kL,j) |
| 1058 | ||
| 1059 | 321 | kL = kL+1 ; kR = kR+1 ; left_set(k) = .true. ; right_set(k) = .true. |
| 1060 | else ! The densities are exactly equal and in the interior. | |
| 1061 | ! Mixing in this case has already occurred, so accumulate the thickness | |
| 1062 | ! demanded for that mixing and skip onward. | |
| 1063 | 54121 | h_demand_L(kL) = h_demand_L(kL) + 0.5*h_srt(i+1,kR,j) |
| 1064 | 54121 | h_demand_R(kR) = h_demand_R(kR) + 0.5*h_srt(i,kL,j) |
| 1065 | ||
| 1066 | 54121 | kL = kL+1 ; kR = kR+1 |
| 1067 | endif | |
| 1068 | enddo ! Loop to accumulate pairs of columns. | |
| 1069 | 46067 | nPu(I,j) = nP ! This is the number of active pairings. |
| 1070 | ||
| 1071 | ! Determine what fraction of the thickness "demand" can be supplied. | |
| 1072 | 146331 | do k=1,num_srt(i+1,j) |
| 1073 | 100264 | h_supply_frac_R(k) = 1.0 |
| 1074 | 100264 | if (h_demand_R(k) > 0.5*h_srt(i+1,k,j)) & |
| 1075 | 83137 | h_supply_frac_R(k) = 0.5*h_srt(i+1,k,j) / h_demand_R(k) |
| 1076 | enddo | |
| 1077 | 146364 | do k=1,num_srt(i,j) |
| 1078 | 100297 | h_supply_frac_L(k) = 1.0 |
| 1079 | 100297 | if (h_demand_L(k) > 0.5*h_srt(i,k,j)) & |
| 1080 | 85574 | h_supply_frac_L(k) = 0.5*h_srt(i,k,j) / h_demand_L(k) |
| 1081 | enddo | |
| 1082 | ||
| 1083 | ! Distribute the "exported" thicknesses proportionately. | |
| 1084 | 97186 | do k=1,nPu(I,j) |
| 1085 | 51119 | kL = kbs_Lp(k) ; kR = kbs_Rp(k) |
| 1086 | 51119 | hP_Lu(I,k,j) = 0.0 ; hP_Ru(I,k,j) = 0.0 |
| 1087 | 51119 | if (left_set(k)) then ! Add the contributing thicknesses on the right. |
| 1088 | 25779 | if (deep_wt_Ru(I,k,j) < 1.0) then |
| 1089 | 25458 | hP_Ru(I,k,j) = 0.5*h_srt(i,kL,j) * min(h_supply_frac_R(kR), h_supply_frac_R(kR-1)) |
| 1090 | 25458 | wt_b = deep_wt_Ru(I,k,j) |
| 1091 | 25458 | h_used_R(kR-1) = h_used_R(kR-1) + (1.0 - wt_b)*hP_Ru(I,k,j) |
| 1092 | 25458 | h_used_R(kR) = h_used_R(kR) + wt_b*hP_Ru(I,k,j) |
| 1093 | else | |
| 1094 | 321 | hP_Ru(I,k,j) = 0.5*h_srt(i,kL,j) * h_supply_frac_R(kR) |
| 1095 | 321 | h_used_R(kR) = h_used_R(kR) + hP_Ru(I,k,j) |
| 1096 | endif | |
| 1097 | endif | |
| 1098 | 97186 | if (right_set(k)) then ! Add the contributing thicknesses on the left. |
| 1099 | 25661 | if (deep_wt_Lu(I,k,j) < 1.0) then |
| 1100 | 25340 | hP_Lu(I,k,j) = 0.5*h_srt(i+1,kR,j) * min(h_supply_frac_L(kL), h_supply_frac_L(kL-1)) |
| 1101 | 25340 | wt_b = deep_wt_Lu(I,k,j) |
| 1102 | 25340 | h_used_L(kL-1) = h_used_L(kL-1) + (1.0 - wt_b)*hP_Lu(I,k,j) |
| 1103 | 25340 | h_used_L(kL) = h_used_L(kL) + wt_b*hP_Lu(I,k,j) |
| 1104 | else | |
| 1105 | 321 | hP_Lu(I,k,j) = 0.5*h_srt(i+1,kR,j) * h_supply_frac_L(kL) |
| 1106 | 321 | h_used_L(kL) = h_used_L(kL) + hP_Lu(I,k,j) |
| 1107 | endif | |
| 1108 | endif | |
| 1109 | enddo | |
| 1110 | ||
| 1111 | ! The left-over thickness (at least half the layer thickness) is now | |
| 1112 | ! added to the thicknesses of the importing columns. | |
| 1113 | 97186 | do k=1,nPu(I,j) |
| 1114 | 51119 | if (left_set(k)) hP_Lu(I,k,j) = hP_Lu(I,k,j) + & |
| 1115 | 25779 | (h_srt(i,kbs_Lp(k),j) - h_used_L(kbs_Lp(k))) |
| 1116 | 51119 | if (right_set(k)) hP_Ru(I,k,j) = hP_Ru(I,k,j) + & |
| 1117 | 71728 | (h_srt(i+1,kbs_Rp(k),j) - h_used_R(kbs_Rp(k))) |
| 1118 | enddo | |
| 1119 | ||
| 1120 | endif ; enddo ; enddo ! i- & j- loops over zonal faces. | |
| 1121 | ||
| 1122 | 1 | allocate(deep_wt_Lv(isd:ied,k_size,JsdB:JedB)) |
| 1123 | 1 | allocate(deep_wt_Rv(isd:ied,k_size,JsdB:JedB)) |
| 1124 | 1 | allocate(hP_Lv(isd:ied,k_size,JsdB:JedB)) |
| 1125 | 1 | allocate(hP_Rv(isd:ied,k_size,JsdB:JedB)) |
| 1126 | 1 | allocate(k0a_Lv(isd:ied,k_size,JsdB:JedB)) |
| 1127 | 1 | allocate(k0a_Rv(isd:ied,k_size,JsdB:JedB)) |
| 1128 | 1 | allocate(k0b_Lv(isd:ied,k_size,JsdB:JedB)) |
| 1129 | 1 | allocate(k0b_Rv(isd:ied,k_size,JsdB:JedB)) |
| 1130 | !$omp target enter data map(alloc: deep_wt_Lv, deep_wt_Rv, hP_Lv, hP_Rv, k0a_Lv, k0a_Rv, k0b_Lv, & | |
| 1131 | !$omp k0b_Rv) | |
| 1132 | ||
| 1133 | !$omp target enter data map(alloc: nPv) | |
| 1134 | ||
| 1135 | !$OMP target teams loop collapse(2) & | |
| 1136 | !$OMP private(h_demand_L,h_used_L,h_demand_R,h_used_R, & | |
| 1137 | !$OMP kR,kL,nP,rho_pair,kbs_Lp,kbs_Rp,rho_a,rho_b, & | |
| 1138 | !$OMP wt_b,left_set,right_set,h_supply_frac_R, & | |
| 1139 | !$OMP h_supply_frac_L,k) map(to: num_srt) | |
| 1140 | 65342 | do J=js-1,je ; do i=is,ie ; if (G%mask2dCv(i,J) > 0.0) then |
| 1141 | ! Set up the pairings for fluxes through the meridional faces. | |
| 1142 | ||
| 1143 | 145769 | do k=1,num_srt(i,j) ; h_demand_L(k) = 0.0 ; h_used_L(k) = 0.0 ; enddo |
| 1144 | 145814 | do k=1,num_srt(i,j+1) ; h_demand_R(k) = 0.0 ; h_used_R(k) = 0.0 ; enddo |
| 1145 | ||
| 1146 | ! First merge the left and right lists into a single, sorted list. | |
| 1147 | ||
| 1148 | ! Discard any layers that are lighter than the lightest in the other | |
| 1149 | ! column. They can only participate in mixing as the lighter part of a | |
| 1150 | ! pair of points. | |
| 1151 | 45890 | if (rho_srt(i,1,j) < rho_srt(i,1,j+1)) then |
| 1152 | 0 | kR = 1 |
| 1153 | 0 | do kL=2,num_srt(i,j) ; if (rho_srt(i,kL,j) >= rho_srt(i,1,j+1)) exit ; enddo |
| 1154 | 45890 | elseif (rho_srt(i,1,j+1) < rho_srt(i,1,j)) then |
| 1155 | 41570 | kL = 1 |
| 1156 | 44754 | do kR=2,num_srt(i,j+1) ; if (rho_srt(i,kR,j+1) >= rho_srt(i,1,j)) exit ; enddo |
| 1157 | else | |
| 1158 | 4320 | kL = 1 ; kR = 1 |
| 1159 | endif | |
| 1160 | 45890 | nP = 0 |
| 1161 | 98120 | do ! Loop to accumulate pairs of columns. |
| 1162 | 144010 | if ((kL > num_srt(i,j)) .or. (kR > num_srt(i,j+1))) exit |
| 1163 | ||
| 1164 | 98120 | if (rho_srt(i,kL,j) > rho_srt(i,kR,j+1)) then |
| 1165 | ! The right point is lighter and defines the density for this trio. | |
| 1166 | 4486 | nP = nP+1 ; k = nP |
| 1167 | 4486 | rho_pair = rho_srt(i,kR,j+1) |
| 1168 | ||
| 1169 | 4486 | k0b_Lv(i,k,J) = k0_srt(i,kL,j) ; k0b_Rv(i,k,J) = k0_srt(i,kR,j+1) |
| 1170 | 4486 | k0a_Lv(i,k,J) = k0_srt(i,kL-1,j) ; k0a_Rv(i,k,J) = k0b_Rv(i,k,J) |
| 1171 | 4486 | kbs_Lp(k) = kL ; kbs_Rp(k) = kR |
| 1172 | ||
| 1173 | 4486 | rho_a = rho_srt(i,kL-1,j) ; rho_b = rho_srt(i,kL,j) |
| 1174 | 4486 | wt_b = 1.0 ; if (abs(rho_a - rho_b) > abs(rho_pair - rho_a)) & |
| 1175 | 4486 | wt_b = (rho_pair - rho_a) / (rho_b - rho_a) |
| 1176 | 4486 | deep_wt_Lv(I,k,J) = wt_b ; deep_wt_Rv(I,k,J) = 1.0 |
| 1177 | ||
| 1178 | 4486 | h_demand_L(kL) = h_demand_L(kL) + 0.5*h_srt(i,kR,j+1) * wt_b |
| 1179 | 4486 | h_demand_L(kL-1) = h_demand_L(kL-1) + 0.5*h_srt(i,kR,j+1) * (1.0-wt_b) |
| 1180 | ||
| 1181 | 4486 | kR = kR+1 ; left_set(k) = .false. ; right_set(k) = .true. |
| 1182 | 93634 | elseif (rho_srt(i,kL,j) < rho_srt(i,kR,j+1)) then |
| 1183 | ! The left point is lighter and defines the density for this trio. | |
| 1184 | 43582 | nP = nP+1 ; k = nP |
| 1185 | 43582 | rho_pair = rho_srt(i,kL,j) |
| 1186 | 43582 | k0b_Lv(i,k,J) = k0_srt(i,kL,j) ; k0b_Rv(i,k,J) = k0_srt(i,kR,j+1) |
| 1187 | 43582 | k0a_Lv(i,k,J) = k0b_Lv(i,k,J) ; k0a_Rv(i,k,J) = k0_srt(i,kR-1,j+1) |
| 1188 | ||
| 1189 | 43582 | kbs_Lp(k) = kL ; kbs_Rp(k) = kR |
| 1190 | ||
| 1191 | 43582 | rho_a = rho_srt(i,kR-1,j+1) ; rho_b = rho_srt(i,kR,j+1) |
| 1192 | 43582 | wt_b = 1.0 ; if (abs(rho_a - rho_b) > abs(rho_pair - rho_a)) & |
| 1193 | 43582 | wt_b = (rho_pair - rho_a) / (rho_b - rho_a) |
| 1194 | 43582 | deep_wt_Lv(I,k,J) = 1.0 ; deep_wt_Rv(I,k,J) = wt_b |
| 1195 | ||
| 1196 | 43582 | h_demand_R(kR) = h_demand_R(kR) + 0.5*h_srt(i,kL,j) * wt_b |
| 1197 | 43582 | h_demand_R(kR-1) = h_demand_R(kR-1) + 0.5*h_srt(i,kL,j) * (1.0-wt_b) |
| 1198 | ||
| 1199 | 43582 | kL = kL+1 ; left_set(k) = .true. ; right_set(k) = .false. |
| 1200 | 50052 | elseif ((k0_srt(i,kL,j) <= nkmb) .or. (k0_srt(i,kR,j+1) <= nkmb)) then |
| 1201 | ! The densities are exactly equal and one layer is above the interior. | |
| 1202 | 0 | nP = nP+1 ; k = nP |
| 1203 | 0 | k0b_Lv(i,k,J) = k0_srt(i,kL,j) ; k0b_Rv(i,k,J) = k0_srt(i,kR,j+1) |
| 1204 | 0 | k0a_Lv(i,k,J) = k0b_Lv(i,k,J) ; k0a_Rv(i,k,J) = k0b_Rv(i,k,J) |
| 1205 | 0 | kbs_Lp(k) = kL ; kbs_Rp(k) = kR |
| 1206 | 0 | deep_wt_Lv(I,k,J) = 1.0 ; deep_wt_Rv(I,k,J) = 1.0 |
| 1207 | ||
| 1208 | 0 | h_demand_L(kL) = h_demand_L(kL) + 0.5*h_srt(i,kR,j+1) |
| 1209 | 0 | h_demand_R(kR) = h_demand_R(kR) + 0.5*h_srt(i,kL,j) |
| 1210 | ||
| 1211 | 0 | kL = kL+1 ; kR = kR+1 ; left_set(k) = .true. ; right_set(k) = .true. |
| 1212 | else ! The densities are exactly equal and in the interior. | |
| 1213 | ! Mixing in this case has already occurred, so accumulate the thickness | |
| 1214 | ! demanded for that mixing and skip onward. | |
| 1215 | 50052 | h_demand_L(kL) = h_demand_L(kL) + 0.5*h_srt(i,kR,j+1) |
| 1216 | 50052 | h_demand_R(kR) = h_demand_R(kR) + 0.5*h_srt(i,kL,j) |
| 1217 | ||
| 1218 | 50052 | kL = kL+1 ; kR = kR+1 |
| 1219 | endif | |
| 1220 | enddo ! Loop to accumulate pairs of columns. | |
| 1221 | 45890 | nPv(i,J) = nP ! This is the number of active pairings. |
| 1222 | ||
| 1223 | ! Determine what fraction of the thickness "demand" can be supplied. | |
| 1224 | 145814 | do k=1,num_srt(i,j+1) |
| 1225 | 99924 | h_supply_frac_R(k) = 1.0 |
| 1226 | 99924 | if (h_demand_R(k) > 0.5*h_srt(i,k,j+1)) & |
| 1227 | 70814 | h_supply_frac_R(k) = 0.5*h_srt(i,k,j+1) / h_demand_R(k) |
| 1228 | enddo | |
| 1229 | 145769 | do k=1,num_srt(i,j) |
| 1230 | 99879 | h_supply_frac_L(k) = 1.0 |
| 1231 | 99879 | if (h_demand_L(k) > 0.5*h_srt(i,k,j)) & |
| 1232 | 95711 | h_supply_frac_L(k) = 0.5*h_srt(i,k,j) / h_demand_L(k) |
| 1233 | enddo | |
| 1234 | ||
| 1235 | ! Distribute the "exported" thicknesses proportionately. | |
| 1236 | 93958 | do k=1,nPv(i,J) |
| 1237 | 48068 | kL = kbs_Lp(k) ; kR = kbs_Rp(k) |
| 1238 | 48068 | hP_Lv(I,k,J) = 0.0 ; hP_Rv(I,k,J) = 0.0 |
| 1239 | 48068 | if (left_set(k)) then ! Add the contributing thicknesses on the right. |
| 1240 | 43582 | if (deep_wt_Rv(I,k,J) < 1.0) then |
| 1241 | 43582 | hP_Rv(I,k,J) = 0.5*h_srt(i,kL,j) * min(h_supply_frac_R(kR), h_supply_frac_R(kR-1)) |
| 1242 | 43582 | wt_b = deep_wt_Rv(I,k,J) |
| 1243 | 43582 | h_used_R(kR-1) = h_used_R(kR-1) + (1.0 - wt_b) * hP_Rv(I,k,J) |
| 1244 | 43582 | h_used_R(kR) = h_used_R(kR) + wt_b * hP_Rv(I,k,J) |
| 1245 | else | |
| 1246 | 0 | hP_Rv(I,k,J) = 0.5*h_srt(i,kL,j) * h_supply_frac_R(kR) |
| 1247 | 0 | h_used_R(kR) = h_used_R(kR) + hP_Rv(I,k,J) |
| 1248 | endif | |
| 1249 | endif | |
| 1250 | 93958 | if (right_set(k)) then ! Add the contributing thicknesses on the left. |
| 1251 | 4486 | if (deep_wt_Lv(I,k,J) < 1.0) then |
| 1252 | 4486 | hP_Lv(I,k,J) = 0.5*h_srt(i,kR,j+1) * min(h_supply_frac_L(kL), h_supply_frac_L(kL-1)) |
| 1253 | 4486 | wt_b = deep_wt_Lv(I,k,J) |
| 1254 | 4486 | h_used_L(kL-1) = h_used_L(kL-1) + (1.0 - wt_b) * hP_Lv(I,k,J) |
| 1255 | 4486 | h_used_L(kL) = h_used_L(kL) + wt_b * hP_Lv(I,k,J) |
| 1256 | else | |
| 1257 | 0 | hP_Lv(I,k,J) = 0.5*h_srt(i,kR,j+1) * h_supply_frac_L(kL) |
| 1258 | 0 | h_used_L(kL) = h_used_L(kL) + hP_Lv(I,k,J) |
| 1259 | endif | |
| 1260 | endif | |
| 1261 | enddo | |
| 1262 | ||
| 1263 | ! The left-over thickness (at least half the layer thickness) is now | |
| 1264 | ! added to the thicknesses of the importing columns. | |
| 1265 | 93958 | do k=1,nPv(i,J) |
| 1266 | 48068 | if (left_set(k)) hP_Lv(I,k,J) = hP_Lv(I,k,J) + & |
| 1267 | 43582 | (h_srt(i,kbs_Lp(k),j) - h_used_L(kbs_Lp(k))) |
| 1268 | 48068 | if (right_set(k)) hP_Rv(I,k,J) = hP_Rv(I,k,J) + & |
| 1269 | 50376 | (h_srt(i,kbs_Rp(k),j+1) - h_used_R(kbs_Rp(k))) |
| 1270 | enddo | |
| 1271 | ||
| 1272 | ||
| 1273 | endif ; enddo ; enddo ! i- & j- loops over meridional faces. | |
| 1274 | !$omp target exit data map(release: h_srt, k0_srt) | |
| 1275 | ! The tracer-specific calculations start here. | |
| 1276 | ||
| 1277 | !$omp target enter data map(alloc: Tr_flux_3d, Tr_adj_vert_L, Tr_adj_vert_R, tr_flux_N, & | |
| 1278 | !$omp tr_flux_S, tr_flux_E, tr_flux_W, tr_flux_conv) | |
| 1279 | ||
| 1280 | 2 | do itt=1,max_itt |
| 1281 | ||
| 1282 | 1 | if (itt > 1) then ! The halos have already been filled if itt==1. |
| 1283 | 0 | call do_group_pass(CS%pass_t, G%Domain, clock=id_clock_pass, omp_offload=.true.) |
| 1284 | endif | |
| 1285 | 5 | do m=1,ntr |
| 1286 | ! Zero out tracer tendencies. | |
| 1287 | 3 | if (CS%answer_date <= 20240330) then |
| 1288 | 0 | do concurrent (k=1:nz, j=jsd:jed, i=isd:ied) |
| 1289 | 0 | tr_flux_conv(i,j,k) = 0.0 |
| 1290 | enddo | |
| 1291 | else | |
| 1292 | 3 | do concurrent (k=1:nz, j=jsd:jed, i=isd:ied) |
| 1293 | 4566144 | tr_flux_N(i,j,k) = 0.0 ; tr_flux_S(i,j,k) = 0.0 |
| 1294 | 4774803 | tr_flux_E(i,j,k) = 0.0 ; tr_flux_W(i,j,k) = 0.0 |
| 1295 | enddo | |
| 1296 | endif | |
| 1297 | 3 | do concurrent (k=1:2*nz, J=JsdB:JedB, i=isd:ied) |
| 1298 | 9180864 | tr_flux_3d(i,j,k) = 0.0 |
| 1299 | 9390627 | tr_adj_vert_R(i,j,k) = 0.0 ; tr_adj_vert_L(i,j,k) = 0.0 |
| 1300 | enddo | |
| 1301 | ||
| 1302 | ! collapse(2) is reproducible for CS%answer_date > 20240330 but not <= 20240330. | |
| 1303 | ! Do concurrent around j-loop doesn't seem to do the right thing. | |
| 1304 | !$omp target teams loop collapse(2) private(Tr_min_face,Tr_max_face,kLa,kLb,kRa,kRb, & | |
| 1305 | !$omp Tr_La,Tr_Lb,Tr_Ra,Tr_Rb,Tr_av_L,wt_b,Tr_av_R, & | |
| 1306 | !$omp h_L,h_R,Tr_flux,Tr_adj_vert,wt_a,vol) | |
| 1307 | 195483 | do j=js,je ; do I=is-1,ie ; if (G%mask2dCu(I,j) > 0.0) then |
| 1308 | ! Determine the fluxes through the zonal faces. | |
| 1309 | ||
| 1310 | ! Find the acceptable range of tracer concentration around this face. | |
| 1311 | 138201 | if (nPu(I,j) >= 1) then |
| 1312 | 131520 | Tr_min_face = min(Tr(m)%t(i,j,1), Tr(m)%t(i+1,j,1)) |
| 1313 | 131520 | Tr_max_face = max(Tr(m)%t(i,j,1), Tr(m)%t(i+1,j,1)) |
| 1314 | 526080 | do k=2,nkmb |
| 1315 | 394560 | Tr_min_face = min(Tr_min_face, Tr(m)%t(i,j,k), Tr(m)%t(i+1,j,k)) |
| 1316 | 526080 | Tr_max_face = max(Tr_max_face, Tr(m)%t(i,j,k), Tr(m)%t(i+1,j,k)) |
| 1317 | enddo | |
| 1318 | ||
| 1319 | ! Include the next two layers denser than the densest buffer layer. | |
| 1320 | 131520 | kLa = nkmb+1 ; if (max_kRho(i,j) < nz+1) kLa = max_kRho(i,j) |
| 1321 | 131520 | kLb = kLa ; if (max_kRho(i,j) < nz) kLb = max_kRho(i,j)+1 |
| 1322 | 131520 | kRa = nkmb+1 ; if (max_kRho(i+1,j) < nz+1) kRa = max_kRho(i+1,j) |
| 1323 | 131520 | kRb = kRa ; if (max_kRho(i+1,j) < nz) kRb = max_kRho(i+1,j)+1 |
| 1324 | 131520 | Tr_La = Tr_min_face ; Tr_Lb = Tr_La ; Tr_Ra = Tr_La ; Tr_Rb = Tr_La |
| 1325 | 131520 | if (h(i,j,kLa) > h_exclude) Tr_La = Tr(m)%t(i,j,kLa) |
| 1326 | 131520 | if ((CS%answer_date <= 20240330) .and. CS%limit_bug) then |
| 1327 | 0 | if (h(i,j,kLb) > h_exclude) Tr_La = Tr(m)%t(i,j,kLb) |
| 1328 | else | |
| 1329 | 131520 | if (h(i,j,kLb) > h_exclude) Tr_Lb = Tr(m)%t(i,j,kLb) |
| 1330 | endif | |
| 1331 | 131520 | if (h(i+1,j,kRa) > h_exclude) Tr_Ra = Tr(m)%t(i+1,j,kRa) |
| 1332 | 131520 | if (h(i+1,j,kRb) > h_exclude) Tr_Rb = Tr(m)%t(i+1,j,kRb) |
| 1333 | 131520 | Tr_min_face = min(Tr_min_face, Tr_La, Tr_Lb, Tr_Ra, Tr_Rb) |
| 1334 | 131520 | Tr_max_face = max(Tr_max_face, Tr_La, Tr_Lb, Tr_Ra, Tr_Rb) |
| 1335 | ||
| 1336 | ! Include all points in diffusive pairings at this face. | |
| 1337 | 284877 | do k=1,nPu(I,j) |
| 1338 | 153357 | Tr_Lb = Tr(m)%t(i,j,k0b_Lu(I,k,j)) |
| 1339 | 153357 | Tr_Rb = Tr(m)%t(i+1,j,k0b_Ru(I,k,j)) |
| 1340 | 153357 | Tr_La = Tr_Lb ; Tr_Ra = Tr_Rb |
| 1341 | 153357 | if (deep_wt_Lu(I,k,j) < 1.0) Tr_La = Tr(m)%t(i,j,k0a_Lu(I,k,j)) |
| 1342 | 153357 | if (deep_wt_Ru(I,k,j) < 1.0) Tr_Ra = Tr(m)%t(i+1,j,k0a_Ru(I,k,j)) |
| 1343 | 153357 | Tr_min_face = min(Tr_min_face, Tr_La, Tr_Lb, Tr_Ra, Tr_Rb) |
| 1344 | 284877 | Tr_max_face = max(Tr_max_face, Tr_La, Tr_Lb, Tr_Ra, Tr_Rb) |
| 1345 | enddo | |
| 1346 | endif | |
| 1347 | ||
| 1348 | 291558 | do k=1,nPu(I,j) |
| 1349 | 153357 | kLb = k0b_Lu(I,k,j) ; Tr_Lb = Tr(m)%t(i,j,kLb) ; Tr_av_L = Tr_Lb |
| 1350 | 153357 | if (deep_wt_Lu(I,k,j) < 1.0) then |
| 1351 | 76020 | kLa = k0a_Lu(I,k,j) ; Tr_La = Tr(m)%t(i,j,kLa) |
| 1352 | 76020 | wt_b = deep_wt_Lu(I,k,j) |
| 1353 | 76020 | Tr_av_L = wt_b*Tr_Lb + (1.0-wt_b)*Tr_La |
| 1354 | endif | |
| 1355 | ||
| 1356 | 153357 | kRb = k0b_Ru(I,k,j) ; Tr_Rb = Tr(m)%t(i+1,j,kRb) ; Tr_av_R = Tr_Rb |
| 1357 | 153357 | if (deep_wt_Ru(I,k,j) < 1.0) then |
| 1358 | 76374 | kRa = k0a_Ru(I,k,j) ; Tr_Ra = Tr(m)%t(i+1,j,kRa) |
| 1359 | 76374 | wt_b = deep_wt_Ru(I,k,j) |
| 1360 | 76374 | Tr_av_R = wt_b*Tr_Rb + (1.0-wt_b)*Tr_Ra |
| 1361 | endif | |
| 1362 | ||
| 1363 | 153357 | h_L = hP_Lu(I,k,j) ; h_R = hP_Ru(I,k,j) |
| 1364 | 153357 | if (CS%answer_date <= 20240330) then |
| 1365 | Tr_flux = I_maxitt * khdt_epi_x(I,j) * (Tr_av_L - Tr_av_R) * & | |
| 1366 | 0 | ((2.0 * h_L * h_R) / (h_L + h_R)) |
| 1367 | else | |
| 1368 | Tr_flux = I_maxitt * ((2.0 * h_L * h_R) / (h_L + h_R)) * & | |
| 1369 | 153357 | khdt_epi_x(I,j) * (Tr_av_L - Tr_av_R) |
| 1370 | endif | |
| 1371 | ||
| 1372 | 153357 | if (deep_wt_Lu(I,k,j) >= 1.0) then |
| 1373 | 77337 | if (CS%answer_date <= 20240330) then |
| 1374 | 0 | tr_flux_conv(i,j,kLb) = tr_flux_conv(i,j,kLb) - Tr_flux |
| 1375 | else | |
| 1376 | 77337 | tr_flux_E(i,j,kLb) = tr_flux_E(i,j,kLb) + Tr_flux |
| 1377 | endif | |
| 1378 | else | |
| 1379 | 76020 | Tr_adj_vert = 0.0 |
| 1380 | 76020 | wt_b = deep_wt_Lu(I,k,j) ; wt_a = 1.0 - wt_b |
| 1381 | 76020 | vol = hP_Lu(I,k,j) * G%areaT(i,j) |
| 1382 | ||
| 1383 | ! Ensure that the tracer flux does not drive the tracer values | |
| 1384 | ! outside of the range Tr_min_face <= Tr <= Tr_max_face, or if it | |
| 1385 | ! does that the concentration in both contributing pieces exceed | |
| 1386 | ! this range equally. With down-gradient fluxes and the initial tracer | |
| 1387 | ! concentrations determining the valid range, the latter condition | |
| 1388 | ! only enters for large values of the effective diffusive CFL number. | |
| 1389 | 76020 | if (Tr_flux > 0.0) then |
| 1390 | 7752 | if (Tr_La < Tr_Lb) then ; if (vol*(Tr_La-Tr_min_face) < Tr_flux) & |
| 1391 | Tr_adj_vert = -wt_a * min(Tr_flux - vol * (Tr_La-Tr_min_face), & | |
| 1392 | 299 | (vol*wt_b) * (Tr_Lb - Tr_La)) |
| 1393 | 6292 | else ; if (vol*(Tr_Lb-Tr_min_face) < Tr_flux) & |
| 1394 | Tr_adj_vert = wt_b * min(Tr_flux - vol * (Tr_Lb-Tr_min_face), & | |
| 1395 | 1494 | (vol*wt_a) * (Tr_La - Tr_Lb)) |
| 1396 | endif | |
| 1397 | 68268 | elseif (Tr_flux < 0.0) then |
| 1398 | 33058 | if (Tr_La > Tr_Lb) then ; if (vol * (Tr_max_face-Tr_La) < -Tr_flux) & |
| 1399 | Tr_adj_vert = wt_a * min(-Tr_flux - vol * (Tr_max_face-Tr_La), & | |
| 1400 | 21839 | (vol*wt_b) * (Tr_La - Tr_Lb)) |
| 1401 | 6306 | else ; if (vol*(Tr_max_face-Tr_Lb) < -Tr_flux) & |
| 1402 | Tr_adj_vert = -wt_b * min(-Tr_flux - vol * (Tr_max_face-Tr_Lb), & | |
| 1403 | 1584 | (vol*wt_a)*(Tr_Lb - Tr_La)) |
| 1404 | endif | |
| 1405 | endif | |
| 1406 | ||
| 1407 | 76020 | if (CS%answer_date <= 20240330) then |
| 1408 | 0 | tr_flux_conv(i,j,kLa) = tr_flux_conv(i,j,kLa) - (wt_a*Tr_flux + Tr_adj_vert) |
| 1409 | 0 | tr_flux_conv(i,j,kLb) = tr_flux_conv(i,j,kLb) - (wt_b*Tr_flux - Tr_adj_vert) |
| 1410 | else | |
| 1411 | 76020 | tr_flux_E(i,j,kLa) = tr_flux_E(i,j,kLa) + (wt_a*Tr_flux + Tr_adj_vert) |
| 1412 | 76020 | tr_flux_E(i,j,kLb) = tr_flux_E(i,j,kLb) + (wt_b*Tr_flux - Tr_adj_vert) |
| 1413 | endif | |
| 1414 | endif | |
| 1415 | ||
| 1416 | 153357 | if (deep_wt_Ru(I,k,j) >= 1.0) then |
| 1417 | 76983 | if (CS%answer_date <= 20240330) then |
| 1418 | 0 | tr_flux_conv(i+1,j,kRb) = tr_flux_conv(i+1,j,kRb) + Tr_flux |
| 1419 | else | |
| 1420 | 76983 | tr_flux_W(i+1,j,kRb) = tr_flux_W(i+1,j,kRb) + Tr_flux |
| 1421 | endif | |
| 1422 | else | |
| 1423 | 76374 | Tr_adj_vert = 0.0 |
| 1424 | 76374 | wt_b = deep_wt_Ru(I,k,j) ; wt_a = 1.0 - wt_b |
| 1425 | 76374 | vol = hP_Ru(I,k,j) * G%areaT(i+1,j) |
| 1426 | ||
| 1427 | ! Ensure that the tracer flux does not drive the tracer values | |
| 1428 | ! outside of the range Tr_min_face <= Tr <= Tr_max_face, or if it | |
| 1429 | ! does that the concentration in both contributing pieces exceed | |
| 1430 | ! this range equally. With down-gradient fluxes and the initial tracer | |
| 1431 | ! concentrations determining the valid range, the latter condition | |
| 1432 | ! only enters for large values of the effective diffusive CFL number. | |
| 1433 | 76374 | if (Tr_flux < 0.0) then |
| 1434 | 7808 | if (Tr_Ra < Tr_Rb) then ; if (vol * (Tr_Ra-Tr_min_face) < -Tr_flux) & |
| 1435 | Tr_adj_vert = -wt_a * min(-Tr_flux - vol * (Tr_Ra-Tr_min_face), & | |
| 1436 | 283 | (vol*wt_b) * (Tr_Rb - Tr_Ra)) |
| 1437 | 6384 | else ; if (vol*(Tr_Rb-Tr_min_face) < (-Tr_flux)) & |
| 1438 | Tr_adj_vert = wt_b * min(-Tr_flux - vol * (Tr_Rb-Tr_min_face), & | |
| 1439 | 1499 | (vol*wt_a) * (Tr_Ra - Tr_Rb)) |
| 1440 | endif | |
| 1441 | 68566 | elseif (Tr_flux > 0.0) then |
| 1442 | 33228 | if (Tr_Ra > Tr_Rb) then ; if (vol * (Tr_max_face-Tr_Ra) < Tr_flux) & |
| 1443 | Tr_adj_vert = wt_a * min(Tr_flux - vol * (Tr_max_face-Tr_Ra), & | |
| 1444 | 21891 | (vol*wt_b) * (Tr_Ra - Tr_Rb)) |
| 1445 | 6432 | else ; if (vol*(Tr_max_face-Tr_Rb) < Tr_flux) & |
| 1446 | Tr_adj_vert = -wt_b * min(Tr_flux - vol * (Tr_max_face-Tr_Rb), & | |
| 1447 | 1602 | (vol*wt_a)*(Tr_Rb - Tr_Ra)) |
| 1448 | endif | |
| 1449 | endif | |
| 1450 | ||
| 1451 | 76374 | if (CS%answer_date <= 20240330) then |
| 1452 | 0 | tr_flux_conv(i+1,j,kRa) = tr_flux_conv(i+1,j,kRa) + (wt_a*Tr_flux - Tr_adj_vert) |
| 1453 | 0 | tr_flux_conv(i+1,j,kRb) = tr_flux_conv(i+1,j,kRb) + (wt_b*Tr_flux + Tr_adj_vert) |
| 1454 | else | |
| 1455 | 76374 | tr_flux_W(i+1,j,kRa) = tr_flux_W(i+1,j,kRa) + (wt_a*Tr_flux - Tr_adj_vert) |
| 1456 | 76374 | tr_flux_W(i+1,j,kRb) = tr_flux_W(i+1,j,kRb) + (wt_b*Tr_flux + Tr_adj_vert) |
| 1457 | endif | |
| 1458 | endif | |
| 1459 | 153357 | if (associated(Tr(m)%df2d_x)) & |
| 1460 | 138201 | Tr(m)%df2d_x(I,j) = Tr(m)%df2d_x(I,j) + Tr_flux * Idt |
| 1461 | enddo ! Loop over pairings at faces. | |
| 1462 | endif ; enddo ; enddo ! i- & j- loops over zonal faces. | |
| 1463 | ||
| 1464 | ! this gives wrong result when using do concurrent on NVHPC 25.9 | |
| 1465 | !$omp target teams loop collapse(2) private(Tr_min_face,Tr_max_face,kLa,kLb,kRa,kRb, & | |
| 1466 | !$omp Tr_La,Tr_Lb,Tr_Ra,Tr_Rb,Tr_av_L,wt_b,Tr_av_R, & | |
| 1467 | !$omp h_L,h_R,Tr_flux,Tr_adj_vert,wt_a,vol) | |
| 1468 | 196026 | do J=js-1,je ; do i=is,ie ; if (G%mask2dCv(i,J) > 0.0) then |
| 1469 | ! Determine the fluxes through the meridional faces. | |
| 1470 | ||
| 1471 | ! Find the acceptable range of tracer concentration around this face. | |
| 1472 | 137670 | if (nPv(i,J) >= 1) then |
| 1473 | 130746 | Tr_min_face = min(Tr(m)%t(i,j,1), Tr(m)%t(i,j+1,1)) |
| 1474 | 130746 | Tr_max_face = max(Tr(m)%t(i,j,1), Tr(m)%t(i,j+1,1)) |
| 1475 | 522984 | do k=2,nkmb |
| 1476 | 392238 | Tr_min_face = min(Tr_min_face, Tr(m)%t(i,j,k), Tr(m)%t(i,j+1,k)) |
| 1477 | 522984 | Tr_max_face = max(Tr_max_face, Tr(m)%t(i,j,k), Tr(m)%t(i,j+1,k)) |
| 1478 | enddo | |
| 1479 | ||
| 1480 | ! Include the next two layers denser than the densest buffer layer. | |
| 1481 | 130746 | kLa = nkmb+1 ; if (max_kRho(i,j) < nz+1) kLa = max_kRho(i,j) |
| 1482 | 130746 | kLb = kLa ; if (max_kRho(i,j) < nz) kLb = max_kRho(i,j)+1 |
| 1483 | 130746 | kRa = nkmb+1 ; if (max_kRho(i,j+1) < nz+1) kRa = max_kRho(i,j+1) |
| 1484 | 130746 | kRb = kRa ; if (max_kRho(i,j+1) < nz) kRb = max_kRho(i,j+1)+1 |
| 1485 | 130746 | Tr_La = Tr_min_face ; Tr_Lb = Tr_La ; Tr_Ra = Tr_La ; Tr_Rb = Tr_La |
| 1486 | 130746 | if (h(i,j,kLa) > h_exclude) Tr_La = Tr(m)%t(i,j,kLa) |
| 1487 | 130746 | if ((CS%answer_date <= 20240330) .and. CS%limit_bug) then |
| 1488 | 0 | if (h(i,j,kLb) > h_exclude) Tr_La = Tr(m)%t(i,j,kLb) |
| 1489 | else | |
| 1490 | 130746 | if (h(i,j,kLb) > h_exclude) Tr_Lb = Tr(m)%t(i,j,kLb) |
| 1491 | endif | |
| 1492 | 130746 | if (h(i,j+1,kRa) > h_exclude) Tr_Ra = Tr(m)%t(i,j+1,kRa) |
| 1493 | 130746 | if (h(i,j+1,kRb) > h_exclude) Tr_Rb = Tr(m)%t(i,j+1,kRb) |
| 1494 | 130746 | Tr_min_face = min(Tr_min_face, Tr_La, Tr_Lb, Tr_Ra, Tr_Rb) |
| 1495 | 130746 | Tr_max_face = max(Tr_max_face, Tr_La, Tr_Lb, Tr_Ra, Tr_Rb) |
| 1496 | ||
| 1497 | ! Include all points in diffusive pairings at this face. | |
| 1498 | 274950 | do k=1,nPv(i,J) |
| 1499 | 144204 | Tr_Lb = Tr(m)%t(i,j,k0b_Lv(i,k,J)) ; Tr_Rb = Tr(m)%t(i,j+1,k0b_Rv(i,k,J)) |
| 1500 | 144204 | Tr_La = Tr_Lb ; Tr_Ra = Tr_Rb |
| 1501 | 144204 | if (deep_wt_Lv(I,k,J) < 1.0) Tr_La = Tr(m)%t(i,j,k0a_Lv(i,k,J)) |
| 1502 | 144204 | if (deep_wt_Rv(I,k,J) < 1.0) Tr_Ra = Tr(m)%t(i,j+1,k0a_Rv(i,k,J)) |
| 1503 | 144204 | Tr_min_face = min(Tr_min_face, Tr_La, Tr_Lb, Tr_Ra, Tr_Rb) |
| 1504 | 274950 | Tr_max_face = max(Tr_max_face, Tr_La, Tr_Lb, Tr_Ra, Tr_Rb) |
| 1505 | enddo | |
| 1506 | endif | |
| 1507 | ||
| 1508 | 281874 | do k=1,nPv(i,J) |
| 1509 | 144204 | kLb = k0b_Lv(i,k,J) ; Tr_Lb = Tr(m)%t(i,j,kLb) ; Tr_av_L = Tr_Lb |
| 1510 | 144204 | if (deep_wt_Lv(I,k,J) < 1.0) then |
| 1511 | 13458 | kLa = k0a_Lv(i,k,J) ; Tr_La = Tr(m)%t(i,j,kLa) |
| 1512 | 13458 | wt_b = deep_wt_Lv(I,k,J) |
| 1513 | 13458 | Tr_av_L = wt_b * Tr_Lb + (1.0-wt_b) * Tr_La |
| 1514 | endif | |
| 1515 | ||
| 1516 | 144204 | kRb = k0b_Rv(i,k,J) ; Tr_Rb = Tr(m)%t(i,j+1,kRb) ; Tr_av_R = Tr_Rb |
| 1517 | 144204 | if (deep_wt_Rv(I,k,J) < 1.0) then |
| 1518 | 130746 | kRa = k0a_Rv(i,k,J) ; Tr_Ra = Tr(m)%t(i,j+1,kRa) |
| 1519 | 130746 | wt_b = deep_wt_Rv(I,k,J) |
| 1520 | 130746 | Tr_av_R = wt_b * Tr_Rb + (1.0-wt_b) * Tr_Ra |
| 1521 | endif | |
| 1522 | ||
| 1523 | 144204 | h_L = hP_Lv(I,k,J) ; h_R = hP_Rv(I,k,J) |
| 1524 | Tr_flux = I_maxitt * ((2.0 * h_L * h_R) / (h_L + h_R)) * & | |
| 1525 | 144204 | khdt_epi_y(i,J) * (Tr_av_L - Tr_av_R) |
| 1526 | 144204 | Tr_flux_3d(i,J,k) = Tr_flux |
| 1527 | ||
| 1528 | 144204 | if (deep_wt_Lv(I,k,J) < 1.0) then |
| 1529 | 13458 | Tr_adj_vert = 0.0 |
| 1530 | 13458 | wt_b = deep_wt_Lv(I,k,J) ; wt_a = 1.0 - wt_b |
| 1531 | 13458 | vol = hP_Lv(I,k,J) * G%areaT(i,j) |
| 1532 | ||
| 1533 | ! Ensure that the tracer flux does not drive the tracer values | |
| 1534 | ! outside of the range Tr_min_face <= Tr <= Tr_max_face. | |
| 1535 | 13458 | if (Tr_flux > 0.0) then |
| 1536 | 1214 | if (Tr_La < Tr_Lb) then ; if (vol * (Tr_La-Tr_min_face) < Tr_flux) & |
| 1537 | Tr_adj_vert = -wt_a * min(Tr_flux - vol * (Tr_La-Tr_min_face), & | |
| 1538 | 120 | (vol*wt_b) * (Tr_Lb - Tr_La)) |
| 1539 | 861 | else ; if (vol*(Tr_Lb-Tr_min_face) < Tr_flux) & |
| 1540 | Tr_adj_vert = wt_b * min(Tr_flux - vol * (Tr_Lb-Tr_min_face), & | |
| 1541 | 129 | (vol*wt_a) * (Tr_La - Tr_Lb)) |
| 1542 | endif | |
| 1543 | 12244 | elseif (Tr_flux < 0.0) then |
| 1544 | 5705 | if (Tr_La > Tr_Lb) then ; if (vol * (Tr_max_face-Tr_La) < -Tr_flux) & |
| 1545 | Tr_adj_vert = wt_a * min(-Tr_flux - vol * (Tr_max_face-Tr_La), & | |
| 1546 | 4437 | (vol*wt_b) * (Tr_La - Tr_Lb)) |
| 1547 | 868 | else ; if (vol*(Tr_max_face-Tr_Lb) < -Tr_flux) & |
| 1548 | Tr_adj_vert = -wt_b * min(-Tr_flux - vol * (Tr_max_face-Tr_Lb), & | |
| 1549 | 110 | (vol*wt_a)*(Tr_Lb - Tr_La)) |
| 1550 | endif | |
| 1551 | endif | |
| 1552 | 13458 | Tr_adj_vert_L(i,J,k) = Tr_adj_vert |
| 1553 | endif | |
| 1554 | ||
| 1555 | 144204 | if (deep_wt_Rv(I,k,J) < 1.0) then |
| 1556 | 130746 | Tr_adj_vert = 0.0 |
| 1557 | 130746 | wt_b = deep_wt_Rv(I,k,J) ; wt_a = 1.0 - wt_b |
| 1558 | 130746 | vol = hP_Rv(I,k,J) * G%areaT(i,j+1) |
| 1559 | ||
| 1560 | ! Ensure that the tracer flux does not drive the tracer values | |
| 1561 | ! outside of the range Tr_min_face <= Tr <= Tr_max_face. | |
| 1562 | 130746 | if (Tr_flux < 0.0) then |
| 1563 | 12828 | if (Tr_Ra < Tr_Rb) then ; if (vol * (Tr_Ra-Tr_min_face) < -Tr_flux) & |
| 1564 | Tr_adj_vert = -wt_a * min(-Tr_flux - vol * (Tr_Ra-Tr_min_face), & | |
| 1565 | 478 | (vol*wt_b) * (Tr_Rb - Tr_Ra)) |
| 1566 | 10461 | else ; if (vol*(Tr_Rb-Tr_min_face) < (-Tr_flux)) & |
| 1567 | Tr_adj_vert = wt_b * min(-Tr_flux - vol * (Tr_Rb-Tr_min_face), & | |
| 1568 | 2201 | (vol*wt_a) * (Tr_Ra - Tr_Rb)) |
| 1569 | endif | |
| 1570 | 117918 | elseif (Tr_flux > 0.0) then |
| 1571 | 56086 | if (Tr_Ra > Tr_Rb) then ; if (vol * (Tr_max_face-Tr_Ra) < Tr_flux) & |
| 1572 | Tr_adj_vert = wt_a * min(Tr_flux - vol * (Tr_max_face-Tr_Ra), & | |
| 1573 | 36185 | (vol*wt_b) * (Tr_Ra - Tr_Rb)) |
| 1574 | 10215 | else ; if (vol*(Tr_max_face-Tr_Rb) < Tr_flux) & |
| 1575 | Tr_adj_vert = -wt_b * min(Tr_flux - vol * (Tr_max_face-Tr_Rb), & | |
| 1576 | 2223 | (vol*wt_a)*(Tr_Rb - Tr_Ra)) |
| 1577 | endif | |
| 1578 | endif | |
| 1579 | 130746 | Tr_adj_vert_R(i,J,k) = Tr_adj_vert |
| 1580 | endif | |
| 1581 | 144204 | if (associated(Tr(m)%df2d_y)) & |
| 1582 | 137670 | Tr(m)%df2d_y(i,J) = Tr(m)%df2d_y(i,J) + Tr_flux * Idt |
| 1583 | enddo ! Loop over pairings at faces. | |
| 1584 | endif ; enddo ; enddo ! i- & j- loops over meridional faces. | |
| 1585 | ||
| 1586 | ! The non-stride-1 loop order here is to facilitate openMP threading. However, it might be | |
| 1587 | ! suboptimal when openMP threading is not used, at which point it might be better to fuse | |
| 1588 | ! this loop with those that precede it and thereby eliminate the need for three 3-d arrays. | |
| 1589 | 3 | if (CS%answer_date <= 20240330) then |
| 1590 | ! KRa/b aren't guaranteed to be unique and update of tr_flux_conv(:, j/j+1, :) means this | |
| 1591 | ! loop must be serial in both j and k. i can be paralellised. | |
| 1592 | 0 | do concurrent (i=is:ie) DO_LOCALITY(local(j, k, kLb, kRb, kLa, wt_b, wt_a, kRa)) |
| 1593 | 0 | do J=js-1,je ; if (G%mask2dCv(i,J) > 0.0) then |
| 1594 | 0 | do k=1,nPv(i,J) |
| 1595 | 0 | kLb = k0b_Lv(i,k,J) ; kRb = k0b_Rv(i,k,J) |
| 1596 | 0 | if (deep_wt_Lv(I,k,J) >= 1.0) then |
| 1597 | 0 | tr_flux_conv(i,j,kLb) = tr_flux_conv(i,j,kLb) - Tr_flux_3d(i,J,k) |
| 1598 | else | |
| 1599 | 0 | kLa = k0a_Lv(i,k,J) |
| 1600 | 0 | wt_b = deep_wt_Lv(I,k,J) ; wt_a = 1.0 - wt_b |
| 1601 | 0 | tr_flux_conv(i,j,kLa) = tr_flux_conv(i,j,kLa) - (wt_a*Tr_flux_3d(i,J,k) + Tr_adj_vert_L(i,J,k)) |
| 1602 | 0 | tr_flux_conv(i,j,kLb) = tr_flux_conv(i,j,kLb) - (wt_b*Tr_flux_3d(i,J,k) - Tr_adj_vert_L(i,J,k)) |
| 1603 | endif | |
| 1604 | 0 | if (deep_wt_Rv(I,k,J) >= 1.0) then |
| 1605 | 0 | tr_flux_conv(i,j+1,kRb) = tr_flux_conv(i,j+1,kRb) + Tr_flux_3d(i,J,k) |
| 1606 | else | |
| 1607 | 0 | kRa = k0a_Rv(i,k,J) |
| 1608 | 0 | wt_b = deep_wt_Rv(I,k,J) ; wt_a = 1.0 - wt_b |
| 1609 | tr_flux_conv(i,j+1,kRa) = tr_flux_conv(i,j+1,kRa) + & | |
| 1610 | 0 | (wt_a*Tr_flux_3d(i,J,k) - Tr_adj_vert_R(i,J,k)) |
| 1611 | tr_flux_conv(i,j+1,kRb) = tr_flux_conv(i,j+1,kRb) + & | |
| 1612 | 0 | (wt_b*Tr_flux_3d(i,J,k) + Tr_adj_vert_R(i,J,k)) |
| 1613 | endif | |
| 1614 | enddo | |
| 1615 | endif ; enddo | |
| 1616 | enddo | |
| 1617 | else | |
| 1618 | ! Update of tr_flux_N/S can be done independently so both i and j can be parallelised. | |
| 1619 | 196563 | do concurrent (J=js-1:je, i=is:ie, G%mask2dCv(i,J) > 0.0) DO_LOCALITY(local(k, kLb, kRb, kLa, wt_b, wt_a, kRa)) |
| 1620 | 340767 | do k=1,nPv(i,J) |
| 1621 | 144204 | kLb = k0b_Lv(i,k,J) ; kRb = k0b_Rv(i,k,J) |
| 1622 | 144204 | if (deep_wt_Lv(I,k,J) >= 1.0) then |
| 1623 | 130746 | tr_flux_N(i,j,kLb) = tr_flux_N(i,j,kLb) + Tr_flux_3d(i,J,k) |
| 1624 | else | |
| 1625 | 13458 | kLa = k0a_Lv(i,k,J) |
| 1626 | 13458 | wt_b = deep_wt_Lv(I,k,J) ; wt_a = 1.0 - wt_b |
| 1627 | 13458 | tr_flux_N(i,j,kLa) = tr_flux_N(i,j,kLa) + (wt_a*Tr_flux_3d(i,J,k) + Tr_adj_vert_L(i,J,k)) |
| 1628 | 13458 | tr_flux_N(i,j,kLb) = tr_flux_N(i,j,kLb) + (wt_b*Tr_flux_3d(i,J,k) - Tr_adj_vert_L(i,J,k)) |
| 1629 | endif | |
| 1630 | 281874 | if (deep_wt_Rv(I,k,J) >= 1.0) then |
| 1631 | 13458 | tr_flux_S(i,j+1,kRb) = tr_flux_S(i,j+1,kRb) + Tr_flux_3d(i,J,k) |
| 1632 | else | |
| 1633 | 130746 | kRa = k0a_Rv(i,k,J) |
| 1634 | 130746 | wt_b = deep_wt_Rv(I,k,J) ; wt_a = 1.0 - wt_b |
| 1635 | 130746 | tr_flux_S(i,j+1,kRa) = tr_flux_S(i,j+1,kRa) + (wt_a*Tr_flux_3d(i,J,k) - Tr_adj_vert_R(i,J,k)) |
| 1636 | 130746 | tr_flux_S(i,j+1,kRb) = tr_flux_S(i,j+1,kRb) + (wt_b*Tr_flux_3d(i,J,k) + Tr_adj_vert_R(i,J,k)) |
| 1637 | endif | |
| 1638 | enddo | |
| 1639 | enddo | |
| 1640 | endif | |
| 1641 | ||
| 1642 | 3 | if (CS%answer_date >= 20240331) then |
| 1643 | 3 | do concurrent (k=1:PEmax_kRho, j=js:je, i=is:ie) |
| 1644 | tr_flux_conv(i,j,k) = ((tr_flux_W(i,j,k) - tr_flux_E(i,j,k)) + & | |
| 1645 | 4083483 | (tr_flux_S(i,j,k) - tr_flux_N(i,j,k))) |
| 1646 | enddo | |
| 1647 | endif | |
| 1648 | ||
| 1649 | 4083483 | do concurrent (k=1:PEmax_kRho, j=js:je, i=is:ie, (G%mask2dT(i,j) > 0.0) .and. (h(i,j,k) > 0.0)) |
| 1650 | 4083483 | Tr(m)%t(i,j,k) = Tr(m)%t(i,j,k) + tr_flux_conv(i,j,k) / (h(i,j,k)*G%areaT(i,j)) |
| 1651 | enddo | |
| 1652 | ||
| 1653 | ! Do user controlled underflow of the tracer concentrations. | |
| 1654 | 4 | if (Tr(m)%conc_underflow > 0.0) then |
| 1655 | !$omp target update from(Tr(m)%t) | |
| 1656 | !$OMP parallel do default(shared) | |
| 1657 | 0 | do k=1,nz ; do j=js,je ; do i=is,ie |
| 1658 | 0 | if (abs(Tr(m)%t(i,j,k)) < Tr(m)%conc_underflow) Tr(m)%t(i,j,k) = 0.0 |
| 1659 | enddo ; enddo ; enddo | |
| 1660 | !$omp target update to(Tr(m)%t) | |
| 1661 | endif | |
| 1662 | ||
| 1663 | enddo ! Loop over tracers | |
| 1664 | enddo ! Loop over iterations | |
| 1665 | ||
| 1666 | !$omp target exit data map(release: Tr_flux_3d, Tr_adj_vert_L, Tr_adj_vert_R, tr_flux_N, & | |
| 1667 | !$omp tr_flux_S, tr_flux_E, tr_flux_W, tr_flux_conv, nPv, nPu, max_kRho, rho_srt, num_srt, & | |
| 1668 | !$omp max_srt) | |
| 1669 | ||
| 1670 | !$omp target exit data map(release: deep_wt_Lu, deep_wt_Ru, hP_Lu, hP_Ru, k0a_Lu, k0a_Ru, & | |
| 1671 | !$omp k0b_Lu, k0b_Ru) | |
| 1672 | ||
| 1673 | !$omp target exit data map(release: deep_wt_Lv, deep_wt_Rv, hP_Lv, hP_Rv, k0a_Lv, k0a_Rv, & | |
| 1674 | !$Omp k0b_Lv, k0b_Rv) | |
| 1675 | ||
| 1676 | 1 | end subroutine tracer_epipycnal_ML_diff |
| 1677 | ||
| 1678 | ||
| 1679 | !> Initialize lateral tracer diffusion module | |
| 1680 | 1 | subroutine tracer_hor_diff_init(Time, G, GV, US, param_file, diag, EOS, diabatic_CSp, CS) |
| 1681 | type(time_type), target, intent(in) :: Time !< current model time | |
| 1682 | type(ocean_grid_type), intent(in) :: G !< ocean grid structure | |
| 1683 | type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure | |
| 1684 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 1685 | type(diag_ctrl), target, intent(inout) :: diag !< diagnostic control | |
| 1686 | type(EOS_type), target, intent(in) :: EOS !< Equation of state CS | |
| 1687 | type(diabatic_CS), pointer, intent(in) :: diabatic_CSp !< Equation of state CS | |
| 1688 | type(param_file_type), intent(in) :: param_file !< parameter file | |
| 1689 | type(tracer_hor_diff_CS), pointer :: CS !< horz diffusion control structure | |
| 1690 | ||
| 1691 | ! This include declares and sets the variable "version". | |
| 1692 | # include "version_variable.h" | |
| 1693 | character(len=40) :: mdl = "MOM_tracer_hor_diff" ! This module's name. | |
| 1694 | integer :: default_answer_date | |
| 1695 | ||
| 1696 | 1 | if (associated(CS)) then |
| 1697 | 0 | call MOM_error(WARNING, "tracer_hor_diff_init called with associated control structure.") |
| 1698 | 0 | return |
| 1699 | endif | |
| 1700 | 1 | allocate(CS) |
| 1701 | ||
| 1702 | 1 | CS%diag => diag |
| 1703 | 1 | CS%show_call_tree = callTree_showQuery() |
| 1704 | ||
| 1705 | ! Read all relevant parameters and write them to the model log. | |
| 1706 | 1 | call log_version(param_file, mdl, version, "") |
| 1707 | call get_param(param_file, mdl, "KHTR", CS%KhTr, & | |
| 1708 | "The background along-isopycnal tracer diffusivity.", & | |
| 1709 | 1 | units="m2 s-1", default=0.0, scale=US%m_to_L**2*US%T_to_s) |
| 1710 | ! call get_param(param_file, mdl, "KHTR_USE_EBT_STRUCT", CS%KhTh_use_ebt_struct, & | |
| 1711 | ! "If true, uses the equivalent barotropic structure "//& | |
| 1712 | ! "as the vertical structure of the tracer diffusivity.",& | |
| 1713 | ! default=.false.) | |
| 1714 | call get_param(param_file, mdl, "KHTR_SLOPE_CFF", CS%KhTr_Slope_Cff, & | |
| 1715 | "The scaling coefficient for along-isopycnal tracer "//& | |
| 1716 | "diffusivity using a shear-based (Visbeck-like) "//& | |
| 1717 | "parameterization. A non-zero value enables this param.", & | |
| 1718 | 1 | units="nondim", default=0.0) |
| 1719 | call get_param(param_file, mdl, "KHTR_MIN", CS%KhTr_Min, & | |
| 1720 | "The minimum along-isopycnal tracer diffusivity.", & | |
| 1721 | 1 | units="m2 s-1", default=0.0, scale=US%m_to_L**2*US%T_to_s) |
| 1722 | 1 | if (CS%KhTr_Min > 0.0) then |
| 1723 | call get_param(param_file, mdl, "FULL_DEPTH_KHTR_MIN", CS%full_depth_khtr_min, & | |
| 1724 | "If true, KHTR_MIN is enforced throughout the whole water column. "//& | |
| 1725 | "Otherwise, KHTR_MIN is only enforced at the surface. This parameter "//& | |
| 1726 | "is only available when KHTR_USE_EBT_STRUCT=True and KHTR_MIN>0.", & | |
| 1727 | 0 | default=.false.) |
| 1728 | endif | |
| 1729 | call get_param(param_file, mdl, "KHTR_MAX", CS%KhTr_Max, & | |
| 1730 | "The maximum along-isopycnal tracer diffusivity.", & | |
| 1731 | 1 | units="m2 s-1", default=0.0, scale=US%m_to_L**2*US%T_to_s) |
| 1732 | call get_param(param_file, mdl, "KHTR_PASSIVITY_COEFF", CS%KhTr_passivity_coeff, & | |
| 1733 | "The coefficient that scales deformation radius over "//& | |
| 1734 | "grid-spacing in passivity, where passivity is the ratio "//& | |
| 1735 | "between along isopycnal mixing of tracers to thickness mixing. "//& | |
| 1736 | "A non-zero value enables this parameterization.", & | |
| 1737 | 1 | units="nondim", default=0.0) |
| 1738 | call get_param(param_file, mdl, "KHTR_PASSIVITY_MIN", CS%KhTr_passivity_min, & | |
| 1739 | "The minimum passivity which is the ratio between "//& | |
| 1740 | "along isopycnal mixing of tracers to thickness mixing.", & | |
| 1741 | 1 | units="nondim", default=0.5) |
| 1742 | call get_param(param_file, mdl, "DIFFUSE_ML_TO_INTERIOR", CS%Diffuse_ML_interior, & | |
| 1743 | "If true, enable epipycnal mixing between the surface "//& | |
| 1744 | 1 | "boundary layer and the interior.", default=.false.) |
| 1745 | call get_param(param_file, mdl, "CHECK_DIFFUSIVE_CFL", CS%check_diffusive_CFL, & | |
| 1746 | "If true, use enough iterations the diffusion to ensure "//& | |
| 1747 | "that the diffusive equivalent of the CFL limit is not "//& | |
| 1748 | "violated. If false, always use the greater of 1 or "//& | |
| 1749 | 1 | "MAX_TR_DIFFUSION_CFL iteration.", default=.false.) |
| 1750 | call get_param(param_file, mdl, "MAX_TR_DIFFUSION_CFL", CS%max_diff_CFL, & | |
| 1751 | "If positive, locally limit the along-isopycnal tracer "//& | |
| 1752 | "diffusivity to keep the diffusive CFL locally at or "//& | |
| 1753 | "below this value. The number of diffusive iterations "//& | |
| 1754 | "is often this value or the next greater integer.", & | |
| 1755 | 1 | units="nondim", default=-1.0) |
| 1756 | call get_param(param_File, mdl, "RECALC_NEUTRAL_SURF", CS%recalc_neutral_surf, & | |
| 1757 | "If true, then recalculate the neutral surfaces if the \n"//& | |
| 1758 | "diffusive CFL is exceeded. If false, assume that the \n"//& | |
| 1759 | 1 | "positions of the surfaces do not change \n", default=.false.) |
| 1760 | call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, & | |
| 1761 | "This sets the default value for the various _ANSWER_DATE parameters.", & | |
| 1762 | 1 | default=99991231, do_not_log=.true.) |
| 1763 | call get_param(param_file, mdl, "HOR_DIFF_ANSWER_DATE", CS%answer_date, & | |
| 1764 | "The vintage of the order of arithmetic to use for the tracer diffusion. "//& | |
| 1765 | "Values of 20240330 or below recover the answers from the original form of the "//& | |
| 1766 | "along-isopycnal mixed layer to interior mixing code, while higher values use "//& | |
| 1767 | "mathematically equivalent expressions that recover rotational symmetry "//& | |
| 1768 | "when DIFFUSE_ML_TO_INTERIOR is true.", & | |
| 1769 | 1 | default=default_answer_date, do_not_log=.not.CS%Diffuse_ML_interior) |
| 1770 | call get_param(param_file, mdl, "HOR_DIFF_LIMIT_BUG", CS%limit_bug, & | |
| 1771 | "If true and the answer date is 20240330 or below, use a rotational symmetry "//& | |
| 1772 | "breaking bug when limiting the tracer properties in tracer_epipycnal_ML_diff.", & | |
| 1773 | 1 | default=.false., do_not_log=((.not.CS%Diffuse_ML_interior).or.(CS%answer_date>=20240331))) |
| 1774 | 1 | CS%ML_KhTR_scale = 1.0 |
| 1775 | 1 | if (CS%Diffuse_ML_interior) then |
| 1776 | call get_param(param_file, mdl, "ML_KHTR_SCALE", CS%ML_KhTR_scale, & | |
| 1777 | "With Diffuse_ML_interior, the ratio of the truly "//& | |
| 1778 | "horizontal diffusivity in the mixed layer to the "//& | |
| 1779 | "epipycnal diffusivity. The valid range is 0 to 1.", & | |
| 1780 | 1 | units="nondim", default=1.0) |
| 1781 | endif | |
| 1782 | ||
| 1783 | CS%use_neutral_diffusion = neutral_diffusion_init(Time, G, GV, US, param_file, diag, EOS, & | |
| 1784 | 1 | diabatic_CSp, CS%neutral_diffusion_CSp ) |
| 1785 | 1 | if (CS%use_neutral_diffusion .and. CS%Diffuse_ML_interior) call MOM_error(FATAL, "MOM_tracer_hor_diff: "// & |
| 1786 | 0 | "USE_NEUTRAL_DIFFUSION and DIFFUSE_ML_TO_INTERIOR are mutually exclusive!") |
| 1787 | CS%use_hor_bnd_diffusion = hor_bnd_diffusion_init(Time, G, GV, US, param_file, diag, diabatic_CSp, & | |
| 1788 | 1 | CS%hor_bnd_diffusion_CSp) |
| 1789 | 1 | if (CS%use_hor_bnd_diffusion .and. CS%Diffuse_ML_interior) call MOM_error(FATAL, "MOM_tracer_hor_diff: "// & |
| 1790 | 0 | "USE_HORIZONTAL_BOUNDARY_DIFFUSION and DIFFUSE_ML_TO_INTERIOR are mutually exclusive!") |
| 1791 | ||
| 1792 | 1 | call get_param(param_file, mdl, "DEBUG", CS%debug, default=.false.) |
| 1793 | ||
| 1794 | 1 | id_clock_diffuse = cpu_clock_id('(Ocean diffuse tracer)', grain=CLOCK_MODULE) |
| 1795 | 1 | id_clock_epimix = cpu_clock_id('(Ocean epipycnal diffuse tracer)',grain=CLOCK_MODULE) |
| 1796 | 1 | id_clock_pass = cpu_clock_id('(Ocean tracer halo updates)', grain=CLOCK_ROUTINE) |
| 1797 | 1 | id_clock_sync = cpu_clock_id('(Ocean tracer global synch)', grain=CLOCK_ROUTINE) |
| 1798 | ||
| 1799 | 1 | CS%id_KhTr_u = -1 |
| 1800 | 1 | CS%id_KhTr_v = -1 |
| 1801 | 1 | CS%id_KhTr_h = -1 |
| 1802 | 1 | CS%id_CFL = -1 |
| 1803 | ||
| 1804 | CS%id_KhTr_u = register_diag_field('ocean_model', 'KHTR_u', diag%axesCui, Time, & | |
| 1805 | 1 | 'Epipycnal tracer diffusivity at zonal faces of tracer cell', 'm2 s-1', conversion=US%L_to_m**2*US%s_to_T) |
| 1806 | CS%id_KhTr_v = register_diag_field('ocean_model', 'KHTR_v', diag%axesCvi, Time, & | |
| 1807 | 1 | 'Epipycnal tracer diffusivity at meridional faces of tracer cell', 'm2 s-1', conversion=US%L_to_m**2*US%s_to_T) |
| 1808 | CS%id_KhTr_h = register_diag_field('ocean_model', 'KHTR_h', diag%axesTi, Time, & | |
| 1809 | 'Epipycnal tracer diffusivity at tracer cell center', 'm2 s-1', conversion=US%L_to_m**2*US%s_to_T, & | |
| 1810 | cmor_field_name='diftrelo', & | |
| 1811 | cmor_standard_name= 'ocean_tracer_epineutral_laplacian_diffusivity', & | |
| 1812 | 1 | cmor_long_name = 'Ocean Tracer Epineutral Laplacian Diffusivity') |
| 1813 | ||
| 1814 | CS%id_khdt_x = register_diag_field('ocean_model', 'KHDT_x', diag%axesCu1, Time, & | |
| 1815 | 1 | 'Epipycnal tracer diffusivity operator at zonal faces of tracer cell', 'm2', conversion=US%L_to_m**2) |
| 1816 | CS%id_khdt_y = register_diag_field('ocean_model', 'KHDT_y', diag%axesCv1, Time, & | |
| 1817 | 1 | 'Epipycnal tracer diffusivity operator at meridional faces of tracer cell', 'm2', conversion=US%L_to_m**2) |
| 1818 | 1 | if (CS%check_diffusive_CFL) then |
| 1819 | CS%id_CFL = register_diag_field('ocean_model', 'CFL_lateral_diff', diag%axesT1, Time,& | |
| 1820 | 0 | 'Grid CFL number for lateral/neutral tracer diffusion', 'nondim') |
| 1821 | endif | |
| 1822 | ||
| 1823 | ||
| 1824 | end subroutine tracer_hor_diff_init | |
| 1825 | ||
| 1826 | 1 | subroutine tracer_hor_diff_end(CS) |
| 1827 | type(tracer_hor_diff_CS), pointer :: CS !< module control structure | |
| 1828 | ||
| 1829 | 1 | call neutral_diffusion_end(CS%neutral_diffusion_CSp) |
| 1830 | 1 | call hor_bnd_diffusion_end(CS%hor_bnd_diffusion_CSp) |
| 1831 | 1 | if (associated(CS)) deallocate(CS) |
| 1832 | ||
| 1833 | 1 | end subroutine tracer_hor_diff_end |
| 1834 | ||
| 1835 | ||
| 1836 | !> \namespace mom_tracer_hor_diff | |
| 1837 | !! | |
| 1838 | !! \section section_intro Introduction to the module | |
| 1839 | !! | |
| 1840 | !! This module contains subroutines that handle horizontal | |
| 1841 | !! diffusion (i.e., isoneutral or along layer) of tracers. | |
| 1842 | !! | |
| 1843 | !! Each of the tracers are subject to Fickian along-coordinate | |
| 1844 | !! diffusion if Khtr is defined and positive. The tracer diffusion | |
| 1845 | !! can use a suitable number of iterations to guarantee stability | |
| 1846 | !! with an arbitrarily large time step. | |
| 1847 | ||
| 1848 | 0 | end module MOM_tracer_hor_diff |