MOM6
Data Types | Functions/Subroutines
mom_hor_visc Module Reference

Data Types

type  hor_visc_cs
 

Functions/Subroutines

subroutine, public horizontal_viscosity (u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, OBC)
 This subroutine determines the acceleration due to the horizontal viscosity. A combination of biharmonic and Laplacian forms can be used. The coefficient may either be a constant or a shear-dependent form. The biharmonic is determined by twice taking the divergence of an appropriately defined stress tensor. The Laplacian is determined by doing so once. To work, the following fields must be set outside of the usual is to ie range before this subroutine is called: v[is-2,is-1,ie+1,ie+2], u[is-2,is-1,ie+1,ie+2], and h[is-1,ie+1], with a similarly sized halo in the y-direction. More...
 
subroutine, public hor_visc_init (Time, G, param_file, diag, CS)
 This subroutine allocates space for and calculates static variables used by this module. The metrics may be 0, 1, or 2-D arrays, while fields like the background viscosities are 2-D arrays. ALLOC is a macro defined in MOM_memory.h to either allocate for dynamic memory, or do nothing when using static memory. More...
 
subroutine, public hor_visc_end (CS)
 

Function/Subroutine Documentation

◆ hor_visc_end()

subroutine, public mom_hor_visc::hor_visc_end ( type(hor_visc_cs), pointer  CS)

Definition at line 1575 of file MOM_hor_visc.F90.

1575 ! This subroutine deallocates any variables allocated in hor_visc_init.
1576 ! Argument: CS - The control structure returned by a previous call to
1577 ! hor_visc_init.
1578  type(hor_visc_cs), pointer :: cs
1579 
1580  if (cs%Laplacian .or. cs%biharmonic) then
1581  dealloc_(cs%dx2h) ; dealloc_(cs%dx2q) ; dealloc_(cs%dy2h) ; dealloc_(cs%dy2q)
1582  dealloc_(cs%dx_dyT) ; dealloc_(cs%dy_dxT) ; dealloc_(cs%dx_dyBu) ; dealloc_(cs%dy_dxBu)
1583  dealloc_(cs%reduction_xx) ; dealloc_(cs%reduction_xy)
1584  endif
1585 
1586  if (cs%Laplacian) then
1587  dealloc_(cs%Kh_bg_xx) ; dealloc_(cs%Kh_bg_xy)
1588  if (cs%bound_Kh) then
1589  dealloc_(cs%Kh_Max_xx) ; dealloc_(cs%Kh_Max_xy)
1590  endif
1591  if (cs%Smagorinsky_Kh) then
1592  dealloc_(cs%Laplac_Const_xx) ; dealloc_(cs%Laplac_Const_xy)
1593  endif
1594  if (cs%Leith_Kh) then
1595  dealloc_(cs%Laplac3_Const_xx) ; dealloc_(cs%Laplac3_Const_xy)
1596  endif
1597  endif
1598 
1599  if (cs%biharmonic) then
1600  dealloc_(cs%Idx2dyCu) ; dealloc_(cs%Idx2dyCv)
1601  dealloc_(cs%Idxdy2u) ; dealloc_(cs%Idxdy2v)
1602  dealloc_(cs%Ah_bg_xx) ; dealloc_(cs%Ah_bg_xy)
1603  if (cs%bound_Ah) then
1604  dealloc_(cs%Ah_Max_xx) ; dealloc_(cs%Ah_Max_xy)
1605  endif
1606  if (cs%Smagorinsky_Ah) then
1607  dealloc_(cs%Biharm_Const_xx) ; dealloc_(cs%Biharm_Const_xy)
1608  if (cs%bound_Coriolis) then
1609  dealloc_(cs%Biharm_Const2_xx) ; dealloc_(cs%Biharm_Const2_xy)
1610  endif
1611  endif
1612  if (cs%Leith_Ah) then
1613  dealloc_(cs%Biharm5_Const_xx) ; dealloc_(cs%Biharm5_Const_xy)
1614  endif
1615  endif
1616  deallocate(cs)
1617 

◆ hor_visc_init()

subroutine, public mom_hor_visc::hor_visc_init ( type(time_type), intent(in)  Time,
type(ocean_grid_type), intent(inout)  G,
type(param_file_type), intent(in)  param_file,
type(diag_ctrl), intent(inout), target  diag,
type(hor_visc_cs), pointer  CS 
)

This subroutine allocates space for and calculates static variables used by this module. The metrics may be 0, 1, or 2-D arrays, while fields like the background viscosities are 2-D arrays. ALLOC is a macro defined in MOM_memory.h to either allocate for dynamic memory, or do nothing when using static memory.

Parameters
[in]timecurrent model time.
[in,out]gThe ocean's grid structure.
[in]param_fileA structure to parse for run-time parameters.
[in,out]diagstructure to regulate diagnostic output.
cspointer to the control structure for this module

Definition at line 995 of file MOM_hor_visc.F90.

References mom_error_handler::mom_error().

