portedportable, not yet portedexecuted, not portableexecutable, not hit by this run
| 1 | ! This file is part of MOM6, the Modular Ocean Model version 6. | |
| 2 | ! See the LICENSE file for licensing information. | |
| 3 | ! SPDX-License-Identifier: Apache-2.0 | |
| 4 | ||
| 5 | !> Calculates Zanna and Bolton 2020 parameterization | |
| 6 | !! Implemented by Perezhogin P.A. Contact: pperezhogin@gmail.com | |
| 7 | module MOM_Zanna_Bolton | |
| 8 | ||
| 9 | use MOM_grid, only : ocean_grid_type | |
| 10 | use MOM_verticalGrid, only : verticalGrid_type | |
| 11 | use MOM_diag_mediator, only : diag_ctrl, time_type | |
| 12 | use MOM_file_parser, only : get_param, log_version, param_file_type | |
| 13 | use MOM_unit_scaling, only : unit_scale_type | |
| 14 | use MOM_diag_mediator, only : post_data, register_diag_field | |
| 15 | use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type, & | |
| 16 | start_group_pass, complete_group_pass | |
| 17 | use MOM_domains, only : To_North, To_East | |
| 18 | use MOM_domains, only : pass_var, CORNER | |
| 19 | use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end | |
| 20 | use MOM_cpu_clock, only : CLOCK_MODULE, CLOCK_ROUTINE | |
| 21 | use MOM_ANN, only : ANN_init, ANN_apply_array_sio, ANN_end, ANN_CS | |
| 22 | ||
| 23 | implicit none ; private | |
| 24 | ||
| 25 | #include <MOM_memory.h> | |
| 26 | ||
| 27 | public ZB2020_lateral_stress, ZB2020_init, ZB2020_end, ZB2020_copy_gradient_and_thickness | |
| 28 | ||
| 29 | !> Control structure for Zanna-Bolton-2020 parameterization. | |
| 30 | type, public :: ZB2020_CS ; private | |
| 31 | ! Parameters | |
| 32 | real :: amplitude !< The nondimensional scaling factor in ZB model, | |
| 33 | !! typically 0.1 - 10 [nondim]. | |
| 34 | integer :: ZB_type !< Select how to compute the trace part of ZB model: | |
| 35 | !! 0 - both deviatoric and trace components are computed | |
| 36 | !! 1 - only deviatoric component is computed | |
| 37 | !! 2 - only trace component is computed | |
| 38 | integer :: ZB_cons !< Select a discretization scheme for ZB model | |
| 39 | !! 0 - non-conservative scheme | |
| 40 | !! 1 - conservative scheme for deviatoric component | |
| 41 | integer :: HPF_iter !< Number of sharpening passes for the Velocity Gradient (VG) components | |
| 42 | !! in ZB model. | |
| 43 | integer :: Stress_iter !< Number of smoothing passes for the Stress tensor components | |
| 44 | !! in ZB model. | |
| 45 | real :: Klower_R_diss !< Attenuation of | |
| 46 | !! the ZB parameterization in the regions of | |
| 47 | !! geostrophically-unbalanced flows (Klower 2018, Juricke2020,2019) | |
| 48 | !! Subgrid stress is multiplied by 1/(1+(shear/(f*R_diss))) | |
| 49 | !! R_diss=-1: attenuation is not used; typical value R_diss=1.0 [nondim] | |
| 50 | integer :: Klower_shear !< Type of expression for shear in Klower formula | |
| 51 | !! 0: sqrt(sh_xx**2 + sh_xy**2) | |
| 52 | !! 1: sqrt(sh_xx**2 + sh_xy**2 + vort_xy**2) | |
| 53 | integer :: Marching_halo !< The number of filter iterations per a single MPI | |
| 54 | !! exchange | |
| 55 | ||
| 56 | real, dimension(:,:,:), allocatable :: & | |
| 57 | sh_xx, & !< Horizontal tension (du/dx - dv/dy) in h (CENTER) | |
| 58 | !! points including metric terms [T-1 ~> s-1] | |
| 59 | sh_xy, & !< Horizontal shearing strain (du/dy + dv/dx) in q (CORNER) | |
| 60 | !! points including metric terms [T-1 ~> s-1] | |
| 61 | vort_xy, & !< Vertical vorticity (dv/dx - du/dy) in q (CORNER) | |
| 62 | !! points including metric terms [T-1 ~> s-1] | |
| 63 | hq !< Thickness in CORNER points [H ~> m or kg m-2] | |
| 64 | ||
| 65 | real, dimension(:,:,:), allocatable :: & | |
| 66 | Txx, & !< Subgrid stress xx component in h [L2 T-2 ~> m2 s-2] | |
| 67 | Tyy, & !< Subgrid stress yy component in h [L2 T-2 ~> m2 s-2] | |
| 68 | Txy !< Subgrid stress xy component in q [L2 T-2 ~> m2 s-2] | |
| 69 | ||
| 70 | real, dimension(:,:), allocatable :: & | |
| 71 | kappa_h, & !< Scaling coefficient in h points [L2 ~> m2] | |
| 72 | kappa_q !< Scaling coefficient in q points [L2 ~> m2] | |
| 73 | ||
| 74 | real, allocatable :: & | |
| 75 | ICoriolis_h(:,:), & !< Inverse Coriolis parameter at h points [T ~> s] | |
| 76 | c_diss(:,:,:) !< Attenuation parameter at h points | |
| 77 | !! (Klower 2018, Juricke2019,2020) [nondim] | |
| 78 | ||
| 79 | real, dimension(:,:), allocatable :: & | |
| 80 | maskw_h, & !< Mask of land point at h points multiplied by filter weight [nondim] | |
| 81 | maskw_q !< Same mask but for q points [nondim] | |
| 82 | ||
| 83 | logical :: use_ann !< If True, momentum fluxes are inferred with ANN | |
| 84 | integer :: stencil_size !< Default is 3x3 | |
| 85 | type(ANN_CS) :: ann_Tall !< ANN instance for off-diagonal and diagonal stress | |
| 86 | character(len=200) :: ann_file_Tall !< Path to netcdf file with ANN | |
| 87 | real :: subroundoff_shear !< Small dimensional constant for save division by zero [T-1 ~> s-1] | |
| 88 | ||
| 89 | type(diag_ctrl), pointer :: diag => NULL() !< A type that regulates diagnostics output | |
| 90 | !>@{ Diagnostic handles | |
| 91 | integer :: id_ZB2020u = -1, id_ZB2020v = -1, id_KE_ZB2020 = -1 | |
| 92 | integer :: id_Txx = -1 | |
| 93 | integer :: id_Tyy = -1 | |
| 94 | integer :: id_Txy = -1 | |
| 95 | integer :: id_cdiss = -1 | |
| 96 | !>@} | |
| 97 | ||
| 98 | !>@{ CPU time clock IDs | |
| 99 | integer :: id_clock_module | |
| 100 | integer :: id_clock_copy | |
| 101 | integer :: id_clock_cdiss | |
| 102 | integer :: id_clock_stress | |
| 103 | integer :: id_clock_stress_ANN | |
| 104 | integer :: id_clock_divergence | |
| 105 | integer :: id_clock_mpi | |
| 106 | integer :: id_clock_filter | |
| 107 | integer :: id_clock_post | |
| 108 | integer :: id_clock_source | |
| 109 | !>@} | |
| 110 | ||
| 111 | !>@{ MPI group passes | |
| 112 | type(group_pass_type) :: & | |
| 113 | pass_Tq, pass_Th, & !< handles for halo passes of Txy and Txx, Tyy | |
| 114 | pass_xx, pass_xy !< handles for halo passes of sh_xx and sh_xy, vort_xy | |
| 115 | integer :: Stress_halo = -1, & !< The halo size in filter of the stress tensor | |
| 116 | HPF_halo = -1 !< The halo size in filter of the velocity gradient | |
| 117 | !>@} | |
| 118 | ||
| 119 | end type ZB2020_CS | |
| 120 | ||
| 121 | contains | |
| 122 | ||
| 123 | !> Read parameters, allocate and precompute arrays, | |
| 124 | !! register diagnosicts used in Zanna_Bolton_2020(). | |
| 125 | 1 | subroutine ZB2020_init(Time, G, GV, US, param_file, diag, CS, use_ZB2020) |
| 126 | type(time_type), intent(in) :: Time !< The current model time. | |
| 127 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. | |
| 128 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure | |
| 129 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 130 | type(param_file_type), intent(in) :: param_file !< Parameter file parser structure. | |
| 131 | type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics structure. | |
| 132 | type(ZB2020_CS), intent(inout) :: CS !< ZB2020 control structure. | |
| 133 | logical, intent(out) :: use_ZB2020 !< If true, turns on ZB scheme. | |
| 134 | ||
| 135 | real :: subroundoff_Cor ! A negligible parameter which avoids division by zero | |
| 136 | ! but small compared to Coriolis parameter [T-1 ~> s-1] | |
| 137 | ||
| 138 | integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq | |
| 139 | integer :: i, j | |
| 140 | ||
| 141 | ! This include declares and sets the variable "version". | |
| 142 | #include "version_variable.h" | |
| 143 | character(len=40) :: mdl = "MOM_Zanna_Bolton" ! This module's name. | |
| 144 | ||
| 145 | 1 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec |
| 146 | 1 | Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB |
| 147 | ||
| 148 | 1 | call log_version(param_file, mdl, version, "") |
| 149 | ||
| 150 | call get_param(param_file, mdl, "USE_ZB2020", use_ZB2020, & | |
| 151 | "If true, turns on Zanna-Bolton-2020 (ZB) " //& | |
| 152 | 1 | "subgrid momentum parameterization of mesoscale eddies.", default=.false.) |
| 153 | 1 | if (.not. use_ZB2020) return |
| 154 | ||
| 155 | call get_param(param_file, mdl, "ZB2020_USE_ANN", CS%use_ann, & | |
| 156 | 0 | "ANN inference of momentum fluxes", default=.false.) |
| 157 | ||
| 158 | call get_param(param_file, mdl, "ZB2020_ANN_STENCIL_SIZE", CS%stencil_size, & | |
| 159 | 0 | "ANN stencil size", default=3) |
| 160 | ||
| 161 | call get_param(param_file, mdl, "ZB2020_ANN_FILE_TALL", CS%ann_file_Tall, & | |
| 162 | "ANN parameters for prediction of Txy, Txx and Tyy netcdf input", & | |
| 163 | 0 | default="INPUT/EXP1/Tall.nc") |
| 164 | ||
| 165 | call get_param(param_file, mdl, "ZB_SCALING", CS%amplitude, & | |
| 166 | "The nondimensional scaling factor in ZB model, " //& | |
| 167 | 0 | "typically 0.5-2.5", units="nondim", default=0.5) |
| 168 | ||
| 169 | call get_param(param_file, mdl, "ZB_TRACE_MODE", CS%ZB_type, & | |
| 170 | "Select how to compute the trace part of ZB model:\n" //& | |
| 171 | "\t 0 - both deviatoric and trace components are computed\n" //& | |
| 172 | "\t 1 - only deviatoric component is computed\n" //& | |
| 173 | 0 | "\t 2 - only trace component is computed", default=0) |
| 174 | ||
| 175 | call get_param(param_file, mdl, "ZB_SCHEME", CS%ZB_cons, & | |
| 176 | "Select a discretization scheme for ZB model:\n" //& | |
| 177 | "\t 0 - non-conservative scheme\n" //& | |
| 178 | 0 | "\t 1 - conservative scheme for deviatoric component", default=1) |
| 179 | ||
| 180 | call get_param(param_file, mdl, "VG_SHARP_PASS", CS%HPF_iter, & | |
| 181 | "Number of sharpening passes for the Velocity Gradient (VG) components " //& | |
| 182 | 0 | "in ZB model.", default=0) |
| 183 | ||
| 184 | call get_param(param_file, mdl, "STRESS_SMOOTH_PASS", CS%Stress_iter, & | |
| 185 | "Number of smoothing passes for the Stress tensor components " //& | |
| 186 | 0 | "in ZB model.", default=0) |
| 187 | ||
| 188 | call get_param(param_file, mdl, "ZB_KLOWER_R_DISS", CS%Klower_R_diss, & | |
| 189 | "Attenuation of " //& | |
| 190 | "the ZB parameterization in the regions of " //& | |
| 191 | "geostrophically-unbalanced flows (Klower 2018, Juricke2020,2019). " //& | |
| 192 | "Subgrid stress is multiplied by 1/(1+(shear/(f*R_diss))):\n" //& | |
| 193 | "\t R_diss=-1. - attenuation is not used\n\t R_diss= 1. - typical value", & | |
| 194 | 0 | units="nondim", default=-1.) |
| 195 | ||
| 196 | call get_param(param_file, mdl, "ZB_KLOWER_SHEAR", CS%Klower_shear, & | |
| 197 | "Type of expression for shear in Klower formula:\n" //& | |
| 198 | "\t 0: sqrt(sh_xx**2 + sh_xy**2)\n" //& | |
| 199 | "\t 1: sqrt(sh_xx**2 + sh_xy**2 + vort_xy**2)", & | |
| 200 | 0 | default=1, do_not_log=.not.CS%Klower_R_diss>0) |
| 201 | ||
| 202 | call get_param(param_file, mdl, "ZB_MARCHING_HALO", CS%Marching_halo, & | |
| 203 | "The number of filter iterations per single MPI " //& | |
| 204 | 0 | "exchange", default=4, do_not_log=(CS%Stress_iter==0).and.(CS%HPF_iter==0)) |
| 205 | ||
| 206 | ! Register fields for output from this module. | |
| 207 | 0 | CS%diag => diag |
| 208 | ||
| 209 | CS%id_ZB2020u = register_diag_field('ocean_model', 'ZB2020u', diag%axesCuL, Time, & | |
| 210 | 0 | 'Zonal Acceleration from Zanna-Bolton 2020', 'm s-2', conversion=US%L_T2_to_m_s2) |
| 211 | CS%id_ZB2020v = register_diag_field('ocean_model', 'ZB2020v', diag%axesCvL, Time, & | |
| 212 | 0 | 'Meridional Acceleration from Zanna-Bolton 2020', 'm s-2', conversion=US%L_T2_to_m_s2) |
| 213 | CS%id_KE_ZB2020 = register_diag_field('ocean_model', 'KE_ZB2020', diag%axesTL, Time, & | |
| 214 | 'Kinetic Energy Source from Horizontal Viscosity', & | |
| 215 | 0 | 'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) |
| 216 | ||
| 217 | CS%id_Txx = register_diag_field('ocean_model', 'Txx', diag%axesTL, Time, & | |
| 218 | 0 | 'Diagonal term (Txx) in the ZB stress tensor', 'm2 s-2', conversion=US%L_T_to_m_s**2) |
| 219 | ||
| 220 | CS%id_Tyy = register_diag_field('ocean_model', 'Tyy', diag%axesTL, Time, & | |
| 221 | 0 | 'Diagonal term (Tyy) in the ZB stress tensor', 'm2 s-2', conversion=US%L_T_to_m_s**2) |
| 222 | ||
| 223 | CS%id_Txy = register_diag_field('ocean_model', 'Txy', diag%axesBL, Time, & | |
| 224 | 0 | 'Off-diagonal term (Txy) in the ZB stress tensor', 'm2 s-2', conversion=US%L_T_to_m_s**2) |
| 225 | ||
| 226 | 0 | if (CS%Klower_R_diss > 0) then |
| 227 | CS%id_cdiss = register_diag_field('ocean_model', 'c_diss', diag%axesTL, Time, & | |
| 228 | 0 | 'Klower (2018) attenuation coefficient', 'nondim') |
| 229 | endif | |
| 230 | ||
| 231 | ! Clock IDs | |
| 232 | ! Only module is measured with syncronization. While smaller | |
| 233 | ! parts are measured without - because these are nested clocks. | |
| 234 | 0 | CS%id_clock_module = cpu_clock_id('(Ocean Zanna-Bolton-2020)', grain=CLOCK_MODULE) |
| 235 | 0 | CS%id_clock_copy = cpu_clock_id('(ZB2020 copy fields)', grain=CLOCK_ROUTINE, sync=.false.) |
| 236 | 0 | CS%id_clock_cdiss = cpu_clock_id('(ZB2020 compute c_diss)', grain=CLOCK_ROUTINE, sync=.false.) |
| 237 | 0 | CS%id_clock_stress = cpu_clock_id('(ZB2020 compute stress)', grain=CLOCK_ROUTINE, sync=.false.) |
| 238 | 0 | CS%id_clock_stress_ANN = cpu_clock_id('(ZB2020 compute stress ANN)', grain=CLOCK_ROUTINE, sync=.false.) |
| 239 | 0 | CS%id_clock_divergence = cpu_clock_id('(ZB2020 compute divergence)', grain=CLOCK_ROUTINE, sync=.false.) |
| 240 | 0 | CS%id_clock_mpi = cpu_clock_id('(ZB2020 filter MPI exchanges)', grain=CLOCK_ROUTINE, sync=.false.) |
| 241 | 0 | CS%id_clock_filter = cpu_clock_id('(ZB2020 filter no MPI)', grain=CLOCK_ROUTINE, sync=.false.) |
| 242 | 0 | CS%id_clock_post = cpu_clock_id('(ZB2020 post data)', grain=CLOCK_ROUTINE, sync=.false.) |
| 243 | 0 | CS%id_clock_source = cpu_clock_id('(ZB2020 compute energy source)', grain=CLOCK_ROUTINE, sync=.false.) |
| 244 | ||
| 245 | 0 | CS%subroundoff_shear = 1e-30 * US%T_to_s |
| 246 | 0 | if (CS%use_ann) then |
| 247 | 0 | call ANN_init(CS%ann_Tall, CS%ann_file_Tall) |
| 248 | endif | |
| 249 | ||
| 250 | ! Allocate memory | |
| 251 | ! We set the stress tensor and velocity gradient tensor to zero | |
| 252 | ! with full halo because they potentially may be filtered | |
| 253 | ! with marching halo algorithm | |
| 254 | 0 | allocate(CS%sh_xx(SZI_(G),SZJ_(G),SZK_(GV)), source=0.) |
| 255 | 0 | allocate(CS%sh_xy(SZIB_(G),SZJB_(G),SZK_(GV)), source=0.) |
| 256 | 0 | allocate(CS%vort_xy(SZIB_(G),SZJB_(G),SZK_(GV)), source=0.) |
| 257 | 0 | allocate(CS%hq(SZIB_(G),SZJB_(G),SZK_(GV))) |
| 258 | ||
| 259 | 0 | allocate(CS%Txx(SZI_(G),SZJ_(G),SZK_(GV)), source=0.) |
| 260 | 0 | allocate(CS%Tyy(SZI_(G),SZJ_(G),SZK_(GV)), source=0.) |
| 261 | 0 | allocate(CS%Txy(SZIB_(G),SZJB_(G),SZK_(GV)), source=0.) |
| 262 | 0 | allocate(CS%kappa_h(SZI_(G),SZJ_(G))) |
| 263 | 0 | allocate(CS%kappa_q(SZIB_(G),SZJB_(G))) |
| 264 | ||
| 265 | ! Precomputing the scaling coefficient | |
| 266 | ! Mask is included to automatically satisfy B.C. | |
| 267 | 0 | do j=js-2,je+2 ; do i=is-2,ie+2 |
| 268 | 0 | CS%kappa_h(i,j) = -CS%amplitude * G%areaT(i,j) * G%mask2dT(i,j) |
| 269 | enddo ; enddo | |
| 270 | ||
| 271 | 0 | do J=Jsq-2,Jeq+2 ; do I=Isq-2,Ieq+2 |
| 272 | 0 | CS%kappa_q(I,J) = -CS%amplitude * G%areaBu(I,J) * G%mask2dBu(I,J) |
| 273 | enddo ; enddo | |
| 274 | ||
| 275 | 0 | if (CS%Klower_R_diss > 0) then |
| 276 | 0 | allocate(CS%ICoriolis_h(SZI_(G),SZJ_(G))) |
| 277 | 0 | allocate(CS%c_diss(SZI_(G),SZJ_(G),SZK_(GV))) |
| 278 | ||
| 279 | 0 | subroundoff_Cor = 1e-30 * US%T_to_s |
| 280 | ! Precomputing 1/(f * R_diss) | |
| 281 | 0 | do j=js-1,je+1 ; do i=is-1,ie+1 |
| 282 | CS%ICoriolis_h(i,j) = 1. / ((abs(0.25 * ((G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1)) & | |
| 283 | + (G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1)))) + subroundoff_Cor) & | |
| 284 | 0 | * CS%Klower_R_diss) |
| 285 | enddo ; enddo | |
| 286 | endif | |
| 287 | ||
| 288 | 0 | if (CS%Stress_iter > 0 .or. CS%HPF_iter > 0) then |
| 289 | ! Include 1/16. factor to the mask for filter implementation | |
| 290 | 0 | allocate(CS%maskw_h(SZI_(G),SZJ_(G))) ; CS%maskw_h(:,:) = G%mask2dT(:,:) * 0.0625 |
| 291 | 0 | allocate(CS%maskw_q(SZIB_(G),SZJB_(G))) ; CS%maskw_q(:,:) = G%mask2dBu(:,:) * 0.0625 |
| 292 | endif | |
| 293 | ||
| 294 | ! Initialize MPI group passes | |
| 295 | 0 | if (CS%Stress_iter > 0) then |
| 296 | ! reduce size of halo exchange accordingly to | |
| 297 | ! Marching halo, number of iterations and the array size | |
| 298 | ! But let exchange width be at least 1 | |
| 299 | CS%Stress_halo = max(min(CS%Marching_halo, CS%Stress_iter, & | |
| 300 | 0 | G%Domain%nihalo, G%Domain%njhalo), 1) |
| 301 | ||
| 302 | call create_group_pass(CS%pass_Tq, CS%Txy, G%Domain, halo=CS%Stress_halo, & | |
| 303 | 0 | position=CORNER) |
| 304 | 0 | call create_group_pass(CS%pass_Th, CS%Txx, G%Domain, halo=CS%Stress_halo) |
| 305 | 0 | call create_group_pass(CS%pass_Th, CS%Tyy, G%Domain, halo=CS%Stress_halo) |
| 306 | endif | |
| 307 | ||
| 308 | 0 | if (CS%HPF_iter > 0) then |
| 309 | ! The minimum halo size is 2 because it is requirement for the | |
| 310 | ! outputs of function filter_velocity_gradients | |
| 311 | CS%HPF_halo = max(min(CS%Marching_halo, CS%HPF_iter, & | |
| 312 | 0 | G%Domain%nihalo, G%Domain%njhalo), 2) |
| 313 | ||
| 314 | 0 | call create_group_pass(CS%pass_xx, CS%sh_xx, G%Domain, halo=CS%HPF_halo) |
| 315 | call create_group_pass(CS%pass_xy, CS%sh_xy, G%Domain, halo=CS%HPF_halo, & | |
| 316 | 0 | position=CORNER) |
| 317 | call create_group_pass(CS%pass_xy, CS%vort_xy, G%Domain, halo=CS%HPF_halo, & | |
| 318 | 0 | position=CORNER) |
| 319 | endif | |
| 320 | ||
| 321 | end subroutine ZB2020_init | |
| 322 | ||
| 323 | !> Deallocate any variables allocated in ZB_2020_init | |
| 324 | 0 | subroutine ZB2020_end(CS) |
| 325 | type(ZB2020_CS), intent(inout) :: CS !< ZB2020 control structure. | |
| 326 | ||
| 327 | 0 | deallocate(CS%sh_xx) |
| 328 | 0 | deallocate(CS%sh_xy) |
| 329 | 0 | deallocate(CS%vort_xy) |
| 330 | 0 | deallocate(CS%hq) |
| 331 | ||
| 332 | 0 | deallocate(CS%Txx) |
| 333 | 0 | deallocate(CS%Tyy) |
| 334 | 0 | deallocate(CS%Txy) |
| 335 | 0 | deallocate(CS%kappa_h) |
| 336 | 0 | deallocate(CS%kappa_q) |
| 337 | ||
| 338 | 0 | if (CS%Klower_R_diss > 0) then |
| 339 | 0 | deallocate(CS%ICoriolis_h) |
| 340 | 0 | deallocate(CS%c_diss) |
| 341 | endif | |
| 342 | ||
| 343 | 0 | if (CS%Stress_iter > 0 .or. CS%HPF_iter > 0) then |
| 344 | 0 | deallocate(CS%maskw_h) |
| 345 | 0 | deallocate(CS%maskw_q) |
| 346 | endif | |
| 347 | ||
| 348 | 0 | if (CS%use_ann) then |
| 349 | 0 | call ANN_end(CS%ann_Tall) |
| 350 | endif | |
| 351 | ||
| 352 | 0 | end subroutine ZB2020_end |
| 353 | ||
| 354 | !> Save precomputed velocity gradients and thickness | |
| 355 | !! from the horizontal eddy viscosity module | |
| 356 | !! We save as much halo for velocity gradients as possible | |
| 357 | !! In symmetric (preferable) memory model: halo 2 for sh_xx | |
| 358 | !! and halo 1 for sh_xy and vort_xy | |
| 359 | !! We apply zero boundary conditions to velocity gradients | |
| 360 | !! which is required for filtering operations | |
| 361 | 0 | subroutine ZB2020_copy_gradient_and_thickness(sh_xx, sh_xy, vort_xy, hq, & |
| 362 | G, GV, CS, k) | |
| 363 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. | |
| 364 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. | |
| 365 | type(ZB2020_CS), intent(inout) :: CS !< ZB2020 control structure. | |
| 366 | ||
| 367 | real, dimension(SZIB_(G),SZJB_(G)), & | |
| 368 | intent(in) :: sh_xy !< horizontal shearing strain (du/dy + dv/dx) | |
| 369 | !! including metric terms [T-1 ~> s-1] | |
| 370 | real, dimension(SZIB_(G),SZJB_(G)), & | |
| 371 | intent(in) :: vort_xy !< Vertical vorticity (dv/dx - du/dy) | |
| 372 | !! including metric terms [T-1 ~> s-1] | |
| 373 | real, dimension(SZIB_(G),SZJB_(G)), & | |
| 374 | intent(in) :: hq !< harmonic mean of the harmonic means | |
| 375 | !! of the u- & v point thicknesses [H ~> m or kg m-2] | |
| 376 | ||
| 377 | real, dimension(SZI_(G),SZJ_(G)), & | |
| 378 | intent(in) :: sh_xx !< horizontal tension (du/dx - dv/dy) | |
| 379 | !! including metric terms [T-1 ~> s-1] | |
| 380 | ||
| 381 | integer, intent(in) :: k !< The vertical index of the layer to be passed. | |
| 382 | ||
| 383 | integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq | |
| 384 | integer :: i, j | |
| 385 | ||
| 386 | 0 | call cpu_clock_begin(CS%id_clock_copy) |
| 387 | ||
| 388 | 0 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec |
| 389 | 0 | Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB |
| 390 | ||
| 391 | 0 | do J=js-1,Jeq ; do I=is-1,Ieq |
| 392 | 0 | CS%hq(I,J,k) = hq(I,J) |
| 393 | enddo ; enddo | |
| 394 | ||
| 395 | ! No physical B.C. is required for | |
| 396 | ! sh_xx in ZB2020. However, filtering | |
| 397 | ! may require BC | |
| 398 | 0 | do j=Jsq-1,je+2 ; do i=Isq-1,ie+2 |
| 399 | 0 | CS%sh_xx(i,j,k) = sh_xx(i,j) * G%mask2dT(i,j) |
| 400 | enddo ; enddo | |
| 401 | ||
| 402 | ! We multiply by mask to remove | |
| 403 | ! implicit dependence on CS%no_slip | |
| 404 | ! flag in hor_visc module | |
| 405 | 0 | do J=js-2,Jeq+1 ; do I=is-2,Ieq+1 |
| 406 | 0 | CS%sh_xy(I,J,k) = sh_xy(I,J) * G%mask2dBu(I,J) |
| 407 | enddo ; enddo | |
| 408 | ||
| 409 | 0 | do J=js-2,Jeq+1 ; do I=is-2,Ieq+1 |
| 410 | 0 | CS%vort_xy(I,J,k) = vort_xy(I,J) * G%mask2dBu(I,J) |
| 411 | enddo ; enddo | |
| 412 | ||
| 413 | 0 | call cpu_clock_end(CS%id_clock_copy) |
| 414 | ||
| 415 | 0 | end subroutine ZB2020_copy_gradient_and_thickness |
| 416 | ||
| 417 | !> Baroclinic Zanna-Bolton-2020 parameterization, see | |
| 418 | !! eq. 6 in https://laurezanna.github.io/files/Zanna-Bolton-2020.pdf | |
| 419 | !! We compute the lateral stress tensor according to ZB2020 model | |
| 420 | !! and update the acceleration due to eddy viscosity (diffu, diffv) | |
| 421 | !! as follows: | |
| 422 | !! diffu = diffu + ZB2020u | |
| 423 | !! diffv = diffv + ZB2020v | |
| 424 | 0 | subroutine ZB2020_lateral_stress(u, v, h, diffu, diffv, G, GV, CS, & |
| 425 | 0 | dx2h, dy2h, dx2q, dy2q) |
| 426 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. | |
| 427 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. | |
| 428 | type(ZB2020_CS), intent(inout) :: CS !< ZB2020 control structure. | |
| 429 | ||
| 430 | real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & | |
| 431 | intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1]. | |
| 432 | real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & | |
| 433 | intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1]. | |
| 434 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & | |
| 435 | intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]. | |
| 436 | ||
| 437 | real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & | |
| 438 | intent(inout) :: diffu !< Zonal acceleration due to eddy viscosity. | |
| 439 | !! It is updated with ZB closure [L T-2 ~> m s-2] | |
| 440 | real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & | |
| 441 | intent(inout) :: diffv !< Meridional acceleration due to eddy viscosity. | |
| 442 | !! It is updated with ZB closure [L T-2 ~> m s-2] | |
| 443 | ||
| 444 | real, dimension(SZI_(G),SZJ_(G)), intent(in) :: dx2h !< dx^2 at h points [L2 ~> m2] | |
| 445 | real, dimension(SZI_(G),SZJ_(G)), intent(in) :: dy2h !< dy^2 at h points [L2 ~> m2] | |
| 446 | ||
| 447 | real, dimension(SZIB_(G),SZJB_(G)), intent(in) :: dx2q !< dx^2 at q points [L2 ~> m2] | |
| 448 | real, dimension(SZIB_(G),SZJB_(G)), intent(in) :: dy2q !< dy^2 at q points [L2 ~> m2] | |
| 449 | ||
| 450 | 0 | call cpu_clock_begin(CS%id_clock_module) |
| 451 | ||
| 452 | ! Compute attenuation if specified | |
| 453 | 0 | call compute_c_diss(G, GV, CS) |
| 454 | ||
| 455 | ! Sharpen velocity gradients if specified | |
| 456 | 0 | call filter_velocity_gradients(G, GV, CS) |
| 457 | ||
| 458 | ! Compute the stress tensor given the | |
| 459 | ! (optionally sharpened) velocity gradients | |
| 460 | 0 | if (CS%use_ann) then |
| 461 | 0 | call compute_stress_ANN_collocated(G, GV, CS) |
| 462 | else | |
| 463 | 0 | call compute_stress(G, GV, CS) |
| 464 | endif | |
| 465 | ||
| 466 | ! Smooth the stress tensor if specified | |
| 467 | 0 | call filter_stress(G, GV, CS) |
| 468 | ||
| 469 | ! Update the acceleration due to eddy viscosity (diffu, diffv) | |
| 470 | ! with the ZB2020 lateral parameterization | |
| 471 | call compute_stress_divergence(u, v, h, diffu, diffv, & | |
| 472 | dx2h, dy2h, dx2q, dy2q, & | |
| 473 | 0 | G, GV, CS) |
| 474 | ||
| 475 | 0 | call cpu_clock_begin(CS%id_clock_post) |
| 476 | 0 | if (CS%id_Txx>0) call post_data(CS%id_Txx, CS%Txx, CS%diag) |
| 477 | 0 | if (CS%id_Tyy>0) call post_data(CS%id_Tyy, CS%Tyy, CS%diag) |
| 478 | 0 | if (CS%id_Txy>0) call post_data(CS%id_Txy, CS%Txy, CS%diag) |
| 479 | ||
| 480 | 0 | if (CS%id_cdiss>0) call post_data(CS%id_cdiss, CS%c_diss, CS%diag) |
| 481 | 0 | call cpu_clock_end(CS%id_clock_post) |
| 482 | ||
| 483 | 0 | call cpu_clock_end(CS%id_clock_module) |
| 484 | ||
| 485 | 0 | end subroutine ZB2020_lateral_stress |
| 486 | ||
| 487 | !> Compute the attenuation parameter similarly | |
| 488 | !! to Klower2018, Juricke2019,2020: c_diss = 1/(1+(shear/(f*R_diss))) | |
| 489 | !! where shear = sqrt(sh_xx**2 + sh_xy**2) or shear = sqrt(sh_xx**2 + sh_xy**2 + vort_xy**2) | |
| 490 | !! In symmetric memory model, components of velocity gradient tensor | |
| 491 | !! should have halo 1 and zero boundary conditions. The result: c_diss having halo 1. | |
| 492 | 0 | subroutine compute_c_diss(G, GV, CS) |
| 493 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. | |
| 494 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure | |
| 495 | type(ZB2020_CS), intent(inout) :: CS !< ZB2020 control structure. | |
| 496 | ||
| 497 | integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz | |
| 498 | integer :: i, j, k | |
| 499 | ||
| 500 | real :: shear ! Shear in Klower2018 formula at h points [T-1 ~> s-1] | |
| 501 | ||
| 502 | 0 | if (.not. CS%Klower_R_diss > 0) & |
| 503 | 0 | return |
| 504 | ||
| 505 | 0 | call cpu_clock_begin(CS%id_clock_cdiss) |
| 506 | ||
| 507 | 0 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke |
| 508 | 0 | Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB |
| 509 | ||
| 510 | 0 | do k=1,nz |
| 511 | ||
| 512 | ! sqrt(sh_xx**2 + sh_xy**2) | |
| 513 | 0 | if (CS%Klower_shear == 0) then |
| 514 | 0 | do j=js-1,je+1 ; do i=is-1,ie+1 |
| 515 | shear = sqrt(CS%sh_xx(i,j,k)**2 + 0.25 * ( & | |
| 516 | ((CS%sh_xy(I-1,J-1,k)**2) + (CS%sh_xy(I,J ,k)**2)) & | |
| 517 | + ((CS%sh_xy(I-1,J ,k)**2) + (CS%sh_xy(I,J-1,k)**2)) & | |
| 518 | 0 | )) |
| 519 | 0 | CS%c_diss(i,j,k) = 1. / (1. + shear * CS%ICoriolis_h(i,j)) |
| 520 | enddo ; enddo | |
| 521 | ||
| 522 | ! sqrt(sh_xx**2 + sh_xy**2 + vort_xy**2) | |
| 523 | 0 | elseif (CS%Klower_shear == 1) then |
| 524 | 0 | do j=js-1,je+1 ; do i=is-1,ie+1 |
| 525 | shear = sqrt(CS%sh_xx(i,j,k)**2 + 0.25 * ( & | |
| 526 | ((CS%sh_xy(I-1,J-1,k)**2 + CS%vort_xy(I-1,J-1,k)**2) & | |
| 527 | + (CS%sh_xy(I,J,k)**2 + CS%vort_xy(I,J,k)**2)) & | |
| 528 | + ((CS%sh_xy(I-1,J,k)**2 + CS%vort_xy(I-1,J,k)**2) & | |
| 529 | + (CS%sh_xy(I,J-1,k)**2 + CS%vort_xy(I,J-1,k)**2)) & | |
| 530 | 0 | )) |
| 531 | 0 | CS%c_diss(i,j,k) = 1. / (1. + shear * CS%ICoriolis_h(i,j)) |
| 532 | enddo ; enddo | |
| 533 | endif | |
| 534 | ||
| 535 | enddo ! end of k loop | |
| 536 | ||
| 537 | 0 | call cpu_clock_end(CS%id_clock_cdiss) |
| 538 | ||
| 539 | end subroutine compute_c_diss | |
| 540 | ||
| 541 | !> Compute stress tensor T = | |
| 542 | !! (Txx, Txy; | |
| 543 | !! Txy, Tyy) | |
| 544 | !! Which consists of the deviatoric and trace components, respectively: | |
| 545 | !! T = (-vort_xy * sh_xy, vort_xy * sh_xx; | |
| 546 | !! vort_xy * sh_xx, vort_xy * sh_xy) + | |
| 547 | !! 1/2 * (vort_xy^2 + sh_xy^2 + sh_xx^2, 0; | |
| 548 | !! 0, vort_xy^2 + sh_xy^2 + sh_xx^2) | |
| 549 | !! This stress tensor is multiplied by precomputed kappa=-CS%amplitude * G%area: | |
| 550 | !! T -> T * kappa | |
| 551 | !! The sign of the stress tensor is such that (neglecting h): | |
| 552 | !! (du/dt, dv/dt) = div(T) | |
| 553 | !! In symmetric memory model: sh_xy and vort_xy should have halo 1 | |
| 554 | !! and zero B.C.; sh_xx should have halo 2 and zero B.C. | |
| 555 | !! Result: Txx, Tyy, Txy with halo 1 and zero B.C. | |
| 556 | 0 | subroutine compute_stress(G, GV, CS) |
| 557 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. | |
| 558 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure | |
| 559 | type(ZB2020_CS), intent(inout) :: CS !< ZB2020 control structure. | |
| 560 | ||
| 561 | real :: & | |
| 562 | vort_xy_h, & ! Vorticity interpolated to h point [T-1 ~> s-1] | |
| 563 | sh_xy_h ! Shearing strain interpolated to h point [T-1 ~> s-1] | |
| 564 | ||
| 565 | real :: & | |
| 566 | sh_xx_q ! Horizontal tension interpolated to q point [T-1 ~> s-1] | |
| 567 | ||
| 568 | ! Local variables | |
| 569 | real :: sum_sq ! 1/2*(vort_xy^2 + sh_xy^2 + sh_xx^2) in h point [T-2 ~> s-2] | |
| 570 | real :: vort_sh ! vort_xy*sh_xy in h point [T-2 ~> s-2] | |
| 571 | ||
| 572 | integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz | |
| 573 | integer :: i, j, k | |
| 574 | ||
| 575 | logical :: sum_sq_flag ! Flag to compute trace | |
| 576 | logical :: vort_sh_scheme_0, vort_sh_scheme_1 ! Flags to compute diagonal trace-free part | |
| 577 | ||
| 578 | 0 | call cpu_clock_begin(CS%id_clock_stress) |
| 579 | ||
| 580 | 0 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke |
| 581 | 0 | Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB |
| 582 | ||
| 583 | 0 | sum_sq = 0. |
| 584 | 0 | vort_sh = 0. |
| 585 | ||
| 586 | 0 | sum_sq_flag = CS%ZB_type /= 1 |
| 587 | 0 | vort_sh_scheme_0 = CS%ZB_type /= 2 .and. CS%ZB_cons == 0 |
| 588 | 0 | vort_sh_scheme_1 = CS%ZB_type /= 2 .and. CS%ZB_cons == 1 |
| 589 | ||
| 590 | 0 | do k=1,nz |
| 591 | ||
| 592 | ! compute Txx, Tyy tensor | |
| 593 | 0 | do j=js-1,je+1 ; do i=is-1,ie+1 |
| 594 | ! It is assumed that B.C. is applied to sh_xy and vort_xy | |
| 595 | sh_xy_h = 0.25 * ( (CS%sh_xy(I-1,J-1,k) + CS%sh_xy(I,J,k)) & | |
| 596 | 0 | + (CS%sh_xy(I-1,J,k) + CS%sh_xy(I,J-1,k)) ) |
| 597 | ||
| 598 | vort_xy_h = 0.25 * ( (CS%vort_xy(I-1,J-1,k) + CS%vort_xy(I,J,k)) & | |
| 599 | 0 | + (CS%vort_xy(I-1,J,k) + CS%vort_xy(I,J-1,k)) ) |
| 600 | ||
| 601 | 0 | if (sum_sq_flag) then |
| 602 | sum_sq = 0.5 * & | |
| 603 | ((vort_xy_h * vort_xy_h & | |
| 604 | + sh_xy_h * sh_xy_h) & | |
| 605 | + CS%sh_xx(i,j,k) * CS%sh_xx(i,j,k) & | |
| 606 | 0 | ) |
| 607 | endif | |
| 608 | ||
| 609 | 0 | if (vort_sh_scheme_0) & |
| 610 | 0 | vort_sh = vort_xy_h * sh_xy_h |
| 611 | ||
| 612 | 0 | if (vort_sh_scheme_1) then |
| 613 | ! It is assumed that B.C. is applied to sh_xy and vort_xy | |
| 614 | vort_sh = 0.25 * ( & | |
| 615 | (((G%areaBu(I-1,J-1) * CS%vort_xy(I-1,J-1,k)) * CS%sh_xy(I-1,J-1,k)) + & | |
| 616 | ((G%areaBu(I ,J ) * CS%vort_xy(I ,J ,k)) * CS%sh_xy(I ,J ,k))) + & | |
| 617 | (((G%areaBu(I-1,J ) * CS%vort_xy(I-1,J ,k)) * CS%sh_xy(I-1,J ,k)) + & | |
| 618 | ((G%areaBu(I ,J-1) * CS%vort_xy(I ,J-1,k)) * CS%sh_xy(I ,J-1,k))) & | |
| 619 | 0 | ) * G%IareaT(i,j) |
| 620 | endif | |
| 621 | ||
| 622 | ! B.C. is already applied in kappa_h | |
| 623 | 0 | CS%Txx(i,j,k) = CS%kappa_h(i,j) * (- vort_sh + sum_sq) |
| 624 | 0 | CS%Tyy(i,j,k) = CS%kappa_h(i,j) * (+ vort_sh + sum_sq) |
| 625 | ||
| 626 | enddo ; enddo | |
| 627 | ||
| 628 | ! Here we assume that Txy is initialized to zero | |
| 629 | 0 | if (CS%ZB_type /= 2) then |
| 630 | 0 | do J=Jsq-1,Jeq+1 ; do I=Isq-1,Ieq+1 |
| 631 | sh_xx_q = 0.25 * ( (CS%sh_xx(i+1,j+1,k) + CS%sh_xx(i,j,k)) & | |
| 632 | 0 | + (CS%sh_xx(i+1,j,k) + CS%sh_xx(i,j+1,k))) |
| 633 | ! B.C. is already applied in kappa_q | |
| 634 | 0 | CS%Txy(I,J,k) = CS%kappa_q(I,J) * (CS%vort_xy(I,J,k) * sh_xx_q) |
| 635 | ||
| 636 | enddo ; enddo | |
| 637 | endif | |
| 638 | ||
| 639 | enddo ! end of k loop | |
| 640 | ||
| 641 | 0 | call cpu_clock_end(CS%id_clock_stress) |
| 642 | ||
| 643 | 0 | end subroutine compute_stress |
| 644 | ||
| 645 | !> Compute stress tensor T = | |
| 646 | !! (Txx, Txy; | |
| 647 | !! Txy, Tyy) | |
| 648 | !! with ANN in non-dimensional form: | |
| 649 | !! T = dx^2 * |grad V|^2 * ANN(grad V / |grad V|) | |
| 650 | !! The sign of the stress tensor is such that: | |
| 651 | !! (du/dt, dv/dt) = 1/h * div(h * T) | |
| 652 | !! Algorithm: | |
| 653 | !! 1) Interpolate input features (sh_xy, sh_xx, vort_xy) to grid centers | |
| 654 | !! 2) Compute norm of velocity gradients on a stencil | |
| 655 | !! 3) Non-dimensionalize input features | |
| 656 | !! 4) Make ANN inference in grid centers | |
| 657 | !! 5) Restore physical dimensionality and interpolate Txy back to corners | |
| 658 | 0 | subroutine compute_stress_ANN_collocated(G, GV, CS) |
| 659 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. | |
| 660 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure | |
| 661 | type(ZB2020_CS), intent(inout) :: CS !< ZB2020 control structure. | |
| 662 | ||
| 663 | integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz | |
| 664 | integer :: i, j, k, m | |
| 665 | integer :: ii, jj | |
| 666 | integer :: nij | |
| 667 | ||
| 668 | 0 | real, allocatable :: x(:,:) ! Vector of non-dimensional input features |
| 669 | ! number of horizontal grid points x | |
| 670 | ! (sh_xy, sh_xx, vort_xy) on a stencil [nondim] | |
| 671 | 0 | real, allocatable :: y(:,:) ! Vector of nondimensional |
| 672 | ! output features number of horizontal grid points x | |
| 673 | ! (Txy,Txx,Tyy) [nondim] | |
| 674 | real :: yy(3) ! Vector of dimensional | |
| 675 | ! output features (Txy,Txx,Tyy) [L2 T-2 ~> m2 s-2] | |
| 676 | real :: tmp ! Temporal value of squared norm [T-2 ~> s-2] | |
| 677 | integer :: offset ! Half the stencil size. Used for selection | |
| 678 | integer :: stencil_points ! The number of points after flattening | |
| 679 | ||
| 680 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: & | |
| 681 | 0 | sh_xy_h, & ! sh_xy interpolated to the center [T-1 ~> s-1] |
| 682 | 0 | vort_xy_h, & ! vort_xy interpolated to the center [T-1 ~> s-1] |
| 683 | 0 | norm_h ! Norm of input feautres in center points [T-1 ~> s-1] |
| 684 | ||
| 685 | real, dimension(SZI_(G),SZJ_(G)) :: & | |
| 686 | 0 | sqr_h, & ! Squared norm of velocity gradients in center points [T-2 ~> s-2] |
| 687 | 0 | Txy ! Predicted Txy in center points [T-1 ~> s-1] |
| 688 | ||
| 689 | 0 | call cpu_clock_begin(CS%id_clock_stress_ANN) |
| 690 | ||
| 691 | 0 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke |
| 692 | 0 | Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB |
| 693 | ||
| 694 | ! Number of horizontal grid points in ANN inference loop below | |
| 695 | 0 | nij = (ie - is + 5) * (je - js + 5) |
| 696 | 0 | allocate(x(nij, 3 * CS%stencil_size**2)) |
| 697 | 0 | allocate(y(nij, 3)) |
| 698 | ||
| 699 | 0 | sh_xy_h = 0. |
| 700 | 0 | vort_xy_h = 0. |
| 701 | 0 | norm_h = 0. |
| 702 | ||
| 703 | 0 | call pass_var(CS%sh_xy, G%Domain, clock=CS%id_clock_mpi, position=CORNER) |
| 704 | 0 | call pass_var(CS%sh_xx, G%Domain, clock=CS%id_clock_mpi) |
| 705 | 0 | call pass_var(CS%vort_xy, G%Domain, clock=CS%id_clock_mpi, position=CORNER) |
| 706 | ||
| 707 | 0 | offset = (CS%stencil_size-1)/2 |
| 708 | 0 | stencil_points = CS%stencil_size**2 |
| 709 | ||
| 710 | ! Interpolate input features | |
| 711 | 0 | do k=1,nz |
| 712 | 0 | do j=js-2,je+2 ; do i=is-2,ie+2 |
| 713 | ! It is assumed that B.C. is applied to sh_xy and vort_xy | |
| 714 | sh_xy_h(i,j,k) = 0.25 * ( (CS%sh_xy(I-1,J-1,k) + CS%sh_xy(I,J,k)) & | |
| 715 | 0 | + (CS%sh_xy(I-1,J,k) + CS%sh_xy(I,J-1,k)) ) |
| 716 | ||
| 717 | vort_xy_h(i,j,k) = 0.25 * ( (CS%vort_xy(I-1,J-1,k) + CS%vort_xy(I,J,k)) & | |
| 718 | 0 | + (CS%vort_xy(I-1,J,k) + CS%vort_xy(I,J-1,k)) ) |
| 719 | ||
| 720 | 0 | sqr_h(i,j) = (((CS%sh_xx(i,j,k)**2) + (sh_xy_h(i,j,k)**2)) + (vort_xy_h(i,j,k)**2)) * G%mask2dT(i,j) |
| 721 | enddo ; enddo | |
| 722 | ||
| 723 | 0 | do j=js,je ; do i=is,ie |
| 724 | 0 | tmp = 0.0 |
| 725 | 0 | do jj=j-offset,j+offset ; do ii=i-offset,i+offset |
| 726 | 0 | tmp = tmp + sqr_h(ii,jj) |
| 727 | enddo ; enddo | |
| 728 | 0 | norm_h(i,j,k) = sqrt(tmp) |
| 729 | enddo ; enddo | |
| 730 | enddo | |
| 731 | ||
| 732 | 0 | call pass_var(sh_xy_h, G%Domain, clock=CS%id_clock_mpi) |
| 733 | 0 | call pass_var(vort_xy_h, G%Domain, clock=CS%id_clock_mpi) |
| 734 | 0 | call pass_var(norm_h, G%Domain, clock=CS%id_clock_mpi) |
| 735 | ||
| 736 | 0 | do k=1,nz |
| 737 | 0 | m = 0 |
| 738 | 0 | do j=js-2,je+2 ; do i=is-2,ie+2 |
| 739 | 0 | m = m + 1 |
| 740 | x(m,1:stencil_points) = & | |
| 741 | RESHAPE(sh_xy_h(i-offset:i+offset, & | |
| 742 | 0 | j-offset:j+offset,k), (/stencil_points/)) |
| 743 | x(m,stencil_points+1:2*stencil_points) = & | |
| 744 | RESHAPE(CS%sh_xx(i-offset:i+offset, & | |
| 745 | 0 | j-offset:j+offset,k), (/stencil_points/)) |
| 746 | x(m,2*stencil_points+1:3*stencil_points) = & | |
| 747 | RESHAPE(vort_xy_h(i-offset:i+offset, & | |
| 748 | 0 | j-offset:j+offset,k), (/stencil_points/)) |
| 749 | ||
| 750 | 0 | x(m,:) = x(m,:) / (norm_h(i,j,k) + CS%subroundoff_shear) |
| 751 | enddo ; enddo | |
| 752 | ||
| 753 | 0 | call ANN_apply_array_sio(nij, x, y, CS%ann_Tall) |
| 754 | ||
| 755 | 0 | m = 0 |
| 756 | 0 | do j=js-2,je+2 ; do i=is-2,ie+2 |
| 757 | 0 | m = m+1 |
| 758 | 0 | yy(:) = y(m, :) * norm_h(i,j,k) * norm_h(i,j,k) * CS%kappa_h(i,j) |
| 759 | ||
| 760 | 0 | Txy(i,j) = yy(1) |
| 761 | 0 | CS%Txx(i,j,k) = yy(2) |
| 762 | 0 | CS%Tyy(i,j,k) = yy(3) |
| 763 | enddo ; enddo | |
| 764 | ||
| 765 | 0 | do J=Jsq-1,Jeq+1 ; do I=Isq-1,Ieq+1 |
| 766 | CS%Txy(I,J,k) = 0.25 * ( (Txy(i+1,j+1) + Txy(i,j)) & | |
| 767 | 0 | + (Txy(i+1,j) + Txy(i,j+1))) * G%mask2dBu(I,J) |
| 768 | enddo ; enddo | |
| 769 | ||
| 770 | enddo ! end of k loop | |
| 771 | ||
| 772 | 0 | call pass_var(CS%Txy, G%Domain, clock=CS%id_clock_mpi, position=CORNER) |
| 773 | 0 | call pass_var(CS%Txx, G%Domain, clock=CS%id_clock_mpi) |
| 774 | 0 | call pass_var(CS%Tyy, G%Domain, clock=CS%id_clock_mpi) |
| 775 | ||
| 776 | 0 | deallocate(x) |
| 777 | 0 | deallocate(y) |
| 778 | ||
| 779 | 0 | call cpu_clock_end(CS%id_clock_stress_ANN) |
| 780 | ||
| 781 | 0 | end subroutine compute_stress_ANN_collocated |
| 782 | ||
| 783 | !> Compute the divergence of subgrid stress | |
| 784 | !! weighted with thickness, i.e. | |
| 785 | !! (fx,fy) = 1/h Div(h * [Txx, Txy; Txy, Tyy]) | |
| 786 | !! and update the acceleration due to eddy viscosity as | |
| 787 | !! diffu = diffu + dx; diffv = diffv + dy | |
| 788 | !! Optionally, before computing the divergence, we attenuate the stress | |
| 789 | !! according to the Klower formula. | |
| 790 | !! In symmetric memory model: Txx, Tyy, Txy, c_diss should have halo 1 | |
| 791 | !! with applied zero B.C. | |
| 792 | 0 | subroutine compute_stress_divergence(u, v, h, diffu, diffv, dx2h, dy2h, dx2q, dy2q, G, GV, CS) |
| 793 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. | |
| 794 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure | |
| 795 | type(ZB2020_CS), intent(in) :: CS !< ZB2020 control structure. | |
| 796 | real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & | |
| 797 | intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1]. | |
| 798 | real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & | |
| 799 | intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1]. | |
| 800 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & | |
| 801 | intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]. | |
| 802 | real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & | |
| 803 | intent(out) :: diffu !< Zonal acceleration due to convergence of | |
| 804 | !! along-coordinate stress tensor [L T-2 ~> m s-2] | |
| 805 | real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & | |
| 806 | intent(out) :: diffv !< Meridional acceleration due to convergence | |
| 807 | !! of along-coordinate stress tensor [L T-2 ~> m s-2] | |
| 808 | real, dimension(SZI_(G),SZJ_(G)), & | |
| 809 | intent(in) :: dx2h !< dx^2 at h points [L2 ~> m2] | |
| 810 | real, dimension(SZI_(G),SZJ_(G)), & | |
| 811 | intent(in) :: dy2h !< dy^2 at h points [L2 ~> m2] | |
| 812 | real, dimension(SZIB_(G),SZJB_(G)), & | |
| 813 | intent(in) :: dx2q !< dx^2 at q points [L2 ~> m2] | |
| 814 | real, dimension(SZIB_(G),SZJB_(G)), & | |
| 815 | intent(in) :: dy2q !< dy^2 at q points [L2 ~> m2] | |
| 816 | ||
| 817 | ! Local variables | |
| 818 | real, dimension(SZI_(G),SZJ_(G)) :: & | |
| 819 | 0 | Mxx, & ! Subgrid stress Txx multiplied by thickness and dy^2 [H L4 T-2 ~> m5 s-2] |
| 820 | 0 | Myy ! Subgrid stress Tyy multiplied by thickness and dx^2 [H L4 T-2 ~> m5 s-2] |
| 821 | ||
| 822 | real, dimension(SZIB_(G),SZJB_(G)) :: & | |
| 823 | 0 | Mxy ! Subgrid stress Txy multiplied by thickness [H L2 T-2 ~> m3 s-2] |
| 824 | real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: & | |
| 825 | 0 | ZB2020u !< Zonal acceleration due to convergence of |
| 826 | !! along-coordinate stress tensor for ZB model | |
| 827 | !! [L T-2 ~> m s-2] | |
| 828 | real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: & | |
| 829 | 0 | ZB2020v !< Meridional acceleration due to convergence |
| 830 | !! of along-coordinate stress tensor for ZB model | |
| 831 | !! [L T-2 ~> m s-2] | |
| 832 | ||
| 833 | real :: h_u ! Thickness interpolated to u points [H ~> m or kg m-2]. | |
| 834 | real :: h_v ! Thickness interpolated to v points [H ~> m or kg m-2]. | |
| 835 | real :: fx ! Zonal acceleration [L T-2 ~> m s-2] | |
| 836 | real :: fy ! Meridional acceleration [L T-2 ~> m s-2] | |
| 837 | ||
| 838 | real :: h_neglect ! Thickness so small it can be lost in | |
| 839 | ! roundoff and so neglected [H ~> m or kg m-2] | |
| 840 | ||
| 841 | integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz | |
| 842 | integer :: i, j, k | |
| 843 | logical :: save_ZB2020u, save_ZB2020v ! Save the acceleration due to ZB2020 model | |
| 844 | ||
| 845 | 0 | call cpu_clock_begin(CS%id_clock_divergence) |
| 846 | ||
| 847 | 0 | save_ZB2020u = (CS%id_ZB2020u > 0) .or. (CS%id_KE_ZB2020 > 0) |
| 848 | 0 | save_ZB2020v = (CS%id_ZB2020v > 0) .or. (CS%id_KE_ZB2020 > 0) |
| 849 | ||
| 850 | 0 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke |
| 851 | 0 | Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB |
| 852 | ||
| 853 | 0 | h_neglect = GV%H_subroundoff |
| 854 | ||
| 855 | 0 | do k=1,nz |
| 856 | 0 | if (CS%Klower_R_diss > 0) then |
| 857 | 0 | do J=js-1,Jeq ; do I=is-1,Ieq |
| 858 | Mxy(I,J) = (CS%Txy(I,J,k) * & | |
| 859 | (0.25 * ( (CS%c_diss(i,j ,k) + CS%c_diss(i+1,j+1,k)) & | |
| 860 | + (CS%c_diss(i,j+1,k) + CS%c_diss(i+1,j ,k))) & | |
| 861 | ) & | |
| 862 | 0 | ) * CS%hq(I,J,k) |
| 863 | enddo ; enddo | |
| 864 | else | |
| 865 | 0 | do J=js-1,Jeq ; do I=is-1,Ieq |
| 866 | 0 | Mxy(I,J) = CS%Txy(I,J,k) * CS%hq(I,J,k) |
| 867 | enddo ; enddo | |
| 868 | endif | |
| 869 | ||
| 870 | 0 | if (CS%Klower_R_diss > 0) then |
| 871 | 0 | do j=js-1,je+1 ; do i=is-1,ie+1 |
| 872 | 0 | Mxx(i,j) = ((CS%Txx(i,j,k) * CS%c_diss(i,j,k)) * h(i,j,k)) * dy2h(i,j) |
| 873 | 0 | Myy(i,j) = ((CS%Tyy(i,j,k) * CS%c_diss(i,j,k)) * h(i,j,k)) * dx2h(i,j) |
| 874 | enddo ; enddo | |
| 875 | else | |
| 876 | 0 | do j=js-1,je+1 ; do i=is-1,ie+1 |
| 877 | 0 | Mxx(i,j) = ((CS%Txx(i,j,k)) * h(i,j,k)) * dy2h(i,j) |
| 878 | 0 | Myy(i,j) = ((CS%Tyy(i,j,k)) * h(i,j,k)) * dx2h(i,j) |
| 879 | enddo ; enddo | |
| 880 | endif | |
| 881 | ||
| 882 | ! Evaluate du/dt=1/h x.Div(h T) (Line 1495 of MOM_hor_visc.F90) | |
| 883 | 0 | do j=js,je ; do I=Isq,Ieq |
| 884 | 0 | h_u = 0.5 * (G%mask2dT(i,j)*h(i,j,k) + G%mask2dT(i+1,j)*h(i+1,j,k)) + h_neglect |
| 885 | fx = ((G%IdyCu(I,j)*(Mxx(i+1,j) - Mxx(i,j)) + & | |
| 886 | G%IdxCu(I,j)*((dx2q(I,J)*Mxy(I,J)) - (dx2q(I,J-1)*Mxy(I,J-1)))) * & | |
| 887 | 0 | G%IareaCu(I,j)) / h_u |
| 888 | 0 | diffu(I,j,k) = diffu(I,j,k) + fx |
| 889 | 0 | if (save_ZB2020u) & |
| 890 | 0 | ZB2020u(I,j,k) = fx |
| 891 | enddo ; enddo | |
| 892 | ||
| 893 | ! Evaluate dv/dt=1/h y.Div(h T) (Line 1517 of MOM_hor_visc.F90) | |
| 894 | 0 | do J=Jsq,Jeq ; do i=is,ie |
| 895 | 0 | h_v = 0.5 * (G%mask2dT(i,j)*h(i,j,k) + G%mask2dT(i,j+1)*h(i,j+1,k)) + h_neglect |
| 896 | fy = ((G%IdxCv(i,J)*(Myy(i,j+1) - Myy(i,j)) + & | |
| 897 | G%IdyCv(i,J)*((dy2q(I,J)*Mxy(I,J)) - (dy2q(I-1,J)*Mxy(I-1,J)))) * & | |
| 898 | 0 | G%IareaCv(i,J)) / h_v |
| 899 | 0 | diffv(i,J,k) = diffv(i,J,k) + fy |
| 900 | 0 | if (save_ZB2020v) & |
| 901 | 0 | ZB2020v(i,J,k) = fy |
| 902 | enddo ; enddo | |
| 903 | ||
| 904 | enddo ! end of k loop | |
| 905 | ||
| 906 | 0 | call cpu_clock_end(CS%id_clock_divergence) |
| 907 | ||
| 908 | 0 | call cpu_clock_begin(CS%id_clock_post) |
| 909 | 0 | if (CS%id_ZB2020u>0) call post_data(CS%id_ZB2020u, ZB2020u, CS%diag) |
| 910 | 0 | if (CS%id_ZB2020v>0) call post_data(CS%id_ZB2020v, ZB2020v, CS%diag) |
| 911 | 0 | call cpu_clock_end(CS%id_clock_post) |
| 912 | ||
| 913 | 0 | call compute_energy_source(u, v, h, ZB2020u, ZB2020v, G, GV, CS) |
| 914 | ||
| 915 | 0 | end subroutine compute_stress_divergence |
| 916 | ||
| 917 | !> Filtering of the velocity gradients sh_xx, sh_xy, vort_xy. | |
| 918 | !! Here instead of smoothing we do sharpening, i.e. | |
| 919 | !! return (initial - smoothed) fields. | |
| 920 | !! The algorithm: marching halo with non-blocking grouped MPI | |
| 921 | !! exchanges. The input array sh_xx should have halo 2 with | |
| 922 | !! applied zero B.C. The arrays sh_xy and vort_xy should have | |
| 923 | !! halo 1 with applied B.C. The output have the same halo and B.C. | |
| 924 | 0 | subroutine filter_velocity_gradients(G, GV, CS) |
| 925 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. | |
| 926 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure | |
| 927 | type(ZB2020_CS), intent(inout) :: CS !< ZB2020 control structure. | |
| 928 | ||
| 929 | real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: & | |
| 930 | 0 | sh_xx ! Copy of CS%sh_xx [T-1 ~> s-1] |
| 931 | real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)) :: & | |
| 932 | 0 | sh_xy, vort_xy ! Copy of CS%sh_xy and CS%vort_xy [T-1 ~> s-1] |
| 933 | ||
| 934 | integer :: xx_halo, xy_halo, vort_halo ! currently available halo for gradient components | |
| 935 | integer :: xx_iter, xy_iter, vort_iter ! remaining number of iterations | |
| 936 | integer :: niter ! required number of iterations | |
| 937 | ||
| 938 | integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz | |
| 939 | integer :: i, j, k | |
| 940 | ||
| 941 | 0 | niter = CS%HPF_iter |
| 942 | ||
| 943 | 0 | if (niter == 0) return |
| 944 | ||
| 945 | 0 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke |
| 946 | 0 | Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB |
| 947 | ||
| 948 | 0 | if (.not. G%symmetric) & |
| 949 | call do_group_pass(CS%pass_xx, G%Domain, & | |
| 950 | 0 | clock=CS%id_clock_mpi) |
| 951 | ||
| 952 | ! This is just copy of the array | |
| 953 | 0 | call cpu_clock_begin(CS%id_clock_filter) |
| 954 | 0 | do k=1,nz |
| 955 | ! Halo of size 2 is valid | |
| 956 | 0 | do j=js-2,je+2 ; do i=is-2,ie+2 |
| 957 | 0 | sh_xx(i,j,k) = CS%sh_xx(i,j,k) |
| 958 | enddo ; enddo | |
| 959 | ! Only halo of size 1 is valid | |
| 960 | 0 | do J=Jsq-1,Jeq+1 ; do I=Isq-1,Ieq+1 |
| 961 | 0 | sh_xy(I,J,k) = CS%sh_xy(I,J,k) |
| 962 | 0 | vort_xy(I,J,k) = CS%vort_xy(I,J,k) |
| 963 | enddo ; enddo | |
| 964 | enddo | |
| 965 | 0 | call cpu_clock_end(CS%id_clock_filter) |
| 966 | ||
| 967 | 0 | xx_halo = 2 ; xy_halo = 1 ; vort_halo = 1 |
| 968 | 0 | xx_iter = niter ; xy_iter = niter ; vort_iter = niter |
| 969 | ||
| 970 | do while & | |
| 971 | (xx_iter > 0 .or. xy_iter > 0 .or. & ! filter iterations remain to be done | |
| 972 | 0 | xx_halo < 2 .or. xy_halo < 1) ! there is no halo for VG tensor |
| 973 | ||
| 974 | ! ---------- filtering sh_xx --------- | |
| 975 | 0 | if (xx_halo < 2) then |
| 976 | 0 | call complete_group_pass(CS%pass_xx, G%Domain, clock=CS%id_clock_mpi) |
| 977 | 0 | xx_halo = CS%HPF_halo |
| 978 | endif | |
| 979 | ||
| 980 | 0 | call filter_hq(G, GV, CS, xx_halo, xx_iter, h=CS%sh_xx) |
| 981 | ||
| 982 | 0 | if (xx_halo < 2) & |
| 983 | 0 | call start_group_pass(CS%pass_xx, G%Domain, clock=CS%id_clock_mpi) |
| 984 | ||
| 985 | ! ------ filtering sh_xy, vort_xy ---- | |
| 986 | 0 | if (xy_halo < 1) then |
| 987 | 0 | call complete_group_pass(CS%pass_xy, G%Domain, clock=CS%id_clock_mpi) |
| 988 | 0 | xy_halo = CS%HPF_halo ; vort_halo = CS%HPF_halo |
| 989 | endif | |
| 990 | ||
| 991 | 0 | call filter_hq(G, GV, CS, xy_halo, xy_iter, q=CS%sh_xy) |
| 992 | 0 | call filter_hq(G, GV, CS, vort_halo, vort_iter, q=CS%vort_xy) |
| 993 | ||
| 994 | 0 | if (xy_halo < 1) & |
| 995 | 0 | call start_group_pass(CS%pass_xy, G%Domain, clock=CS%id_clock_mpi) |
| 996 | ||
| 997 | enddo | |
| 998 | ||
| 999 | ! We implement sharpening by computing residual | |
| 1000 | ! B.C. are already applied to all fields | |
| 1001 | 0 | call cpu_clock_begin(CS%id_clock_filter) |
| 1002 | 0 | do k=1,nz |
| 1003 | 0 | do j=js-2,je+2 ; do i=is-2,ie+2 |
| 1004 | 0 | CS%sh_xx(i,j,k) = sh_xx(i,j,k) - CS%sh_xx(i,j,k) |
| 1005 | enddo ; enddo | |
| 1006 | 0 | do J=Jsq-1,Jeq+1 ; do I=Isq-1,Ieq+1 |
| 1007 | 0 | CS%sh_xy(I,J,k) = sh_xy(I,J,k) - CS%sh_xy(I,J,k) |
| 1008 | 0 | CS%vort_xy(I,J,k) = vort_xy(I,J,k) - CS%vort_xy(I,J,k) |
| 1009 | enddo ; enddo | |
| 1010 | enddo | |
| 1011 | 0 | call cpu_clock_end(CS%id_clock_filter) |
| 1012 | ||
| 1013 | 0 | if (.not. G%symmetric) & |
| 1014 | call do_group_pass(CS%pass_xy, G%Domain, & | |
| 1015 | 0 | clock=CS%id_clock_mpi) |
| 1016 | ||
| 1017 | end subroutine filter_velocity_gradients | |
| 1018 | ||
| 1019 | !> Filtering of the stress tensor Txx, Tyy, Txy. | |
| 1020 | !! The algorithm: marching halo with non-blocking grouped MPI | |
| 1021 | !! exchanges. The input arrays (Txx, Tyy, Txy) must have halo 1 | |
| 1022 | !! with zero B.C. applied. The output have the same halo and B.C. | |
| 1023 | 0 | subroutine filter_stress(G, GV, CS) |
| 1024 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. | |
| 1025 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure | |
| 1026 | type(ZB2020_CS), intent(inout) :: CS !< ZB2020 control structure. | |
| 1027 | ||
| 1028 | integer :: Txx_halo, Tyy_halo, Txy_halo ! currently available halo for stress components | |
| 1029 | integer :: Txx_iter, Tyy_iter, Txy_iter ! remaining number of iterations | |
| 1030 | integer :: niter ! required number of iterations | |
| 1031 | ||
| 1032 | 0 | niter = CS%Stress_iter |
| 1033 | ||
| 1034 | 0 | if (niter == 0) return |
| 1035 | ||
| 1036 | 0 | Txx_halo = 1 ; Tyy_halo = 1 ; Txy_halo = 1 ; ! these are required halo for Txx, Tyy, Txy |
| 1037 | 0 | Txx_iter = niter ; Tyy_iter = niter ; Txy_iter = niter |
| 1038 | ||
| 1039 | do while & | |
| 1040 | (Txx_iter > 0 .or. Txy_iter > 0 .or. & ! filter iterations remain to be done | |
| 1041 | 0 | Txx_halo < 1 .or. Txy_halo < 1) ! there is no halo for Txx or Txy |
| 1042 | ||
| 1043 | ! ---------- filtering Txy ----------- | |
| 1044 | 0 | if (Txy_halo < 1) then |
| 1045 | 0 | call complete_group_pass(CS%pass_Tq, G%Domain, clock=CS%id_clock_mpi) |
| 1046 | 0 | Txy_halo = CS%Stress_halo |
| 1047 | endif | |
| 1048 | ||
| 1049 | 0 | call filter_hq(G, GV, CS, Txy_halo, Txy_iter, q=CS%Txy) |
| 1050 | ||
| 1051 | 0 | if (Txy_halo < 1) & |
| 1052 | 0 | call start_group_pass(CS%pass_Tq, G%Domain, clock=CS%id_clock_mpi) |
| 1053 | ||
| 1054 | ! ------- filtering Txx, Tyy --------- | |
| 1055 | 0 | if (Txx_halo < 1) then |
| 1056 | 0 | call complete_group_pass(CS%pass_Th, G%Domain, clock=CS%id_clock_mpi) |
| 1057 | 0 | Txx_halo = CS%Stress_halo ; Tyy_halo = CS%Stress_halo |
| 1058 | endif | |
| 1059 | ||
| 1060 | 0 | call filter_hq(G, GV, CS, Txx_halo, Txx_iter, h=CS%Txx) |
| 1061 | 0 | call filter_hq(G, GV, CS, Tyy_halo, Tyy_iter, h=CS%Tyy) |
| 1062 | ||
| 1063 | 0 | if (Txx_halo < 1) & |
| 1064 | 0 | call start_group_pass(CS%pass_Th, G%Domain, clock=CS%id_clock_mpi) |
| 1065 | ||
| 1066 | enddo | |
| 1067 | ||
| 1068 | end subroutine filter_stress | |
| 1069 | ||
| 1070 | !> Wrapper for filter_3D function. The border indices for q and h | |
| 1071 | !! arrays are substituted. | |
| 1072 | 0 | subroutine filter_hq(G, GV, CS, current_halo, remaining_iterations, q, h) |
| 1073 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. | |
| 1074 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure | |
| 1075 | type(ZB2020_CS), intent(in) :: CS !< ZB2020 control structure. | |
| 1076 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), optional, & | |
| 1077 | intent(inout) :: h !< Input/output array in h points [arbitrary] | |
| 1078 | real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)), optional, & | |
| 1079 | intent(inout) :: q !< Input/output array in q points [arbitrary] | |
| 1080 | integer, intent(inout) :: current_halo !< Currently available halo points | |
| 1081 | integer, intent(inout) :: remaining_iterations !< The number of iterations to perform | |
| 1082 | ||
| 1083 | logical :: direction ! The direction of the first 1D filter | |
| 1084 | ||
| 1085 | 0 | direction = (MOD(G%first_direction,2) == 0) |
| 1086 | ||
| 1087 | 0 | call cpu_clock_begin(CS%id_clock_filter) |
| 1088 | ||
| 1089 | 0 | if (present(h)) then |
| 1090 | call filter_3D(h, CS%maskw_h, & | |
| 1091 | G%isd, G%ied, G%jsd, G%jed, & | |
| 1092 | G%isc, G%iec, G%jsc, G%jec, GV%ke, & | |
| 1093 | current_halo, remaining_iterations, & | |
| 1094 | 0 | direction) |
| 1095 | endif | |
| 1096 | ||
| 1097 | 0 | if (present(q)) then |
| 1098 | call filter_3D(q, CS%maskw_q, & | |
| 1099 | G%IsdB, G%IedB, G%JsdB, G%JedB, & | |
| 1100 | G%IscB, G%IecB, G%JscB, G%JecB, GV%ke, & | |
| 1101 | current_halo, remaining_iterations, & | |
| 1102 | 0 | direction) |
| 1103 | endif | |
| 1104 | ||
| 1105 | 0 | call cpu_clock_end(CS%id_clock_filter) |
| 1106 | 0 | end subroutine filter_hq |
| 1107 | ||
| 1108 | !> Spatial lateral filter applied to 3D array. The lateral filter is given | |
| 1109 | !! by the convolutional kernel: | |
| 1110 | !! [1 2 1] | |
| 1111 | !! C = |2 4 2| * 1/16 | |
| 1112 | !! [1 2 1] | |
| 1113 | !! The fast algorithm decomposes the 2D filter into two 1D filters as follows: | |
| 1114 | !! [1] | |
| 1115 | !! C = |2| * [1 2 1] * 1/16 | |
| 1116 | !! [1] | |
| 1117 | !! The input array must have zero B.C. applied. B.C. is applied for output array. | |
| 1118 | !! Note that maskw contains both land mask and 1/16 factor. | |
| 1119 | !! Filter implements marching halo. The available halo is specified and as many | |
| 1120 | !! filter iterations as possible and as needed are performed. | |
| 1121 | 0 | subroutine filter_3D(x, maskw, isd, ied, jsd, jed, is, ie, js, je, nz, & |
| 1122 | current_halo, remaining_iterations, & | |
| 1123 | direction) | |
| 1124 | integer, intent(in) :: isd !< Indices of array size | |
| 1125 | integer, intent(in) :: ied !< Indices of array size | |
| 1126 | integer, intent(in) :: jsd !< Indices of array size | |
| 1127 | integer, intent(in) :: jed !< Indices of array size | |
| 1128 | integer, intent(in) :: is !< Indices of owned points | |
| 1129 | integer, intent(in) :: ie !< Indices of owned points | |
| 1130 | integer, intent(in) :: js !< Indices of owned points | |
| 1131 | integer, intent(in) :: je !< Indices of owned points | |
| 1132 | integer, intent(in) :: nz !< Vertical array size | |
| 1133 | real, dimension(isd:ied,jsd:jed,nz), & | |
| 1134 | intent(inout) :: x !< Input/output array [arbitrary] | |
| 1135 | real, dimension(isd:ied,jsd:jed), & | |
| 1136 | intent(in) :: maskw !< Mask array of land points divided by 16 [nondim] | |
| 1137 | integer, intent(inout) :: current_halo !< Currently available halo points | |
| 1138 | integer, intent(inout) :: remaining_iterations !< The number of iterations to perform | |
| 1139 | logical, intent(in) :: direction !< The direction of the first 1D filter | |
| 1140 | ||
| 1141 | real, parameter :: weight = 2. ! Filter weight [nondim] | |
| 1142 | integer :: i, j, k, iter, niter, halo | |
| 1143 | ||
| 1144 | 0 | real :: tmp(isd:ied, jsd:jed) ! Array with temporary results [arbitrary] |
| 1145 | ||
| 1146 | ! Do as many iterations as needed and possible | |
| 1147 | 0 | niter = min(current_halo, remaining_iterations) |
| 1148 | 0 | if (niter == 0) return ! nothing to do |
| 1149 | ||
| 1150 | ! Update remaining iterations | |
| 1151 | 0 | remaining_iterations = remaining_iterations - niter |
| 1152 | ! Update halo information | |
| 1153 | 0 | current_halo = current_halo - niter |
| 1154 | ||
| 1155 | 0 | do k=1,Nz |
| 1156 | halo = niter-1 + & | |
| 1157 | 0 | current_halo ! Save as many halo points as possible |
| 1158 | 0 | do iter=1,niter |
| 1159 | ||
| 1160 | 0 | if (direction) then |
| 1161 | 0 | do j = js-halo, je+halo ; do i = is-halo-1, ie+halo+1 |
| 1162 | 0 | tmp(i,j) = weight * x(i,j,k) + (x(i,j-1,k) + x(i,j+1,k)) |
| 1163 | enddo ; enddo | |
| 1164 | ||
| 1165 | 0 | do j = js-halo, je+halo ; do i = is-halo, ie+halo |
| 1166 | 0 | x(i,j,k) = (weight * tmp(i,j) + (tmp(i-1,j) + tmp(i+1,j))) * maskw(i,j) |
| 1167 | enddo ; enddo | |
| 1168 | else | |
| 1169 | 0 | do j = js-halo-1, je+halo+1 ; do i = is-halo, ie+halo |
| 1170 | 0 | tmp(i,j) = weight * x(i,j,k) + (x(i-1,j,k) + x(i+1,j,k)) |
| 1171 | enddo ; enddo | |
| 1172 | ||
| 1173 | 0 | do j = js-halo, je+halo ; do i = is-halo, ie+halo |
| 1174 | 0 | x(i,j,k) = (weight * tmp(i,j) + (tmp(i,j-1) + tmp(i,j+1))) * maskw(i,j) |
| 1175 | enddo ; enddo | |
| 1176 | endif | |
| 1177 | ||
| 1178 | 0 | halo = halo - 1 |
| 1179 | enddo | |
| 1180 | enddo | |
| 1181 | ||
| 1182 | end subroutine filter_3D | |
| 1183 | ||
| 1184 | !> Computes the 3D energy source term for the ZB2020 scheme | |
| 1185 | !! similarly to MOM_diagnostics.F90, specifically 1125 line. | |
| 1186 | 0 | subroutine compute_energy_source(u, v, h, fx, fy, G, GV, CS) |
| 1187 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. | |
| 1188 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. | |
| 1189 | type(ZB2020_CS), intent(in) :: CS !< ZB2020 control structure. | |
| 1190 | ||
| 1191 | real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & | |
| 1192 | intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1]. | |
| 1193 | real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & | |
| 1194 | intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1]. | |
| 1195 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & | |
| 1196 | intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]. | |
| 1197 | ||
| 1198 | real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & | |
| 1199 | intent(in) :: fx !< Zonal acceleration due to convergence of | |
| 1200 | !! along-coordinate stress tensor [L T-2 ~> m s-2] | |
| 1201 | real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & | |
| 1202 | intent(in) :: fy !< Meridional acceleration due to convergence | |
| 1203 | !! of along-coordinate stress tensor [L T-2 ~> m s-2] | |
| 1204 | ||
| 1205 | 0 | real :: KE_term(SZI_(G),SZJ_(G),SZK_(GV)) ! A term in the kinetic energy budget |
| 1206 | ! [H L2 T-3 ~> m3 s-3 or W m-2] | |
| 1207 | 0 | real :: KE_u(SZIB_(G),SZJ_(G)) ! The area integral of a KE term in a layer at u-points |
| 1208 | ! [H L4 T-3 ~> m5 s-3 or kg m2 s-3] | |
| 1209 | 0 | real :: KE_v(SZI_(G),SZJB_(G)) ! The area integral of a KE term in a layer at v-points |
| 1210 | ! [H L4 T-3 ~> m5 s-3 or kg m2 s-3] | |
| 1211 | ||
| 1212 | real :: uh ! Transport through zonal faces = u*h*dy, | |
| 1213 | ! [H L2 T-1 ~> m3 s-1 or kg s-1]. | |
| 1214 | real :: vh ! Transport through meridional faces = v*h*dx, | |
| 1215 | ! [H L2 T-1 ~> m3 s-1 or kg s-1]. | |
| 1216 | ||
| 1217 | type(group_pass_type) :: pass_KE_uv ! A handle used for group halo passes | |
| 1218 | ||
| 1219 | integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz | |
| 1220 | integer :: i, j, k | |
| 1221 | ||
| 1222 | 0 | if (CS%id_KE_ZB2020 > 0) then |
| 1223 | 0 | call cpu_clock_begin(CS%id_clock_source) |
| 1224 | 0 | call create_group_pass(pass_KE_uv, KE_u, KE_v, G%Domain, To_North+To_East) |
| 1225 | ||
| 1226 | 0 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke |
| 1227 | 0 | Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB |
| 1228 | ||
| 1229 | 0 | KE_term(:,:,:) = 0. |
| 1230 | ! Calculate the KE source from Zanna-Bolton2020 [H L2 T-3 ~> m3 s-3]. | |
| 1231 | 0 | do k=1,nz |
| 1232 | 0 | KE_u(:,:) = 0. |
| 1233 | 0 | KE_v(:,:) = 0. |
| 1234 | 0 | do j=js,je ; do I=Isq,Ieq |
| 1235 | uh = u(I,j,k) * 0.5 * (G%mask2dT(i,j)*h(i,j,k) + G%mask2dT(i+1,j)*h(i+1,j,k)) * & | |
| 1236 | 0 | G%dyCu(I,j) |
| 1237 | 0 | KE_u(I,j) = uh * G%dxCu(I,j) * fx(I,j,k) |
| 1238 | enddo ; enddo | |
| 1239 | 0 | do J=Jsq,Jeq ; do i=is,ie |
| 1240 | vh = v(i,J,k) * 0.5 * (G%mask2dT(i,j)*h(i,j,k) + G%mask2dT(i,j+1)*h(i,j+1,k)) * & | |
| 1241 | 0 | G%dxCv(i,J) |
| 1242 | 0 | KE_v(i,J) = vh * G%dyCv(i,J) * fy(i,J,k) |
| 1243 | enddo ; enddo | |
| 1244 | 0 | call do_group_pass(pass_KE_uv, G%domain, clock=CS%id_clock_mpi) |
| 1245 | 0 | do j=js,je ; do i=is,ie |
| 1246 | KE_term(i,j,k) = 0.5 * G%IareaT(i,j) & | |
| 1247 | 0 | * ((KE_u(I,j) + KE_u(I-1,j)) + (KE_v(i,J) + KE_v(i,J-1))) |
| 1248 | enddo ; enddo | |
| 1249 | enddo | |
| 1250 | ||
| 1251 | 0 | call cpu_clock_end(CS%id_clock_source) |
| 1252 | ||
| 1253 | 0 | call cpu_clock_begin(CS%id_clock_post) |
| 1254 | 0 | call post_data(CS%id_KE_ZB2020, KE_term, CS%diag) |
| 1255 | 0 | call cpu_clock_end(CS%id_clock_post) |
| 1256 | endif | |
| 1257 | ||
| 1258 | 0 | end subroutine compute_energy_source |
| 1259 | ||
| 1260 | 0 | end module MOM_Zanna_Bolton |