995  type(time_type), intent(in) :: time !< current model time.
996  type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
997  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
998  !! parameters.
999  type(diag_ctrl), target, intent(inout) :: diag !< structure to regulate diagnostic output.
1000  type(hor_visc_cs), pointer :: cs !< pointer to the control structure for this module
1001 
1002 ! This subroutine allocates space for and calculates static variables
1003 ! used by this module. The metrics may be 0, 1, or 2-D arrays,
1004 ! while fields like the background viscosities are 2-D arrays.
1005 ! ALLOC is a macro defined in MOM_memory.h to either allocate
1006 ! for dynamic memory, or do nothing when using static memory.
1007 !
1008 ! Arguments:
1009 ! (in) Time - current model time
1010 ! (in) G - ocean grid structure
1011 ! (in) param_file - structure to parse for model parameter values
1012 ! (in) diag - structure to regulate diagnostic output
1013 ! (in/out) CS - pointer to the control structure for this module
1014 
1015  real, dimension(SZIB_(G),SZJ_(G)) :: u0u, u0v
1016  real, dimension(SZI_(G),SZJB_(G)) :: v0u, v0v
1017  ! u0v is the Laplacian sensitivities to the v velocities
1018  ! at u points, in m-2, with u0u, v0u, and v0v defined similarly.
1019  real :: grid_sp_h2 ! Harmonic mean of the squares of the grid
1020  real :: grid_sp_h3 ! Harmonic mean of the squares of the grid^(3/2)
1021  real :: grid_sp_q2 ! spacings at h and q points (m2)
1022  real :: grid_sp_q3 ! spacings at h and q points^(3/2) (m3)
1023  real :: kh_limit ! A coefficient (1/s) used, along with the
1024  ! grid spacing, to limit Laplacian viscosity.
1025  real :: fmax ! maximum absolute value of f at the four
1026  ! vorticity points around a thickness point (1/s)
1027  real :: boundcorconst ! constant (s2/m2)
1028  real :: ah_limit ! coefficient (1/s) used, along with the
1029  ! grid spacing, to limit biharmonic viscosity
1030  real :: kh ! Lapacian horizontal viscosity (m2/s)
1031  real :: ah ! biharmonic horizontal viscosity (m4/s)
1032  real :: kh_vel_scale ! this speed (m/s) times grid spacing gives Lap visc
1033  real :: ah_vel_scale ! this speed (m/s) times grid spacing cubed gives bih visc
1034  real :: smag_lap_const ! nondimensional Laplacian Smagorinsky constant
1035  real :: smag_bi_const ! nondimensional biharmonic Smagorinsky constant
1036  real :: leith_lap_const ! nondimensional Laplacian Leith constant
1037  real :: leith_bi_const ! nondimensional biharmonic Leith constant
1038  real :: dt ! dynamics time step (sec)
1039  real :: idt ! inverse of dt (1/s)
1040  real :: denom ! work variable; the denominator of a fraction
1041  real :: maxvel ! largest permitted velocity components (m/s)
1042  real :: bound_cor_vel ! grid-scale velocity variations at which value
1043  ! the quadratically varying biharmonic viscosity
1044  ! balances Coriolis acceleration (m/s)
1045  logical :: bound_cor_def ! parameter setting of BOUND_CORIOLIS
1046  logical :: get_all ! If true, read and log all parameters, regardless of
1047  ! whether they are used, to enable spell-checking of
1048  ! valid parameters.
1049  character(len=64) :: inputdir, filename
1050  integer :: is, ie, js, je, isq, ieq, jsq, jeq, nz
1051  integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
1052  integer :: i, j
1053 
1054 ! This include declares and sets the variable "version".
1055 #include "version_variable.h"
1056  character(len=40) :: mdl = "MOM_hor_visc" ! module name
1057 
1058  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1059  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1060  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1061  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
1062 
1063  if (associated(cs)) then
1064  call mom_error(warning, "hor_visc_init called with an associated "// &
1065  "control structure.")
1066  return
1067  endif
1068  allocate(cs)
1069 
1070  cs%diag => diag
1071 
1072  ! Read parameters and write them to the model log.
1073  call log_version(param_file, mdl, version, "")
1074 
1075  ! It is not clear whether these initialization lines are needed for the
1076  ! cases where the corresponding parameters are not read.
1077  cs%bound_Kh = .false. ; cs%better_bound_Kh = .false. ; cs%Smagorinsky_Kh = .false. ; cs%Leith_Kh = .false.
1078  cs%bound_Ah = .false. ; cs%better_bound_Ah = .false. ; cs%Smagorinsky_Ah = .false. ; cs%Leith_Ah = .false.
1079  cs%bound_Coriolis = .false.
1080  cs%Modified_Leith = .false.
1081 
1082  kh = 0.0 ; ah = 0.0
1083 
1084  ! If GET_ALL_PARAMS is true, all parameters are read in all cases to enable
1085  ! parameter spelling checks.
1086  call get_param(param_file, mdl, "GET_ALL_PARAMS", get_all, default=.false.)
1087 
1088  call get_param(param_file, mdl, "LAPLACIAN", cs%Laplacian, &
1089  "If true, use a Laplacian horizontal viscosity.", &
1090  default=.false.)
1091  if (cs%Laplacian .or. get_all) then
1092  call get_param(param_file, mdl, "KH", kh, &
1093  "The background Laplacian horizontal viscosity.", &
1094  units = "m2 s-1", default=0.0)
1095  call get_param(param_file, mdl, "KH_BG_MIN", cs%Kh_bg_min, &
1096  "The minimum value allowed for Laplacian horizontal viscosity, KH.", &
1097  units = "m2 s-1", default=0.0)
1098  call get_param(param_file, mdl, "KH_VEL_SCALE", kh_vel_scale, &
1099  "The velocity scale which is multiplied by the grid \n"//&
1100  "spacing to calculate the Laplacian viscosity. \n"//&
1101  "The final viscosity is the largest of this scaled \n"//&
1102  "viscosity, the Smagorinsky and Leith viscosities, and KH.", &
1103  units="m s-1", default=0.0)
1104 
1105  call get_param(param_file, mdl, "SMAGORINSKY_KH", cs%Smagorinsky_Kh, &
1106  "If true, use a Smagorinsky nonlinear eddy viscosity.", &
1107  default=.false.)
1108  if (cs%Smagorinsky_Kh .or. get_all) &
1109  call get_param(param_file, mdl, "SMAG_LAP_CONST", smag_lap_const, &
1110  "The nondimensional Laplacian Smagorinsky constant, \n"//&
1111  "often 0.15.", units="nondim", default=0.0, &
1112  fail_if_missing = cs%Smagorinsky_Kh)
1113 
1114  call get_param(param_file, mdl, "LEITH_KH", cs%Leith_Kh, &
1115  "If true, use a Leith nonlinear eddy viscosity.", &
1116  default=.false.)
1117 
1118  call get_param(param_file, mdl, "MODIFIED_LEITH", cs%Modified_Leith, &
1119  "If true, add a term to Leith viscosity which is \n"//&
1120  "proportional to the gradient of divergence.", &
1121  default=.false.)
1122 
1123  if (cs%Leith_Kh .or. get_all) &
1124  call get_param(param_file, mdl, "LEITH_LAP_CONST", leith_lap_const, &
1125  "The nondimensional Laplacian Leith constant, \n"//&
1126  "often ??", units="nondim", default=0.0, &
1127  fail_if_missing = cs%Leith_Kh)
1128 
1129  call get_param(param_file, mdl, "BOUND_KH", cs%bound_Kh, &
1130  "If true, the Laplacian coefficient is locally limited \n"//&
1131  "to be stable.", default=.true.)
1132  call get_param(param_file, mdl, "BETTER_BOUND_KH", cs%better_bound_Kh, &
1133  "If true, the Laplacian coefficient is locally limited \n"//&
1134  "to be stable with a better bounding than just BOUND_KH.", &
1135  default=cs%bound_Kh)
1136  endif
1137 
1138  call get_param(param_file, mdl, "BIHARMONIC", cs%biharmonic, &
1139  "If true, use a biharmonic horizontal viscosity. \n"//&
1140  "BIHARMONIC may be used with LAPLACIAN.", &
1141  default=.true.)
1142  if (cs%biharmonic .or. get_all) then
1143  call get_param(param_file, mdl, "AH", ah, &
1144  "The background biharmonic horizontal viscosity.", &
1145  units = "m4 s-1", default=0.0)
1146  call get_param(param_file, mdl, "AH_VEL_SCALE", ah_vel_scale, &
1147  "The velocity scale which is multiplied by the cube of \n"//&
1148  "the grid spacing to calculate the biharmonic viscosity. \n"//&
1149  "The final viscosity is the largest of this scaled \n"//&
1150  "viscosity, the Smagorinsky and Leith viscosities, and AH.", &
1151  units="m s-1", default=0.0)
1152  call get_param(param_file, mdl, "SMAGORINSKY_AH", cs%Smagorinsky_Ah, &
1153  "If true, use a biharmonic Smagorinsky nonlinear eddy \n"//&
1154  "viscosity.", default=.false.)
1155  call get_param(param_file, mdl, "LEITH_AH", cs%Leith_Ah, &
1156  "If true, use a biharmonic Leith nonlinear eddy \n"//&
1157  "viscosity.", default=.false.)
1158 
1159  call get_param(param_file, mdl, "BOUND_AH", cs%bound_Ah, &
1160  "If true, the biharmonic coefficient is locally limited \n"//&
1161  "to be stable.", default=.true.)
1162  call get_param(param_file, mdl, "BETTER_BOUND_AH", cs%better_bound_Ah, &
1163  "If true, the biharmonic coefficient is locally limited \n"//&
1164  "to be stable with a better bounding than just BOUND_AH.", &
1165  default=cs%bound_Ah)
1166 
1167  if (cs%Smagorinsky_Ah .or. get_all) then
1168  call get_param(param_file, mdl, "SMAG_BI_CONST",smag_bi_const, &
1169  "The nondimensional biharmonic Smagorinsky constant, \n"//&
1170  "typically 0.015 - 0.06.", units="nondim", default=0.0, &
1171  fail_if_missing = cs%Smagorinsky_Ah)
1172 
1173  call get_param(param_file, mdl, "BOUND_CORIOLIS", bound_cor_def, default=.false.)
1174  call get_param(param_file, mdl, "BOUND_CORIOLIS_BIHARM", cs%bound_Coriolis, &
1175  "If true use a viscosity that increases with the square \n"//&
1176  "of the velocity shears, so that the resulting viscous \n"//&
1177  "drag is of comparable magnitude to the Coriolis terms \n"//&
1178  "when the velocity differences between adjacent grid \n"//&
1179  "points is 0.5*BOUND_CORIOLIS_VEL. The default is the \n"//&
1180  "value of BOUND_CORIOLIS (or false).", default=bound_cor_def)
1181  if (cs%bound_Coriolis .or. get_all) then
1182  call get_param(param_file, mdl, "MAXVEL", maxvel, default=3.0e8)
1183  bound_cor_vel = maxvel
1184  call get_param(param_file, mdl, "BOUND_CORIOLIS_VEL", bound_cor_vel, &
1185  "The velocity scale at which BOUND_CORIOLIS_BIHARM causes \n"//&
1186  "the biharmonic drag to have comparable magnitude to the \n"//&
1187  "Coriolis acceleration. The default is set by MAXVEL.", &
1188  units="m s-1", default=maxvel)
1189  endif
1190  endif
1191 
1192  if (cs%Leith_Ah .or. get_all) then
1193  call get_param(param_file, mdl, "LEITH_BI_CONST",leith_bi_const, &
1194  "The nondimensional biharmonic Leith constant, \n"//&
1195  "typical values are thus far undetermined", units="nondim", default=0.0, &
1196  fail_if_missing = cs%Leith_Ah)
1197  endif
1198 
1199  endif
1200 
1201  if (cs%better_bound_Ah .or. cs%better_bound_Kh .or. get_all) &
1202  call get_param(param_file, mdl, "HORVISC_BOUND_COEF", cs%bound_coef, &
1203  "The nondimensional coefficient of the ratio of the \n"//&
1204  "viscosity bounds to the theoretical maximum for \n"//&
1205  "stability without considering other terms.", units="nondim", &
1206  default=0.8)
1207 
1208  call get_param(param_file, mdl, "NOSLIP", cs%no_slip, &
1209  "If true, no slip boundary conditions are used; otherwise \n"//&
1210  "free slip boundary conditions are assumed. The \n"//&
1211  "implementation of the free slip BCs on a C-grid is much \n"//&
1212  "cleaner than the no slip BCs. The use of free slip BCs \n"//&
1213  "is strongly encouraged, and no slip BCs are not used with \n"//&
1214  "the biharmonic viscosity.", default=.false.)
1215 
1216  call get_param(param_file, mdl, "USE_KH_BG_2D", cs%use_Kh_bg_2d, &
1217  "If true, read a file containing 2-d background harmonic \n"//&
1218  "viscosities. The final viscosity is the maximum of the other "//&
1219  "terms and this background value.", default=.false.)
1220 
1221 
1222  if (cs%bound_Kh .or. cs%bound_Ah .or. cs%better_bound_Kh .or. cs%better_bound_Ah) &
1223  call get_param(param_file, mdl, "DT", dt, &
1224  "The (baroclinic) dynamics time step.", units = "s", &
1225  fail_if_missing=.true.)
1226 
1227  if (cs%no_slip .and. cs%biharmonic) &
1228  call mom_error(fatal,"ERROR: NOSLIP and BIHARMONIC cannot be defined "// &
1229  "at the same time in MOM.")
1230 
1231  if (.not.(cs%Laplacian .or. cs%biharmonic)) then
1232  ! Only issue inviscid warning if not in single column mode (usually 2x2 domain)
1233  if ( max(g%domain%niglobal, g%domain%njglobal)>2 ) call mom_error(warning, &
1234  "hor_visc_init: It is usually a very bad idea not to use either "//&
1235  "LAPLACIAN or BIHARMONIC viscosity.")
1236  return ! We are not using either Laplacian or Bi-harmonic lateral viscosity
1237  endif
1238 
1239  alloc_(cs%dx2h(isd:ied,jsd:jed)) ; cs%dx2h(:,:) = 0.0
1240  alloc_(cs%dy2h(isd:ied,jsd:jed)) ; cs%dy2h(:,:) = 0.0
1241  alloc_(cs%dx2q(isdb:iedb,jsdb:jedb)) ; cs%dx2q(:,:) = 0.0
1242  alloc_(cs%dy2q(isdb:iedb,jsdb:jedb)) ; cs%dy2q(:,:) = 0.0
1243  alloc_(cs%dx_dyT(isd:ied,jsd:jed)) ; cs%dx_dyT(:,:) = 0.0
1244  alloc_(cs%dy_dxT(isd:ied,jsd:jed)) ; cs%dy_dxT(:,:) = 0.0
1245  alloc_(cs%dx_dyBu(isdb:iedb,jsdb:jedb)) ; cs%dx_dyBu(:,:) = 0.0
1246  alloc_(cs%dy_dxBu(isdb:iedb,jsdb:jedb)) ; cs%dy_dxBu(:,:) = 0.0
1247 
1248  if (cs%Laplacian) then
1249  alloc_(cs%Kh_bg_xx(isd:ied,jsd:jed)) ; cs%Kh_bg_xx(:,:) = 0.0
1250  alloc_(cs%Kh_bg_xy(isdb:iedb,jsdb:jedb)) ; cs%Kh_bg_xy(:,:) = 0.0
1251  if (cs%bound_Kh .or. cs%better_bound_Kh) then
1252  alloc_(cs%Kh_Max_xx(isdb:iedb,jsdb:jedb)) ; cs%Kh_Max_xx(:,:) = 0.0
1253  alloc_(cs%Kh_Max_xy(isdb:iedb,jsdb:jedb)) ; cs%Kh_Max_xy(:,:) = 0.0
1254  endif
1255  if (cs%Smagorinsky_Kh) then
1256  alloc_(cs%Laplac_Const_xx(isd:ied,jsd:jed)) ; cs%Laplac_Const_xx(:,:) = 0.0
1257  alloc_(cs%Laplac_Const_xy(isdb:iedb,jsdb:jedb)) ; cs%Laplac_Const_xy(:,:) = 0.0
1258  endif
1259  if (cs%Leith_Kh) then
1260  alloc_(cs%Laplac3_Const_xx(isd:ied,jsd:jed)) ; cs%Laplac3_Const_xx(:,:) = 0.0
1261  alloc_(cs%Laplac3_Const_xy(isdb:iedb,jsdb:jedb)) ; cs%Laplac3_Const_xy(:,:) = 0.0
1262  endif
1263 
1264  endif
1265  alloc_(cs%reduction_xx(isd:ied,jsd:jed)) ; cs%reduction_xx(:,:) = 0.0
1266  alloc_(cs%reduction_xy(isdb:iedb,jsdb:jedb)) ; cs%reduction_xy(:,:) = 0.0
1267 
1268  if (cs%use_Kh_bg_2d) then
1269  alloc_(cs%Kh_bg_2d(isd:ied,jsd:jed)) ; cs%Kh_bg_2d(:,:) = 0.0
1270  call get_param(param_file, mdl, "KH_BG_2D_FILENAME", filename, &
1271  'The filename containing a 2d map of "Kh".', &
1272  default='KH_background_2d.nc')
1273  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
1274  inputdir = slasher(inputdir)
1275  call read_data(trim(inputdir)//trim(filename), 'Kh', cs%Kh_bg_2d, &
1276  domain=g%domain%mpp_domain, timelevel=1)
1277  call pass_var(cs%Kh_bg_2d, g%domain)
1278  endif
1279 
1280  if (cs%biharmonic) then
1281  alloc_(cs%Idx2dyCu(isdb:iedb,jsd:jed)) ; cs%Idx2dyCu(:,:) = 0.0
1282  alloc_(cs%Idx2dyCv(isd:ied,jsdb:jedb)) ; cs%Idx2dyCv(:,:) = 0.0
1283  alloc_(cs%Idxdy2u(isdb:iedb,jsd:jed)) ; cs%Idxdy2u(:,:) = 0.0
1284  alloc_(cs%Idxdy2v(isd:ied,jsdb:jedb)) ; cs%Idxdy2v(:,:) = 0.0
1285 
1286  alloc_(cs%Ah_bg_xx(isd:ied,jsd:jed)) ; cs%Ah_bg_xx(:,:) = 0.0
1287  alloc_(cs%Ah_bg_xy(isdb:iedb,jsdb:jedb)) ; cs%Ah_bg_xy(:,:) = 0.0
1288  if (cs%bound_Ah .or. cs%better_bound_Ah) then
1289  alloc_(cs%Ah_Max_xx(isd:ied,jsd:jed)) ; cs%Ah_Max_xx(:,:) = 0.0
1290  alloc_(cs%Ah_Max_xy(isdb:iedb,jsdb:jedb)) ; cs%Ah_Max_xy(:,:) = 0.0
1291  endif
1292  if (cs%Smagorinsky_Ah) then
1293  alloc_(cs%Biharm_Const_xx(isd:ied,jsd:jed)) ; cs%Biharm_Const_xx(:,:) = 0.0
1294  alloc_(cs%Biharm_Const_xy(isdb:iedb,jsdb:jedb)) ; cs%Biharm_Const_xy(:,:) = 0.0
1295  if (cs%bound_Coriolis) then
1296  alloc_(cs%Biharm_Const2_xx(isd:ied,jsd:jed)) ; cs%Biharm_Const2_xx(:,:) = 0.0
1297  alloc_(cs%Biharm_Const2_xy(isdb:iedb,jsdb:jedb)) ; cs%Biharm_Const2_xy(:,:) = 0.0
1298  endif
1299  endif
1300  if (cs%Leith_Ah) then
1301  alloc_(cs%Biharm5_Const_xx(isd:ied,jsd:jed)) ; cs%Biharm5_Const_xx(:,:) = 0.0
1302  alloc_(cs%Biharm5_Const_xy(isdb:iedb,jsdb:jedb)) ; cs%Biharm5_Const_xy(:,:) = 0.0
1303  endif
1304  endif
1305 
1306  do j=js-2,jeq+1 ; do i=is-2,ieq+1
1307  cs%DX2q(i,j) = g%dxBu(i,j)*g%dxBu(i,j) ; cs%DY2q(i,j) = g%dyBu(i,j)*g%dyBu(i,j)
1308  cs%DX_dyBu(i,j) = g%dxBu(i,j)*g%IdyBu(i,j) ; cs%DY_dxBu(i,j) = g%dyBu(i,j)*g%IdxBu(i,j)
1309  enddo ; enddo
1310  do j=jsq-1,jeq+2 ; do i=isq-1,ieq+2
1311  cs%DX2h(i,j) = g%dxT(i,j)*g%dxT(i,j) ; cs%DY2h(i,j) = g%dyT(i,j)*g%dyT(i,j)
1312  cs%DX_dyT(i,j) = g%dxT(i,j)*g%IdyT(i,j) ; cs%DY_dxT(i,j) = g%dyT(i,j)*g%IdxT(i,j)
1313  enddo ; enddo
1314 
1315  do j=jsq,jeq+1 ; do i=isq,ieq+1
1316  cs%reduction_xx(i,j) = 1.0
1317  if ((g%dy_Cu(i,j) > 0.0) .and. (g%dy_Cu(i,j) < g%dyCu(i,j)) .and. &
1318  (g%dy_Cu(i,j) < g%dyCu(i,j) * cs%reduction_xx(i,j))) &
1319  cs%reduction_xx(i,j) = g%dy_Cu(i,j) / g%dyCu(i,j)
1320  if ((g%dy_Cu(i-1,j) > 0.0) .and. (g%dy_Cu(i-1,j) < g%dyCu(i-1,j)) .and. &
1321  (g%dy_Cu(i-1,j) < g%dyCu(i-1,j) * cs%reduction_xx(i,j))) &
1322  cs%reduction_xx(i,j) = g%dy_Cu(i-1,j) / g%dyCu(i-1,j)
1323  if ((g%dx_Cv(i,j) > 0.0) .and. (g%dx_Cv(i,j) < g%dxCv(i,j)) .and. &
1324  (g%dx_Cv(i,j) < g%dxCv(i,j) * cs%reduction_xx(i,j))) &
1325  cs%reduction_xx(i,j) = g%dx_Cv(i,j) / g%dxCv(i,j)
1326  if ((g%dx_Cv(i,j-1) > 0.0) .and. (g%dx_Cv(i,j-1) < g%dxCv(i,j-1)) .and. &
1327  (g%dx_Cv(i,j-1) < g%dxCv(i,j-1) * cs%reduction_xx(i,j))) &
1328  cs%reduction_xx(i,j) = g%dx_Cv(i,j-1) / g%dxCv(i,j-1)
1329  enddo ; enddo
1330 
1331  do j=js-1,jeq ; do i=is-1,ieq
1332  cs%reduction_xy(i,j) = 1.0
1333  if ((g%dy_Cu(i,j) > 0.0) .and. (g%dy_Cu(i,j) < g%dyCu(i,j)) .and. &
1334  (g%dy_Cu(i,j) < g%dyCu(i,j) * cs%reduction_xy(i,j))) &
1335  cs%reduction_xy(i,j) = g%dy_Cu(i,j) / g%dyCu(i,j)
1336  if ((g%dy_Cu(i,j+1) > 0.0) .and. (g%dy_Cu(i,j+1) < g%dyCu(i,j+1)) .and. &
1337  (g%dy_Cu(i,j+1) < g%dyCu(i,j+1) * cs%reduction_xy(i,j))) &
1338  cs%reduction_xy(i,j) = g%dy_Cu(i,j+1) / g%dyCu(i,j+1)
1339  if ((g%dx_Cv(i,j) > 0.0) .and. (g%dx_Cv(i,j) < g%dxCv(i,j)) .and. &
1340  (g%dx_Cv(i,j) < g%dxCv(i,j) * cs%reduction_xy(i,j))) &
1341  cs%reduction_xy(i,j) = g%dx_Cv(i,j) / g%dxCv(i,j)
1342  if ((g%dx_Cv(i+1,j) > 0.0) .and. (g%dx_Cv(i+1,j) < g%dxCv(i+1,j)) .and. &
1343  (g%dx_Cv(i+1,j) < g%dxCv(i+1,j) * cs%reduction_xy(i,j))) &
1344  cs%reduction_xy(i,j) = g%dx_Cv(i+1,j) / g%dxCv(i+1,j)
1345  enddo ; enddo
1346 
1347  if (cs%Laplacian) then
1348  ! The 0.3 below was 0.4 in MOM1.10. The change in hq requires
1349  ! this to be less than 1/3, rather than 1/2 as before.
1350  if (cs%bound_Kh .or. cs%bound_Ah) kh_limit = 0.3 / (dt*4.0)
1351  do j=jsq,jeq+1 ; do i=isq,ieq+1
1352  grid_sp_h2 = (2.0*cs%DX2h(i,j)*cs%DY2h(i,j)) / (cs%DX2h(i,j) + cs%DY2h(i,j))
1353  grid_sp_h3 = grid_sp_h2*sqrt(grid_sp_h2)
1354  if (cs%Smagorinsky_Kh) cs%LAPLAC_CONST_xx(i,j) = smag_lap_const * grid_sp_h2
1355  if (cs%Leith_Kh) cs%LAPLAC3_CONST_xx(i,j) = leith_lap_const * grid_sp_h3
1356 
1357  cs%Kh_bg_xx(i,j) = max(kh, kh_vel_scale * sqrt(grid_sp_h2))
1358 
1359  if (cs%use_Kh_bg_2d) cs%Kh_bg_xx(i,j) = max(cs%Kh_bg_2d(i,j), cs%Kh_bg_xx(i,j))
1360 
1361  if (cs%bound_Kh .and. .not.cs%better_bound_Kh) then
1362  cs%Kh_Max_xx(i,j) = kh_limit * grid_sp_h2
1363  cs%Kh_bg_xx(i,j) = min(cs%Kh_bg_xx(i,j), cs%Kh_Max_xx(i,j))
1364  endif
1365  enddo ; enddo
1366  do j=js-1,jeq ; do i=is-1,ieq
1367  grid_sp_q2 = (2.0*cs%DX2q(i,j)*cs%DY2q(i,j)) / (cs%DX2q(i,j) + cs%DY2q(i,j))
1368  grid_sp_q3 = grid_sp_q2*sqrt(grid_sp_q2)
1369  if (cs%Smagorinsky_Kh) cs%LAPLAC_CONST_xy(i,j) = smag_lap_const * grid_sp_q2
1370  if (cs%Leith_Kh) cs%LAPLAC3_CONST_xy(i,j) = leith_lap_const * grid_sp_q3
1371 
1372  cs%Kh_bg_xy(i,j) = max(kh, kh_vel_scale * sqrt(grid_sp_q2))
1373 
1374  if (cs%use_Kh_bg_2d) cs%Kh_bg_xy(i,j) = max(cs%Kh_bg_2d(i,j), cs%Kh_bg_xy(i,j))
1375 
1376  if (cs%bound_Kh .and. .not.cs%better_bound_Kh) then
1377  cs%Kh_Max_xy(i,j) = kh_limit * grid_sp_q2
1378  cs%Kh_bg_xy(i,j) = min(cs%Kh_bg_xy(i,j), cs%Kh_Max_xy(i,j))
1379  endif
1380  enddo ; enddo
1381  endif
1382 
1383  if (cs%biharmonic) then
1384 
1385  do j=js-1,jeq+1 ; do i=isq-1,ieq+1
1386  cs%IDX2dyCu(i,j) = (g%IdxCu(i,j)*g%IdxCu(i,j)) * g%IdyCu(i,j)
1387  cs%IDXDY2u(i,j) = g%IdxCu(i,j) * (g%IdyCu(i,j)*g%IdyCu(i,j))
1388  enddo ; enddo
1389  do j=jsq-1,jeq+1 ; do i=is-1,ieq+1
1390  cs%IDX2dyCv(i,j) = (g%IdxCv(i,j)*g%IdxCv(i,j)) * g%IdyCv(i,j)
1391  cs%IDXDY2v(i,j) = g%IdxCv(i,j) * (g%IdyCv(i,j)*g%IdyCv(i,j))
1392  enddo ; enddo
1393 
1394  cs%Ah_bg_xy(:,:) = 0.0
1395  ! The 0.3 below was 0.4 in MOM1.10. The change in hq requires
1396  ! this to be less than 1/3, rather than 1/2 as before.
1397  ah_limit = 0.3 / (dt*64.0)
1398  if (cs%Smagorinsky_Ah .and. cs%bound_Coriolis) &
1399  boundcorconst = 1.0 / (5.0*(bound_cor_vel*bound_cor_vel))
1400  do j=jsq,jeq+1 ; do i=isq,ieq+1
1401  grid_sp_h2 = (2.0*cs%DX2h(i,j)*cs%DY2h(i,j)) / (cs%DX2h(i,j)+cs%DY2h(i,j))
1402  grid_sp_h3 = grid_sp_h2*sqrt(grid_sp_h2)
1403 
1404  if (cs%Smagorinsky_Ah) then
1405  cs%BIHARM_CONST_xx(i,j) = smag_bi_const * (grid_sp_h2 * grid_sp_h2)
1406  if (cs%bound_Coriolis) then
1407  fmax = max(abs(g%CoriolisBu(i-1,j-1)), abs(g%CoriolisBu(i,j-1)), &
1408  abs(g%CoriolisBu(i-1,j)), abs(g%CoriolisBu(i,j)))
1409  cs%Biharm_Const2_xx(i,j) = (grid_sp_h2 * grid_sp_h2 * grid_sp_h2) * &
1410  (fmax * boundcorconst)
1411  endif
1412  endif
1413  if (cs%Leith_Ah) then
1414  cs%BIHARM5_CONST_xx(i,j) = leith_bi_const * (grid_sp_h2 * grid_sp_h3)
1415  endif
1416 
1417  cs%Ah_bg_xx(i,j) = max(ah, ah_vel_scale * grid_sp_h2 * sqrt(grid_sp_h2))
1418  if (cs%bound_Ah .and. .not.cs%better_bound_Ah) then
1419  cs%Ah_Max_xx(i,j) = ah_limit * (grid_sp_h2 * grid_sp_h2)
1420  cs%Ah_bg_xx(i,j) = min(cs%Ah_bg_xx(i,j), cs%Ah_Max_xx(i,j))
1421  endif
1422  enddo ; enddo
1423  do j=js-1,jeq ; do i=is-1,ieq
1424  grid_sp_q2 = (2.0*cs%DX2q(i,j)*cs%DY2q(i,j)) / (cs%DX2q(i,j)+cs%DY2q(i,j))
1425  grid_sp_q3 = grid_sp_q2*sqrt(grid_sp_q2)
1426 
1427  if (cs%Smagorinsky_Ah) then
1428  cs%BIHARM_CONST_xy(i,j) = smag_bi_const * (grid_sp_q2 * grid_sp_q2)
1429  if (cs%bound_Coriolis) then
1430  cs%Biharm_Const2_xy(i,j) = (grid_sp_q2 * grid_sp_q2 * grid_sp_q2) * &
1431  (abs(g%CoriolisBu(i,j)) * boundcorconst)
1432  endif
1433  endif
1434 
1435  if (cs%Leith_Ah) then
1436  cs%BIHARM5_CONST_xy(i,j) = leith_bi_const * (grid_sp_q2 * grid_sp_q3)
1437  endif
1438 
1439  cs%Ah_bg_xy(i,j) = max(ah, ah_vel_scale * grid_sp_q2 * sqrt(grid_sp_q2))
1440  if (cs%bound_Ah .and. .not.cs%better_bound_Ah) then
1441  cs%Ah_Max_xy(i,j) = ah_limit * (grid_sp_q2 * grid_sp_q2)
1442  cs%Ah_bg_xy(i,j) = min(cs%Ah_bg_xy(i,j), cs%Ah_Max_xy(i,j))
1443  endif
1444  enddo ; enddo
1445  endif
1446 
1447  ! The Laplacian bounds should avoid overshoots when CS%bound_coef < 1.
1448  if (cs%Laplacian .and. cs%better_bound_Kh) then
1449  idt = 1.0 / dt
1450  do j=jsq,jeq+1 ; do i=isq,ieq+1
1451  denom = max( &
1452  (cs%DY2h(i,j) * cs%DY_dxT(i,j) * (g%IdyCu(i,j) + g%IdyCu(i-1,j)) * &
1453  max(g%IdyCu(i,j)*g%IareaCu(i,j), g%IdyCu(i-1,j)*g%IareaCu(i-1,j)) ), &
1454  (cs%DX2h(i,j) * cs%DX_dyT(i,j) * (g%IdxCv(i,j) + g%IdxCv(i,j-1)) * &
1455  max(g%IdxCv(i,j)*g%IareaCv(i,j), g%IdxCv(i,j-1)*g%IareaCv(i,j-1)) ) )
1456  cs%Kh_Max_xx(i,j) = 0.0
1457  if (denom > 0.0) &
1458  cs%Kh_Max_xx(i,j) = cs%bound_coef * 0.25 * idt / denom
1459  enddo ; enddo
1460  do j=js-1,jeq ; do i=is-1,ieq
1461  denom = max( &
1462  (cs%DX2q(i,j) * cs%DX_dyBu(i,j) * (g%IdxCu(i,j+1) + g%IdxCu(i,j)) * &
1463  max(g%IdxCu(i,j)*g%IareaCu(i,j), g%IdxCu(i,j+1)*g%IareaCu(i,j+1)) ), &
1464  (cs%DY2q(i,j) * cs%DY_dxBu(i,j) * (g%IdyCv(i+1,j) + g%IdyCv(i,j)) * &
1465  max(g%IdyCv(i,j)*g%IareaCv(i,j), g%IdyCv(i+1,j)*g%IareaCv(i+1,j)) ) )
1466  cs%Kh_Max_xy(i,j) = 0.0
1467  if (denom > 0.0) &
1468  cs%Kh_Max_xy(i,j) = cs%bound_coef * 0.25 * idt / denom
1469  enddo ; enddo
1470  endif
1471 
1472  ! The biharmonic bounds should avoid overshoots when CS%bound_coef < 0.5, but
1473  ! empirically work for CS%bound_coef <~ 1.0
1474  if (cs%biharmonic .and. cs%better_bound_Ah) then
1475  idt = 1.0 / dt
1476  do j=js-1,jeq+1 ; do i=isq-1,ieq+1
1477  u0u(i,j) = cs%IDXDY2u(i,j)*(cs%DY2h(i+1,j)*cs%DY_dxT(i+1,j)*(g%IdyCu(i+1,j) + g%IdyCu(i,j)) + &
1478  cs%DY2h(i,j) * cs%DY_dxT(i,j) * (g%IdyCu(i,j) + g%IdyCu(i-1,j)) ) + &
1479  cs%IDX2dyCu(i,j)*(cs%DX2q(i,j) * cs%DX_dyBu(i,j) * (g%IdxCu(i,j+1) + g%IdxCu(i,j)) + &
1480  cs%DX2q(i,j-1)*cs%DX_dyBu(i,j-1)*(g%IdxCu(i,j) + g%IdxCu(i,j-1)) )
1481 
1482  u0v(i,j) = cs%IDXDY2u(i,j)*(cs%DY2h(i+1,j)*cs%DX_dyT(i+1,j)*(g%IdxCv(i+1,j) + g%IdxCv(i+1,j-1)) + &
1483  cs%DY2h(i,j) * cs%DX_dyT(i,j) * (g%IdxCv(i,j) + g%IdxCv(i,j-1)) ) + &
1484  cs%IDX2dyCu(i,j)*(cs%DX2q(i,j) * cs%DY_dxBu(i,j) * (g%IdyCv(i+1,j) + g%IdyCv(i,j)) + &
1485  cs%DX2q(i,j-1)*cs%DY_dxBu(i,j-1)*(g%IdyCv(i+1,j-1) + g%IdyCv(i,j-1)) )
1486  enddo ; enddo
1487  do j=jsq-1,jeq+1 ; do i=is-1,ieq+1
1488  v0u(i,j) = cs%IDXDY2v(i,j)*(cs%DY2q(i,j) * cs%DX_dyBu(i,j) * (g%IdxCu(i,j+1) + g%IdxCu(i,j)) + &
1489  cs%DY2q(i-1,j)*cs%DX_dyBu(i-1,j)*(g%IdxCu(i-1,j+1) + g%IdxCu(i-1,j)) ) + &
1490  cs%IDX2dyCv(i,j)*(cs%DX2h(i,j+1)*cs%DY_dxT(i,j+1)*(g%IdyCu(i,j+1) + g%IdyCu(i-1,j+1)) + &
1491  cs%DX2h(i,j) * cs%DY_dxT(i,j) * (g%IdyCu(i,j) + g%IdyCu(i-1,j)) )
1492 
1493  v0v(i,j) = cs%IDXDY2v(i,j)*(cs%DY2q(i,j) * cs%DY_dxBu(i,j) * (g%IdyCv(i+1,j) + g%IdyCv(i,j)) + &
1494  cs%DY2q(i-1,j)*cs%DY_dxBu(i-1,j)*(g%IdyCv(i,j) + g%IdyCv(i-1,j)) ) + &
1495  cs%IDX2dyCv(i,j)*(cs%DX2h(i,j+1)*cs%DX_dyT(i,j+1)*(g%IdxCv(i,j+1) + g%IdxCv(i,j)) + &
1496  cs%DX2h(i,j) * cs%DX_dyT(i,j) * (g%IdxCv(i,j) + g%IdxCv(i,j-1)) )
1497  enddo ; enddo
1498 
1499  do j=jsq,jeq+1 ; do i=isq,ieq+1
1500  denom = max( &
1501  (cs%DY2h(i,j) * &
1502  (cs%DY_dxT(i,j)*(g%IdyCu(i,j)*u0u(i,j) + g%IdyCu(i-1,j)*u0u(i-1,j)) + &
1503  cs%DX_dyT(i,j)*(g%IdxCv(i,j)*v0u(i,j) + g%IdxCv(i,j-1)*v0u(i,j-1))) * &
1504  max(g%IdyCu(i,j)*g%IareaCu(i,j), g%IdyCu(i-1,j)*g%IareaCu(i-1,j)) ), &
1505  (cs%DX2h(i,j) * &
1506  (cs%DY_dxT(i,j)*(g%IdyCu(i,j)*u0v(i,j) + g%IdyCu(i-1,j)*u0v(i-1,j)) + &
1507  cs%DX_dyT(i,j)*(g%IdxCv(i,j)*v0v(i,j) + g%IdxCv(i,j-1)*v0v(i,j-1))) * &
1508  max(g%IdxCv(i,j)*g%IareaCv(i,j), g%IdxCv(i,j-1)*g%IareaCv(i,j-1)) ) )
1509  cs%Ah_Max_xx(i,j) = 0.0
1510  if (denom > 0.0) &
1511  cs%Ah_Max_xx(i,j) = cs%bound_coef * 0.5 * idt / denom
1512  enddo ; enddo
1513 
1514  do j=js-1,jeq ; do i=is-1,ieq
1515  denom = max( &
1516  (cs%DX2q(i,j) * &
1517  (cs%DX_dyBu(i,j)*(u0u(i,j+1)*g%IdxCu(i,j+1) + u0u(i,j)*g%IdxCu(i,j)) + &
1518  cs%DY_dxBu(i,j)*(v0u(i+1,j)*g%IdyCv(i+1,j) + v0u(i,j)*g%IdyCv(i,j))) * &
1519  max(g%IdxCu(i,j)*g%IareaCu(i,j), g%IdxCu(i,j+1)*g%IareaCu(i,j+1)) ), &
1520  (cs%DY2q(i,j) * &
1521  (cs%DX_dyBu(i,j)*(u0v(i,j+1)*g%IdxCu(i,j+1) + u0v(i,j)*g%IdxCu(i,j)) + &
1522  cs%DY_dxBu(i,j)*(v0v(i+1,j)*g%IdyCv(i+1,j) + v0v(i,j)*g%IdyCv(i,j))) * &
1523  max(g%IdyCv(i,j)*g%IareaCv(i,j), g%IdyCv(i+1,j)*g%IareaCv(i+1,j)) ) )
1524  cs%Ah_Max_xy(i,j) = 0.0
1525  if (denom > 0.0) &
1526  cs%Ah_Max_xy(i,j) = cs%bound_coef * 0.5 * idt / denom
1527  enddo ; enddo
1528  endif
1529 
1530  ! Register fields for output from this module.
1531 
1532  cs%id_diffu = register_diag_field('ocean_model', 'diffu', diag%axesCuL, time, &
1533  'Zonal Acceleration from Horizontal Viscosity', 'meter second-2')
1534 
1535  cs%id_diffv = register_diag_field('ocean_model', 'diffv', diag%axesCvL, time, &
1536  'Meridional Acceleration from Horizontal Viscosity', 'meter second-2')
1537 
1538  if (cs%biharmonic) then
1539  cs%id_Ah_h = register_diag_field('ocean_model', 'Ahh', diag%axesTL, time, &
1540  'Biharmonic Horizontal Viscosity at h Points', 'meter4 second-1', &
1541  cmor_field_name='difmxybo', cmor_units='m4 s-1', &
1542  cmor_long_name='Ocean lateral biharmonic viscosity', &
1543  cmor_standard_name='ocean_momentum_xy_biharmonic_diffusivity')
1544 
1545  cs%id_Ah_q = register_diag_field('ocean_model', 'Ahq', diag%axesBL, time, &
1546  'Biharmonic Horizontal Viscosity at q Points', 'meter4 second-1')
1547  endif
1548 
1549  if (cs%Laplacian) then
1550  cs%id_Kh_h = register_diag_field('ocean_model', 'Khh', diag%axesTL, time, &
1551  'Laplacian Horizontal Viscosity at h Points', 'meter2 second-1', &
1552  cmor_field_name='difmxylo', cmor_units='m2 s-1', &
1553  cmor_long_name='Ocean lateral Laplacian viscosity', &
1554  cmor_standard_name='ocean_momentum_xy_laplacian_diffusivity')
1555 
1556  cs%id_Kh_q = register_diag_field('ocean_model', 'Khq', diag%axesBL, time, &
1557  'Laplacian Horizontal Viscosity at q Points', 'meter2 second-1')
1558  endif
1559 
1560  cs%id_FrictWork = register_diag_field('ocean_model','FrictWork',diag%axesTL,time,&
1561  'Integral work done by lateral friction terms', 'Watt meter-2')
1562 
1563  cs%id_FrictWorkIntz = register_diag_field('ocean_model','FrictWorkIntz',diag%axesT1,time, &
1564  'Depth integrated work done by lateral friction', 'Watt meter-2', &
1565  cmor_field_name='dispkexyfo', cmor_units='W m-2', &
1566  cmor_long_name='Depth integrated ocean kinetic energy dissipation due to lateral friction',&
1567  cmor_standard_name='ocean_kinetic_energy_dissipation_per_unit_area_due_to_xy_friction')
1568 
1569  if (cs%Laplacian .or. get_all) then
1570  endif
1571 
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:19
Here is the call graph for this function:

◆ horizontal_viscosity()

subroutine, public mom_hor_visc::horizontal_viscosity ( real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(in)  u,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(in)  v,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h,
real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(out)  diffu,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(out)  diffv,
type(meke_type), pointer  MEKE,
type(varmix_cs), pointer  VarMix,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(hor_visc_cs), pointer  CS,
type(ocean_obc_type), optional, pointer  OBC 
)

This subroutine determines the acceleration due to the horizontal viscosity. A combination of biharmonic and Laplacian forms can be used. The coefficient may either be a constant or a shear-dependent form. The biharmonic is determined by twice taking the divergence of an appropriately defined stress tensor. The Laplacian is determined by doing so once. To work, the following fields must be set outside of the usual is to ie range before this subroutine is called: v[is-2,is-1,ie+1,ie+2], u[is-2,is-1,ie+1,ie+2], and h[is-1,ie+1], with a similarly sized halo in the y-direction.

Parameters
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]uThe zonal velocity, in m s-1.
[in]vThe meridional velocity, in m s-1.
[in]hLayer thicknesses, in H
[out]diffuZonal acceleration due to convergence of
[out]diffvMeridional acceleration due to convergence
mekePointer to a structure containing fields related to Mesoscale Eddy Kinetic Energy.
varmixPointer to a structure with fields that specify the spatially variable viscosities
csPontrol structure returned by a previous call to hor_visc_init.
obcPointer to an open boundary condition type

Definition at line 232 of file MOM_hor_visc.F90.

References mom_error_handler::mom_error(), mom_open_boundary::obc_direction_n, and mom_open_boundary::obc_direction_s.

Referenced by mom_dynamics_split_rk2::step_mom_dyn_split_rk2(), mom_dynamics_unsplit::step_mom_dyn_unsplit(), and mom_dynamics_unsplit_rk2::step_mom_dyn_unsplit_rk2().

232  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
233  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
234  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
235  intent(in) :: u !< The zonal velocity, in m s-1.
236  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
237  intent(in) :: v !< The meridional velocity, in m s-1.
238  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
239  intent(in) :: h !< Layer thicknesses, in H
240  !! (usually m or kg m-2).
241  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
242  intent(out) :: diffu !< Zonal acceleration due to convergence of
243  !! along-coordinate stress tensor (m/s2)
244  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
245  intent(out) :: diffv !< Meridional acceleration due to convergence
246  !! of along-coordinate stress tensor (m/s2).
247  type(meke_type), pointer :: meke !< Pointer to a structure containing fields
248  !! related to Mesoscale Eddy Kinetic Energy.
249  type(varmix_cs), pointer :: varmix !< Pointer to a structure with fields that
250  !! specify the spatially variable viscosities
251  type(hor_visc_cs), pointer :: cs !< Pontrol structure returned by a previous
252  !! call to hor_visc_init.
253  type(ocean_obc_type), pointer, optional :: obc !< Pointer to an open boundary condition type
254 
255 ! Arguments:
256 ! (in) u - zonal velocity (m/s)
257 ! (in) v - meridional velocity (m/s)
258 ! (in) h - layer thickness (m or kg m-2); h units are referred to as H.
259 ! (out) diffu - zonal acceleration due to convergence of
260 ! along-coordinate stress tensor (m/s2)
261 ! (out) diffv - meridional acceleration due to convergence of
262 ! along-coordinate stress tensor (m/s2)
263 ! (inout) MEKE - pointer to a structure containing fields related to
264 ! Mesoscale Eddy Kinetic Energy
265 ! (in) VarMix - pointer to a structure with fields that specify the
266 ! spatially variable viscosities
267 ! (in) G - ocean grid structure
268 ! (in) GV - The ocean's vertical grid structure.
269 ! (in) CS - control structure returned by a previous call to
270 ! hor_visc_init
271 ! (in) OBC - pointer to an open boundary condition type
272 
273 ! By R. Hallberg, August 1998 - November 1998.
274 ! This subroutine determines the acceleration due to the
275 ! horizontal viscosity. A combination of biharmonic and Laplacian
276 ! forms can be used. The coefficient may either be a constant or
277 ! a shear-dependent form. The biharmonic is determined by twice
278 ! taking the divergence of an appropriately defined stress tensor.
279 ! The Laplacian is determined by doing so once.
280 ! To work, the following fields must be set outside of the usual
281 ! is to ie range before this subroutine is called:
282 ! v[is-2,is-1,ie+1,ie+2], u[is-2,is-1,ie+1,ie+2], and h[is-1,ie+1],
283 ! with a similarly sized halo in the y-direction.
284 
285  real, dimension(SZIB_(G),SZJ_(G)) :: &
286  u0, & ! Laplacian of u (m-1 s-1)
287  h_u ! Thickness interpolated to u points, in H.
288  real, dimension(SZI_(G),SZJB_(G)) :: &
289  v0, & ! Laplacian of v (m-1 s-1)
290  h_v ! Thickness interpolated to v points, in H.
291 
292  real, dimension(SZI_(G),SZJ_(G)) :: &
293  sh_xx, & ! horizontal tension (du/dx - dv/dy) (1/sec) including metric terms
294  str_xx,& ! str_xx is the diagonal term in the stress tensor (H m2 s-2)
295  bhstr_xx,& ! A copy of str_xx that only contains the biharmonic contribution (H m2 s-2)
296  div_xx, & ! horizontal divergence (du/dx + dv/dy) (1/sec) including metric terms
297  frictworkintz ! depth integrated energy dissipated by lateral friction (W/m2)
298 
299  real, dimension(SZIB_(G),SZJB_(G)) :: &
300  dvdx, dudy, & ! components in the shearing strain (s-1)
301  sh_xy, & ! horizontal shearing strain (du/dy + dv/dx) (1/sec) including metric terms
302  str_xy, & ! str_xy is the cross term in the stress tensor (H m2 s-2)
303  bhstr_xy, & ! A copy of str_xy that only contains the biharmonic contribution (H m2 s-2)
304  vort_xy ! vertical vorticity (dv/dx - du/dy) (1/sec) including metric terms
305 
306  real, dimension(SZI_(G),SZJB_(G)) :: &
307  vort_xy_dx, & ! x-derivative of vertical vorticity (d/dx(dv/dx - du/dy)) (m-1 sec-1) including metric terms
308  div_xx_dy ! y-derivative of horizontal divergence (d/dy(du/dx + dv/dy)) (m-1 sec-1) including metric terms
309 
310  real, dimension(SZIB_(G),SZJ_(G)) :: &
311  vort_xy_dy, & ! y-derivative of vertical vorticity (d/dy(dv/dx - du/dy)) (m-1 sec-1) including metric terms
312  div_xx_dx ! x-derivative of horizontal divergence (d/dx(du/dx + dv/dy)) (m-1 sec-1) including metric terms
313 
314  real, dimension(SZIB_(G),SZJB_(G),SZK_(G)) :: &
315  ah_q, & ! biharmonic viscosity at corner points (m4/s)
316  kh_q ! Laplacian viscosity at corner points (m2/s)
317 
318  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
319  ah_h, & ! biharmonic viscosity at thickness points (m4/s)
320  kh_h, & ! Laplacian viscosity at thickness points (m2/s)
321  frictwork ! energy dissipated by lateral friction (W/m2)
322 
323  real :: ah ! biharmonic viscosity (m4/s)
324  real :: kh ! Laplacian viscosity (m2/s)
325  real :: ahsm ! Smagorinsky biharmonic viscosity (m4/s)
326  real :: khsm ! Smagorinsky Laplacian viscosity (m2/s)
327  real :: ahlth ! 2D Leith biharmonic viscosity (m4/s)
328  real :: khlth ! 2D Leith Laplacian viscosity (m2/s)
329  real :: mod_leith ! nondimensional coefficient for divergence part of modified Leith
330  ! viscosity. Here set equal to nondimensional Laplacian Leith constant.
331  ! This is set equal to zero if modified Leith is not used.
332  real :: shear_mag ! magnitude of the shear (1/s)
333  real :: vort_mag ! magnitude of the vorticity (1/s)
334  real :: h2uq, h2vq ! temporary variables in units of H^2 (i.e. m2 or kg2 m-4).
335  real :: hu, hv ! Thicknesses interpolated by arithmetic means to corner
336  ! points; these are first interpolated to u or v velocity
337  ! points where masks are applied, in units of H (i.e. m or kg m-2).
338  real :: hq ! harmonic mean of the harmonic means of the u- & v-
339  ! point thicknesses, in H; This form guarantees that hq/hu < 4.
340  real :: h_neglect ! thickness so small it can be lost in roundoff and so neglected (H)
341  real :: h_neglect3 ! h_neglect^3, in H3
342  real :: hrat_min ! minimum thicknesses at the 4 neighboring
343  ! velocity points divided by the thickness at the stress
344  ! point (h or q point) (nondimensional)
345  real :: visc_bound_rem ! fraction of overall viscous bounds that
346  ! remain to be applied (nondim)
347  real :: kh_scale ! A factor between 0 and 1 by which the horizontal
348  ! Laplacian viscosity is rescaled
349  real :: roscl ! The scaling function for MEKE source term
350  real :: fath ! abs(f) at h-point for MEKE source term (s-1)
351 
352  logical :: rescale_kh
353  logical :: find_frictwork
354  logical :: apply_obc = .false.
355  logical :: use_meke_ku
356  integer :: is, ie, js, je, isq, ieq, jsq, jeq, nz
357  integer :: i, j, k, n
358 
359  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
360  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
361 
362  h_neglect = gv%H_subroundoff
363  h_neglect3 = h_neglect**3
364 
365  if (present(obc)) then ; if (associated(obc)) then ; if (obc%OBC_pe) then
366  apply_obc = obc%Flather_u_BCs_exist_globally .or. obc%Flather_v_BCs_exist_globally
367  apply_obc = .true.
368  endif ; endif ; endif
369 
370  if (.not.associated(cs)) call mom_error(fatal, &
371  "MOM_hor_visc: Module must be initialized before it is used.")
372  if (.not.(cs%Laplacian .or. cs%biharmonic)) return
373 
374  find_frictwork = (cs%id_FrictWork > 0)
375  if (cs%id_FrictWorkIntz > 0) find_frictwork = .true.
376  if (associated(meke)) then
377  if (associated(meke%mom_src)) find_frictwork = .true.
378  endif
379 
380  rescale_kh = .false.
381  if (associated(varmix)) then
382  rescale_kh = varmix%Resoln_scaled_Kh
383  if (rescale_kh .and. &
384  (.not.associated(varmix%Res_fn_h) .or. .not.associated(varmix%Res_fn_q))) &
385  call mom_error(fatal, "MOM_hor_visc: VarMix%Res_fn_h and " //&
386  "VarMix%Res_fn_q both need to be associated with Resoln_scaled_Kh.")
387  endif
388 
389  ! Toggle whether to use a Laplacian viscosity derived from MEKE
390  use_meke_ku = associated(meke%Ku)
391 
392 !$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,CS,G,GV,u,v,is,js,ie,je,h, &
393 !$OMP rescale_Kh,VarMix,h_neglect,h_neglect3, &
394 !$OMP Kh_h,Ah_h,Kh_q,Ah_q,diffu,apply_OBC,OBC,diffv, &
395 !$OMP find_FrictWork,FrictWork,use_MEKE_Ku,MEKE) &
396 !$OMP private(u0, v0, sh_xx, str_xx, visc_bound_rem, &
397 !$OMP sh_xy, str_xy, Ah, Kh, AhSm, KhSm, dvdx, dudy, &
398 !$OMP bhstr_xx, bhstr_xy,FatH,RoScl, hu, hv, h_u, h_v, &
399 !$OMP vort_xy,vort_xy_dx,vort_xy_dy,Vort_mag,AhLth,KhLth, &
400 !$OMP div_xx, div_xx_dx, div_xx_dy, mod_Leith, &
401 !$OMP Shear_mag, h2uq, h2vq, hq, Kh_scale, hrat_min)
402  do k=1,nz
403 
404 ! This code uses boundary conditions that are consistent with
405 ! free slip and no normal flow boundary conditions. The boundary
406 ! conditions for the western boundary, for example, are:
407 ! dv/dx = 0, d^3v/dx^3 = 0, u = 0, d^2u/dx^2 = 0 .
408 ! The overall scheme is second order accurate.
409 ! All of the metric terms are retained, and the repeated use of
410 ! the symmetric stress tensor insures that no stress is applied with
411 ! no flow or solid-body rotation, even with non-constant values of
412 ! of the biharmonic viscosity.
413 
414 ! The following are the forms of the horizontal tension and hori-
415 ! shearing strain advocated by Smagorinsky (1993) and discussed in
416 ! Griffies and Hallberg (MWR, 2000).
417  do j=jsq-1,jeq+2 ; do i=isq-1,ieq+2
418  sh_xx(i,j) = (cs%DY_dxT(i,j)*(g%IdyCu(i,j) * u(i,j,k) - &
419  g%IdyCu(i-1,j) * u(i-1,j,k)) - &
420  cs%DX_dyT(i,j)*(g%IdxCv(i,j) * v(i,j,k) - &
421  g%IdxCv(i,j-1)*v(i,j-1,k)))
422  div_xx(i,j) = 0.5*((g%dyCu(i,j) * u(i,j,k) * (h(i+1,j,k)+h(i,j,k)) - &
423  g%dyCu(i-1,j) * u(i-1,j,k) * (h(i-1,j,k)+h(i,j,k)) ) + &
424  (g%dxCv(i,j) * v(i,j,k) * (h(i,j,k)+h(i,j+1,k)) - &
425  g%dxCv(i,j-1)*v(i,j-1,k)*(h(i,j,k)+h(i,j-1,k))))*g%IareaT(i,j)/ &
426  (h(i,j,k) + h_neglect)
427  enddo ; enddo
428  ! Components for the shearing strain
429  do j=js-2,jeq+1 ; do i=is-2,ieq+1
430  dvdx(i,j) = cs%DY_dxBu(i,j)*(v(i+1,j,k)*g%IdyCv(i+1,j) - v(i,j,k)*g%IdyCv(i,j))
431  dudy(i,j) = cs%DX_dyBu(i,j)*(u(i,j+1,k)*g%IdxCu(i,j+1) - u(i,j,k)*g%IdxCu(i,j))
432  enddo ; enddo
433 
434  ! Interpolate the thicknesses to velocity points.
435  ! The extra wide halos are to accomodate the cross-corner-point projections
436  ! in OBCs, which are not ordinarily be necessary, and might not be necessary
437  ! even with OBCs if the accelerations are zeroed at OBC points, in which
438  ! case the j-loop for h_u could collapse to j=js=1,je+1. -RWH
439  do j=js-2,je+2 ; do i=isq-1,ieq+1
440  h_u(i,j) = 0.5 * (h(i,j,k) + h(i+1,j,k))
441  enddo ; enddo
442  do j=jsq-1,jeq+1 ; do i=is-2,ie+2
443  h_v(i,j) = 0.5 * (h(i,j,k) + h(i,j+1,k))
444  enddo ; enddo
445 
446  ! Adjust contributions to shearing strain and interpolated values of
447  ! thicknesses on open boundaries.
448  if (apply_obc) then ; do n=1,obc%number_of_segments
449  j = obc%segment(n)%HI%JsdB ; i = obc%segment(n)%HI%IsdB
450  if (obc%zero_strain .or. obc%freeslip_strain) then
451  if (obc%segment(n)%is_N_or_S .and. (j >= js-2) .and. (j <= jeq+1)) then
452  do i=obc%segment(n)%HI%IsdB,obc%segment(n)%HI%IedB
453  if (obc%zero_strain) then
454  dvdx(i,j) = 0. ; dudy(i,j) = 0.
455  elseif (obc%freeslip_strain) then
456  dudy(i,j) = 0.
457  endif
458  enddo
459  elseif (obc%segment(n)%is_E_or_W .and. (i >= is-2) .and. (i <= ieq+1)) then
460  do j=obc%segment(n)%HI%JsdB,obc%segment(n)%HI%JedB
461  if (obc%zero_strain) then
462  dvdx(i,j) = 0. ; dudy(i,j) = 0.
463  elseif (obc%freeslip_strain) then
464  dvdx(i,j) = 0.
465  endif
466  enddo
467  endif
468  endif
469  if (obc%segment(n)%direction == obc_direction_n) then
470  ! There are extra wide halos here to accomodate the cross-corner-point
471  ! OBC projections, but they might not be necessary if the accelerations
472  ! are always zeroed out at OBC points, in which case the i-loop below
473  ! becomes do i=is-1,ie+1. -RWH
474  if ((j >= jsq-1) .and. (j <= jeq+1)) then
475  do i = max(is-2,obc%segment(n)%HI%isd), min(ie+2,obc%segment(n)%HI%ied)
476  h_v(i,j) = h(i,j,k)
477  enddo
478  endif
479  elseif (obc%segment(n)%direction == obc_direction_s) then
480  if ((j >= jsq-1) .and. (j <= jeq+1)) then
481  do i = max(is-2,obc%segment(n)%HI%isd), min(ie+2,obc%segment(n)%HI%ied)
482  h_v(i,j) = h(i,j+1,k)
483  enddo
484  endif
485  elseif (obc%segment(n)%direction == obc_direction_e) then
486  if ((i >= isq-1) .and. (i <= ieq+1)) then
487  do j = max(js-2,obc%segment(n)%HI%jsd), min(je+2,obc%segment(n)%HI%jed)
488  h_u(i,j) = h(i,j,k)
489  enddo
490  endif
491  elseif (obc%segment(n)%direction == obc_direction_w) then
492  if ((i >= isq-1) .and. (i <= ieq+1)) then
493  do j = max(js-2,obc%segment(n)%HI%jsd), min(je+2,obc%segment(n)%HI%jed)
494  h_u(i,j) = h(i+1,j,k)
495  enddo
496  endif
497  endif
498  enddo ; endif
499  ! Now project thicknesses across corner points on OBCs.
500  if (apply_obc) then ; do n=1,obc%number_of_segments
501  j = obc%segment(n)%HI%JsdB ; i = obc%segment(n)%HI%IsdB
502  if (obc%segment(n)%direction == obc_direction_n) then
503  if ((j >= js-2) .and. (j <= je)) then
504  do i = max(isq-1,obc%segment(n)%HI%IsdB), min(ieq+1,obc%segment(n)%HI%IedB)
505  h_u(i,j+1) = h_u(i,j)
506  enddo
507  endif
508  elseif (obc%segment(n)%direction == obc_direction_s) then
509  if ((j >= js-1) .and. (j <= je+1)) then
510  do i = max(isq-1,obc%segment(n)%HI%isd), min(ieq+1,obc%segment(n)%HI%ied)
511  h_u(i,j) = h_u(i,j+1)
512  enddo
513  endif
514  elseif (obc%segment(n)%direction == obc_direction_e) then
515  if ((i >= is-2) .and. (i <= ie)) then
516  do j = max(jsq-1,obc%segment(n)%HI%jsd), min(jeq+1,obc%segment(n)%HI%jed)
517  h_v(i+1,j) = h_v(i,j)
518  enddo
519  endif
520  elseif (obc%segment(n)%direction == obc_direction_w) then
521  if ((i >= is-1) .and. (i <= ie+1)) then
522  do j = max(jsq-1,obc%segment(n)%HI%jsd), min(jeq+1,obc%segment(n)%HI%jed)
523  h_v(i,j) = h_v(i+1,j)
524  enddo
525  endif
526  endif
527  enddo ; endif
528 
529  if (cs%no_slip) then
530  do j=js-2,jeq+1 ; do i=is-2,ieq+1
531  sh_xy(i,j) = (2.0-g%mask2dBu(i,j)) * ( dvdx(i,j) + dudy(i,j) )
532  enddo ; enddo
533  else
534  do j=js-2,jeq+1 ; do i=is-2,ieq+1
535  sh_xy(i,j) = g%mask2dBu(i,j) * ( dvdx(i,j) + dudy(i,j) )
536  enddo ; enddo
537  endif
538 
539  if (cs%no_slip) then
540  do j=js-2,jeq+1 ; do i=is-2,ieq+1
541  vort_xy(i,j) = (2.0-g%mask2dBu(i,j)) * ( dvdx(i,j) - dudy(i,j) )
542  enddo ; enddo
543  else
544  do j=js-2,jeq+1 ; do i=is-2,ieq+1
545  vort_xy(i,j) = g%mask2dBu(i,j) * ( dvdx(i,j) - dudy(i,j) )
546  enddo ; enddo
547  endif
548 
549 ! Vorticity gradient
550  do j=js-2,jeq+1 ; do i=is-1,ieq+1
551  vort_xy_dx(i,j) = cs%DY_dxBu(i,j)*(vort_xy(i,j)*g%IdyCu(i,j) - vort_xy(i-1,j)*g%IdyCu(i-1,j))
552  enddo ; enddo
553 
554  do j=js-1,jeq+1 ; do i=is-2,ieq+1
555  vort_xy_dy(i,j) = cs%DX_dyBu(i,j)*(vort_xy(i,j)*g%IdxCv(i,j) - vort_xy(i,j-1)*g%IdxCv(i,j-1))
556  enddo ; enddo
557 
558 ! Divergence gradient
559  do j=jsq-1,jeq+2 ; do i=is-2,ieq+1
560  div_xx_dx(i,j) = g%IdxCu(i,j)*(div_xx(i+1,j) - div_xx(i,j))
561  enddo ; enddo
562 
563  do j=js-2,jeq+1 ; do i=isq-1,ieq+2
564  div_xx_dy(i,j) = g%IdyCv(i,j)*(div_xx(i,j+1) - div_xx(i,j))
565  enddo ; enddo
566 
567 ! Coefficient for modified Leith
568  if (cs%Modified_Leith) then
569  mod_leith = 1.0
570  else
571  mod_leith = 0.0
572  endif
573 
574 ! Evaluate u0 = x.Div(Grad u) and v0 = y.Div( Grad u)
575  if (cs%biharmonic) then
576  do j=js-1,jeq+1 ; do i=isq-1,ieq+1
577  u0(i,j) = cs%IDXDY2u(i,j)*(cs%DY2h(i+1,j)*sh_xx(i+1,j) - cs%DY2h(i,j)*sh_xx(i,j)) + &
578  cs%IDX2dyCu(i,j)*(cs%DX2q(i,j)*sh_xy(i,j) - cs%DX2q(i,j-1)*sh_xy(i,j-1))
579  enddo ; enddo
580  do j=jsq-1,jeq+1 ; do i=is-1,ieq+1
581  v0(i,j) = cs%IDXDY2v(i,j)*(cs%DY2q(i,j)*sh_xy(i,j) - cs%DY2q(i-1,j)*sh_xy(i-1,j)) - &
582  cs%IDX2dyCv(i,j)*(cs%DX2h(i,j+1)*sh_xx(i,j+1) - cs%DX2h(i,j)*sh_xx(i,j))
583  enddo ; enddo
584  if (apply_obc .and. obc%zero_biharmonic) then
585  do n=1,obc%number_of_segments
586  i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
587  if (obc%segment(n)%is_N_or_S .and. (j >= jsq-1) .and. (j <= jeq+1)) then
588  do i=obc%segment(n)%HI%isd,obc%segment(n)%HI%ied
589  v0(i,j) = 0.
590  enddo
591  elseif (obc%segment(n)%is_E_or_W .and. (i >= isq-1) .and. (i <= ieq+1)) then
592  do j=obc%segment(n)%HI%jsd,obc%segment(n)%HI%jed
593  u0(i,j) = 0.
594  enddo
595  endif
596  enddo
597  endif
598  endif
599 
600  do j=jsq,jeq+1 ; do i=isq,ieq+1
601  if ((cs%Smagorinsky_Kh) .or. (cs%Smagorinsky_Ah)) then
602  shear_mag = sqrt(sh_xx(i,j)*sh_xx(i,j) + &
603  0.25*((sh_xy(i-1,j-1)*sh_xy(i-1,j-1) + sh_xy(i,j)*sh_xy(i,j)) + &
604  (sh_xy(i-1,j)*sh_xy(i-1,j) + sh_xy(i,j-1)*sh_xy(i,j-1))))
605  endif
606  if ((cs%Leith_Kh) .or. (cs%Leith_Ah)) then
607  vort_mag = sqrt( &
608  0.5*((vort_xy_dx(i,j-1)*vort_xy_dx(i,j-1) + vort_xy_dx(i,j)*vort_xy_dx(i,j)) + &
609  (vort_xy_dy(i-1,j)*vort_xy_dy(i-1,j) + vort_xy_dy(i,j)*vort_xy_dy(i,j))) + &
610  mod_leith*0.5*((div_xx_dx(i,j)*div_xx_dx(i,j) + div_xx_dx(i-1,j)*div_xx_dx(i-1,j)) + &
611  (div_xx_dy(i,j)*div_xx_dy(i,j) + div_xx_dy(i,j-1)*div_xx_dy(i,j-1))))
612  endif
613  if (cs%better_bound_Ah .or. cs%better_bound_Kh) then
614  hrat_min = min(1.0, min(h_u(i,j), h_u(i-1,j), h_v(i,j), h_v(i,j-1)) / &
615  (h(i,j,k) + h_neglect) )
616  visc_bound_rem = 1.0
617  endif
618 
619  if (cs%Laplacian) then
620  ! Determine the Laplacian viscosity at h points, using the
621  ! largest value from several parameterizations.
622  kh_scale = 1.0 ; if (rescale_kh) kh_scale = varmix%Res_fn_h(i,j)
623  khsm = 0.0; khlth = 0.0
624  if ((cs%Smagorinsky_Kh) .or. (cs%Leith_Kh)) then
625  if (cs%Smagorinsky_Kh) &
626  khsm = cs%LAPLAC_CONST_xx(i,j) * shear_mag
627  if (cs%Leith_Kh) &
628  khlth = cs%LAPLAC3_CONST_xx(i,j) * vort_mag
629  kh = kh_scale * max(khlth, max(cs%Kh_bg_xx(i,j), khsm))
630  if (cs%bound_Kh .and. .not.cs%better_bound_Kh) &
631  kh = min(kh, cs%Kh_Max_xx(i,j))
632  else
633  kh = kh_scale * cs%Kh_bg_xx(i,j)
634  endif
635 
636  if (use_meke_ku) then
637  kh = kh + meke%Ku(i,j)
638  endif
639  if (cs%better_bound_Kh) then
640  if (kh >= hrat_min*cs%Kh_Max_xx(i,j)) then
641  visc_bound_rem = 0.0
642  kh = hrat_min*cs%Kh_Max_xx(i,j)
643  else
644  !visc_bound_rem = 1.0 - abs(Kh) / (hrat_min*CS%Kh_Max_xx(i,j))
645  visc_bound_rem = 1.0 - kh / (hrat_min*cs%Kh_Max_xx(i,j))
646  endif
647  endif
648 
649  if (cs%id_Kh_h>0) kh_h(i,j,k) = kh
650 
651  str_xx(i,j) = -kh * sh_xx(i,j)
652  else ! not Laplacian
653  str_xx(i,j) = 0.0
654  endif ! Laplacian
655 
656  if (cs%biharmonic) then
657 ! Determine the biharmonic viscosity at h points, using the
658 ! largest value from several parameterizations.
659  ahsm = 0.0; ahlth = 0.0
660  if ((cs%Smagorinsky_Ah) .or. (cs%Leith_Ah)) then
661  if (cs%Smagorinsky_Ah) then
662  if (cs%bound_Coriolis) then
663  ahsm = shear_mag * (cs%BIHARM_CONST_xx(i,j) + &
664  cs%Biharm_Const2_xx(i,j)*shear_mag)
665  else
666  ahsm = cs%BIHARM_CONST_xx(i,j) * shear_mag
667  endif
668  endif
669  if (cs%Leith_Ah) &
670  ahlth = vort_mag * (cs%BIHARM_CONST_xx(i,j))
671  ah = max(max(cs%Ah_bg_xx(i,j), ahsm),ahlth)
672  if (cs%bound_Ah .and. .not.cs%better_bound_Ah) &
673  ah = min(ah, cs%Ah_Max_xx(i,j))
674  else
675  ah = cs%Ah_bg_xx(i,j)
676  endif ! Smagorinsky_Ah or Leith_Ah
677 
678  if (cs%better_bound_Ah) then
679  ah = min(ah, visc_bound_rem*hrat_min*cs%Ah_Max_xx(i,j))
680  endif
681 
682  if (cs%id_Ah_h>0) ah_h(i,j,k) = ah
683 
684  str_xx(i,j) = str_xx(i,j) + ah * &
685  (cs%DY_dxT(i,j)*(g%IdyCu(i,j)*u0(i,j) - g%IdyCu(i-1,j)*u0(i-1,j)) - &
686  cs%DX_dyT(i,j) *(g%IdxCv(i,j)*v0(i,j) - g%IdxCv(i,j-1)*v0(i,j-1)))
687 
688  ! Keep a copy of the biharmonic contribution for backscatter parameterization
689  bhstr_xx(i,j) = ah * &
690  (cs%DY_dxT(i,j)*(g%IdyCu(i,j)*u0(i,j) - g%IdyCu(i-1,j)*u0(i-1,j)) - &
691  cs%DX_dyT(i,j) *(g%IdxCv(i,j)*v0(i,j) - g%IdxCv(i,j-1)*v0(i,j-1)))
692  bhstr_xx(i,j) = bhstr_xx(i,j) * (h(i,j,k) * cs%reduction_xx(i,j))
693 
694  endif ! biharmonic
695 
696  str_xx(i,j) = str_xx(i,j) * (h(i,j,k) * cs%reduction_xx(i,j))
697  enddo ; enddo
698 
699  if (cs%biharmonic) then
700  ! Gradient of Laplacian, for use in bi-harmonic term
701  do j=js-1,jeq ; do i=is-1,ieq
702  dvdx(i,j) = cs%DY_dxBu(i,j)*(v0(i+1,j)*g%IdyCv(i+1,j) - v0(i,j)*g%IdyCv(i,j))
703  dudy(i,j) = cs%DX_dyBu(i,j)*(u0(i,j+1)*g%IdxCu(i,j+1) - u0(i,j)*g%IdxCu(i,j))
704  enddo ; enddo
705  ! Adjust contributions to shearing strain on open boundaries.
706  if (apply_obc) then ; if (obc%zero_strain .or. obc%freeslip_strain) then
707  do n=1,obc%number_of_segments
708  j = obc%segment(n)%HI%JsdB ; i = obc%segment(n)%HI%IsdB
709  if (obc%segment(n)%is_N_or_S .and. (j >= js-1) .and. (j <= jeq)) then
710  do i=obc%segment(n)%HI%IsdB,obc%segment(n)%HI%IedB
711  if (obc%zero_strain) then
712  dvdx(i,j) = 0. ; dudy(i,j) = 0.
713  elseif (obc%freeslip_strain) then
714  dudy(i,j) = 0.
715  endif
716  enddo
717  elseif (obc%segment(n)%is_E_or_W .and. (i >= is-1) .and. (i <= ieq)) then
718  do j=obc%segment(n)%HI%JsdB,obc%segment(n)%HI%JedB
719  if (obc%zero_strain) then
720  dvdx(i,j) = 0. ; dudy(i,j) = 0.
721  elseif (obc%freeslip_strain) then
722  dvdx(i,j) = 0.
723  endif
724  enddo
725  endif
726  enddo
727  endif ; endif
728  endif
729 
730  do j=js-1,jeq ; do i=is-1,ieq
731  if ((cs%Smagorinsky_Kh) .or. (cs%Smagorinsky_Ah)) then
732  shear_mag = sqrt(sh_xy(i,j)*sh_xy(i,j) + &
733  0.25*((sh_xx(i,j)*sh_xx(i,j) + sh_xx(i+1,j+1)*sh_xx(i+1,j+1)) + &
734  (sh_xx(i,j+1)*sh_xx(i,j+1) + sh_xx(i+1,j)*sh_xx(i+1,j))))
735  endif
736  if ((cs%Leith_Kh) .or. (cs%Leith_Ah)) &
737  vort_mag = sqrt( &
738  0.5*((vort_xy_dx(i,j)*vort_xy_dx(i,j) + vort_xy_dx(i+1,j)*vort_xy_dx(i+1,j)) + &
739  (vort_xy_dy(i,j)*vort_xy_dy(i,j) + vort_xy_dy(i,j+1)*vort_xy_dy(i,j+1))) + &
740  mod_leith*0.5*((div_xx_dx(i,j)*div_xx_dx(i,j) + div_xx_dx(i,j+1)*div_xx_dx(i,j+1)) + &
741  (div_xx_dy(i,j)*div_xx_dy(i,j) + div_xx_dy(i+1,j)*div_xx_dy(i+1,j))))
742 
743  h2uq = 4.0 * h_u(i,j) * h_u(i,j+1)
744  h2vq = 4.0 * h_v(i,j) * h_v(i+1,j)
745 ! hq = 2.0 * h2uq * h2vq / (h_neglect3 + (h2uq + h2vq) * &
746 ! ((h(i,j,k) + h(i+1,j+1,k)) + (h(i,j+1,k) + h(i+1,j,k))))
747  hq = 2.0 * h2uq * h2vq / (h_neglect3 + (h2uq + h2vq) * &
748  ((h_u(i,j) + h_u(i,j+1)) + (h_v(i,j) + h_v(i+1,j))))
749 
750  if (cs%better_bound_Ah .or. cs%better_bound_Kh) then
751  hrat_min = min(1.0, min(h_u(i,j), h_u(i,j+1), h_v(i,j), h_v(i+1,j)) / &
752  (hq + h_neglect) )
753  visc_bound_rem = 1.0
754  endif
755 
756  if (cs%no_slip .and. (g%mask2dBu(i,j) < 0.5)) then
757  if ((g%mask2dCu(i,j) + g%mask2dCu(i,j+1)) + &
758  (g%mask2dCv(i,j) + g%mask2dCv(i+1,j)) > 0.0) then
759  ! This is a coastal vorticity point, so modify hq and hrat_min.
760 
761  hu = g%mask2dCu(i,j) * h_u(i,j) + g%mask2dCu(i,j+1) * h_u(i,j+1)
762  hv = g%mask2dCv(i,j) * h_v(i,j) + g%mask2dCv(i+1,j) * h_v(i+1,j)
763  if ((g%mask2dCu(i,j) + g%mask2dCu(i,j+1)) * &
764  (g%mask2dCv(i,j) + g%mask2dCv(i+1,j)) == 0.0) then
765  ! Only one of hu and hv is nonzero, so just add them.
766  hq = hu + hv
767  hrat_min = 1.0
768  else
769  ! Both hu and hv are nonzero, so take the harmonic mean.
770  hq = 2.0 * (hu * hv) / ((hu + hv) + h_neglect)
771  hrat_min = min(1.0, min(hu, hv) / (hq + h_neglect) )
772  endif
773  endif
774  endif
775 
776  if (cs%Laplacian) then
777  ! Determine the Laplacian viscosity at q points, using the
778  ! largest value from several parameterizations.
779  kh_scale = 1.0 ; if (rescale_kh) kh_scale = varmix%Res_fn_q(i,j)
780  khsm = 0.0; khlth = 0.0
781  if ((cs%Smagorinsky_Kh) .or. (cs%Leith_Kh)) then
782  if (cs%Smagorinsky_Kh) &
783  khsm = cs%LAPLAC_CONST_xy(i,j) * shear_mag
784  if (cs%Leith_Kh) &
785  khlth = cs%LAPLAC3_CONST_xy(i,j) * vort_mag
786  kh = kh_scale * max(max(cs%Kh_bg_xy(i,j), khsm), khlth)
787  if (cs%bound_Kh .and. .not.cs%better_bound_Kh) &
788  kh = min(kh, cs%Kh_Max_xy(i,j))
789  else
790  kh = kh_scale * cs%Kh_bg_xy(i,j)
791  endif
792  if (use_meke_ku) then
793  kh = kh + 0.25*( (meke%Ku(i,j)+meke%Ku(i+1,j+1)) &
794  +(meke%Ku(i+1,j)+meke%Ku(i,j+1)) )
795  endif
796  ! Place a floor on the viscosity, if desired.
797  kh = max(kh,cs%Kh_bg_min)
798 
799  if (cs%better_bound_Kh) then
800  if (kh >= hrat_min*cs%Kh_Max_xy(i,j)) then
801  visc_bound_rem = 0.0
802  kh = hrat_min*cs%Kh_Max_xy(i,j)
803  elseif (cs%Kh_Max_xy(i,j)>0.) then
804  !visc_bound_rem = 1.0 - abs(Kh) / (hrat_min*CS%Kh_Max_xy(I,J))
805  visc_bound_rem = 1.0 - kh / (hrat_min*cs%Kh_Max_xy(i,j))
806  endif
807  endif
808 
809  if (cs%id_Kh_q>0) kh_q(i,j,k) = kh
810 
811  str_xy(i,j) = -kh * sh_xy(i,j)
812  else ! not Laplacian
813  str_xy(i,j) = 0.0
814  endif ! Laplacian
815 
816  if (cs%biharmonic) then
817  ! Determine the biharmonic viscosity at q points, using the
818  ! largest value from several parameterizations.
819  ahsm = 0.0; ahlth = 0.0
820  if (cs%Smagorinsky_Ah .or. cs%Leith_Ah) then
821  if (cs%Smagorinsky_Ah) then
822  if (cs%bound_Coriolis) then
823  ahsm = shear_mag * (cs%BIHARM_CONST_xy(i,j) + &
824  cs%Biharm_Const2_xy(i,j)*shear_mag)
825  else
826  ahsm = cs%BIHARM_CONST_xy(i,j) * shear_mag
827  endif
828  endif
829  if (cs%Leith_Ah) &
830  ahlth = vort_mag * (cs%BIHARM5_CONST_xy(i,j))
831  ah = max(max(cs%Ah_bg_xy(i,j), ahsm),ahlth)
832  if (cs%bound_Ah .and. .not.cs%better_bound_Ah) &
833  ah = min(ah, cs%Ah_Max_xy(i,j))
834  else
835  ah = cs%Ah_bg_xy(i,j)
836  endif ! Smagorinsky_Ah or Leith_Ah
837  if (cs%better_bound_Ah) then
838  ah = min(ah, visc_bound_rem*hrat_min*cs%Ah_Max_xy(i,j))
839  endif
840 
841  if (cs%id_Ah_q>0) ah_q(i,j,k) = ah
842 
843  str_xy(i,j) = str_xy(i,j) + ah * ( dvdx(i,j) + dudy(i,j) )
844 
845  ! Keep a copy of the biharmonic contribution for backscatter parameterization
846  bhstr_xy(i,j) = ah * ( dvdx(i,j) + dudy(i,j) ) * &
847  (hq * g%mask2dBu(i,j) * cs%reduction_xy(i,j))
848 
849  endif ! biharmonic
850 
851  if (cs%no_slip) then
852  str_xy(i,j) = str_xy(i,j) * (hq * cs%reduction_xy(i,j))
853  else
854  str_xy(i,j) = str_xy(i,j) * (hq * g%mask2dBu(i,j) * cs%reduction_xy(i,j))
855  endif
856  enddo ; enddo
857 
858  do j=js,je ; do i=isq,ieq
859 ! Evaluate 1/h x.Div(h Grad u) or the biharmonic equivalent.
860  diffu(i,j,k) = ((g%IdyCu(i,j)*(cs%DY2h(i,j) *str_xx(i,j) - &
861  cs%DY2h(i+1,j)*str_xx(i+1,j)) + &
862  g%IdxCu(i,j)*(cs%DX2q(i,j-1)*str_xy(i,j-1) - &
863  cs%DX2q(i,j) *str_xy(i,j))) * &
864  g%IareaCu(i,j)) / (h_u(i,j) + h_neglect)
865 
866  enddo ; enddo
867  if (apply_obc) then
868  ! This is not the right boundary condition. If all the masking of tendencies are done
869  ! correctly later then eliminating this block should not change answers.
870  do n=1,obc%number_of_segments
871  if (obc%segment(n)%is_E_or_W) then
872  i = obc%segment(n)%HI%IsdB
873  do j=obc%segment(n)%HI%jsd,obc%segment(n)%HI%jed
874  diffu(i,j,k) = 0.
875  enddo
876  endif
877  enddo
878  endif
879 
880 ! Evaluate 1/h y.Div(h Grad u) or the biharmonic equivalent.
881  do j=jsq,jeq ; do i=is,ie
882  diffv(i,j,k) = ((g%IdyCv(i,j)*(cs%DY2q(i-1,j)*str_xy(i-1,j) - &
883  cs%DY2q(i,j) *str_xy(i,j)) - &
884  g%IdxCv(i,j)*(cs%DX2h(i,j) *str_xx(i,j) - &
885  cs%DX2h(i,j+1)*str_xx(i,j+1))) * &
886  g%IareaCv(i,j)) / (h_v(i,j) + h_neglect)
887  enddo ; enddo
888  if (apply_obc) then
889  ! This is not the right boundary condition. If all the masking of tendencies are done
890  ! correctly later then eliminating this block should not change answers.
891  do n=1,obc%number_of_segments
892  if (obc%segment(n)%is_N_or_S) then
893  j = obc%segment(n)%HI%JsdB
894  do i=obc%segment(n)%HI%isd,obc%segment(n)%HI%ied
895  diffv(i,j,k) = 0.
896  enddo
897  endif
898  enddo
899  endif
900 
901  if (find_frictwork) then ; do j=js,je ; do i=is,ie
902  ! Diagnose str_xx*d_x u - str_yy*d_y v + str_xy*(d_y u + d_x v)
903  frictwork(i,j,k) = gv%H_to_kg_m2 * ( &
904  (str_xx(i,j)*(u(i,j,k)-u(i-1,j,k))*g%IdxT(i,j) &
905  -str_xx(i,j)*(v(i,j,k)-v(i,j-1,k))*g%IdyT(i,j)) &
906  +0.25*((str_xy(i,j)*( &
907  (u(i,j+1,k)-u(i,j,k))*g%IdyBu(i,j) &
908  +(v(i+1,j,k)-v(i,j,k))*g%IdxBu(i,j) ) &
909  +str_xy(i-1,j-1)*( &
910  (u(i-1,j,k)-u(i-1,j-1,k))*g%IdyBu(i-1,j-1) &
911  +(v(i,j-1,k)-v(i-1,j-1,k))*g%IdxBu(i-1,j-1) )) &
912  +(str_xy(i-1,j)*( &
913  (u(i-1,j+1,k)-u(i-1,j,k))*g%IdyBu(i-1,j) &
914  +(v(i,j,k)-v(i-1,j,k))*g%IdxBu(i-1,j) ) &
915  +str_xy(i,j-1)*( &
916  (u(i,j,k)-u(i,j-1,k))*g%IdyBu(i,j-1) &
917  +(v(i+1,j-1,k)-v(i,j-1,k))*g%IdxBu(i,j-1) )) ) )
918  enddo ; enddo ; endif
919 
920  ! Make a similar calculation as for FrictWork above but accumulating into
921  ! the vertically integrated MEKE source term, and adjusting for any
922  ! energy loss seen as a reduction in the [biharmonic] frictional source term.
923  if (find_frictwork .and. associated(meke)) then ; if (associated(meke%mom_src)) then
924  if (k==1) then
925  do j=js,je ; do i=is,ie
926  meke%mom_src(i,j) = 0.
927  enddo ; enddo
928  endif
929  if (meke%backscatter_Ro_c /= 0.) then
930  do j=js,je ; do i=is,ie
931  fath = 0.25*( (abs(g%CoriolisBu(i-1,j-1)) + abs(g%CoriolisBu(i,j))) &
932  +(abs(g%CoriolisBu(i-1,j)) + abs(g%CoriolisBu(i,j-1))) )
933  shear_mag = sqrt(sh_xx(i,j)*sh_xx(i,j) + &
934  0.25*((sh_xy(i-1,j-1)*sh_xy(i-1,j-1) + sh_xy(i,j)*sh_xy(i,j)) + &
935  (sh_xy(i-1,j)*sh_xy(i-1,j) + sh_xy(i,j-1)*sh_xy(i,j-1))))
936  fath = fath ** meke%backscatter_Ro_pow ! f^n
937  shear_mag = ( ( shear_mag ** meke%backscatter_Ro_pow ) + 1.e-30 ) &
938  * meke%backscatter_Ro_c ! c * D^n
939  ! The Rossby number function is g(Ro) = 1/(1+c.Ro^n)
940  ! RoScl = 1 - g(Ro)
941  roscl = shear_mag / ( fath + shear_mag ) ! = 1 - f^n/(f^n+c*D^n)
942  meke%mom_src(i,j) = meke%mom_src(i,j) + gv%H_to_kg_m2 * ( &
943  ((str_xx(i,j)-roscl*bhstr_xx(i,j))*(u(i,j,k)-u(i-1,j,k))*g%IdxT(i,j) &
944  -(str_xx(i,j)-roscl*bhstr_xx(i,j))*(v(i,j,k)-v(i,j-1,k))*g%IdyT(i,j)) &
945  +0.25*(((str_xy(i,j)-roscl*bhstr_xy(i,j))*( &
946  (u(i,j+1,k)-u(i,j,k))*g%IdyBu(i,j) &
947  +(v(i+1,j,k)-v(i,j,k))*g%IdxBu(i,j) ) &
948  +(str_xy(i-1,j-1)-roscl*bhstr_xy(i-1,j-1))*( &
949  (u(i-1,j,k)-u(i-1,j-1,k))*g%IdyBu(i-1,j-1) &
950  +(v(i,j-1,k)-v(i-1,j-1,k))*g%IdxBu(i-1,j-1) )) &
951  +((str_xy(i-1,j)-roscl*bhstr_xy(i-1,j))*( &
952  (u(i-1,j+1,k)-u(i-1,j,k))*g%IdyBu(i-1,j) &
953  +(v(i,j,k)-v(i-1,j,k))*g%IdxBu(i-1,j) ) &
954  +(str_xy(i,j-1)-roscl*bhstr_xy(i,j-1))*( &
955  (u(i,j,k)-u(i,j-1,k))*g%IdyBu(i,j-1) &
956  +(v(i+1,j-1,k)-v(i,j-1,k))*g%IdxBu(i,j-1) )) ) )
957  enddo ; enddo
958  else
959  do j=js,je ; do i=is,ie
960  meke%mom_src(i,j) = meke%mom_src(i,j) + frictwork(i,j,k)
961  enddo ; enddo
962  endif
963  endif ; endif
964 
965  enddo ! end of k loop
966 
967 ! Offer fields for diagnostic averaging.
968  if (cs%id_diffu>0) call post_data(cs%id_diffu, diffu, cs%diag)
969  if (cs%id_diffv>0) call post_data(cs%id_diffv, diffv, cs%diag)
970  if (cs%id_FrictWork>0) call post_data(cs%id_FrictWork, frictwork, cs%diag)
971  if (cs%id_Ah_h>0) call post_data(cs%id_Ah_h, ah_h, cs%diag)
972  if (cs%id_Ah_q>0) call post_data(cs%id_Ah_q, ah_q, cs%diag)
973  if (cs%id_Kh_h>0) call post_data(cs%id_Kh_h, kh_h, cs%diag)
974  if (cs%id_Kh_q>0) call post_data(cs%id_Kh_q, kh_q, cs%diag)
975 
976  if (cs%id_FrictWorkIntz > 0) then
977  do j=js,je
978  do i=is,ie ; frictworkintz(i,j) = frictwork(i,j,1) ; enddo
979  do k=2,nz ; do i=is,ie
980  frictworkintz(i,j) = frictworkintz(i,j) + frictwork(i,j,k)
981  enddo ; enddo
982  enddo
983  call post_data(cs%id_FrictWorkIntz, frictworkintz, cs%diag)
984  endif
985 
986 
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:19
Here is the call graph for this function:
Here is the caller graph for this function: