diff --git a/config_src/mct_driver/MOM_surface_forcing.F90 b/config_src/mct_driver/MOM_surface_forcing.F90
index 5b17cbbcff..47e676a3d3 100644
--- a/config_src/mct_driver/MOM_surface_forcing.F90
+++ b/config_src/mct_driver/MOM_surface_forcing.F90
@@ -215,7 +215,6 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, Time, G, US, CS, &
type(forcing), intent(inout) :: fluxes !< A structure containing pointers to
!! all possible mass, heat or salt flux forcing fields.
!! Unused fields have NULL ptrs.
-
type(time_type), intent(in) :: Time !< The time of the fluxes, used for interpolating the
!! salinity to the right time, when it is being restored.
type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
@@ -244,7 +243,6 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, Time, G, US, CS, &
integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, i0, j0
integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, isr, ier, jsr, jer
- integer :: isc_bnd, iec_bnd, jsc_bnd, jec_bnd
logical :: restore_salinity ! local copy of the argument restore_salt, if it
! is present, or false (no restoring) otherwise.
@@ -395,9 +393,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, Time, G, US, CS, &
enddo; enddo
endif
- !i0 = is - isc_bnd ; j0 = js - jsc_bnd ???
- i0 = 0; j0 = 0 ! TODO: is this right?
-
+ ! obtain fluxes from IOB
+ i0 = 0; j0 = 0
do j=js,je ; do i=is,ie
! liquid precipitation (rain)
if (associated(fluxes%lprec)) &
@@ -456,8 +453,27 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, Time, G, US, CS, &
fluxes%seaice_melt(i,j) = G%mask2dT(i,j) * IOB%seaice_melt(i-i0,j-j0)
! latent heat flux (W/m^2)
- if (associated(fluxes%latent)) &
- fluxes%latent(i,j) = G%mask2dT(i,j) * IOB%latent_flux(i-i0,j-j0)
+ ! old method, latent = IOB%q_flux(i-i0,j-j0)*CS%latent_heat_vapor
+ !if (associated(fluxes%latent)) &
+ ! fluxes%latent(i,j) = G%mask2dT(i,j) * IOB%latent_flux(i-i0,j-j0)
+ ! new method
+ fluxes%latent(i,j) = 0.0
+ ! contribution from frozen ppt
+ if (associated(fluxes%fprec)) then
+ fluxes%latent(i,j) = fluxes%latent(i,j) + IOB%fprec(i-i0,j-j0)*CS%latent_heat_fusion
+ fluxes%latent_fprec_diag(i,j) = G%mask2dT(i,j) * IOB%fprec(i-i0,j-j0)*CS%latent_heat_fusion
+ endif
+ ! contribution from frozen runoff
+ if (associated(fluxes%frunoff)) then
+ fluxes%latent(i,j) = fluxes%latent(i,j) + IOB%rofi_flux(i-i0,j-j0)*CS%latent_heat_fusion
+ fluxes%latent_frunoff_diag(i,j) = G%mask2dT(i,j) * IOB%rofi_flux(i-i0,j-j0)*CS%latent_heat_fusion
+ endif
+ ! contribution from evaporation
+ if (associated(IOB%q_flux)) then
+ fluxes%latent(i,j) = fluxes%latent(i,j) + IOB%q_flux(i-i0,j-j0)*CS%latent_heat_vapor
+ fluxes%latent_evap_diag(i,j) = G%mask2dT(i,j) * IOB%q_flux(i-i0,j-j0)*CS%latent_heat_vapor
+ endif
+ fluxes%latent(i,j) = G%mask2dT(i,j) * fluxes%latent(i,j)
if (associated(IOB%sw_flux_vis_dir)) &
fluxes%sw_vis_dir(i,j) = G%mask2dT(i,j) * IOB%sw_flux_vis_dir(i-i0,j-j0)
@@ -580,8 +596,9 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
!isc_bnd = index_bounds(1) ; iec_bnd = index_bounds(2)
!jsc_bnd = index_bounds(3) ; jec_bnd = index_bounds(4)
+ !if (is_root_pe()) write(*,*)'isc_bnd, jsc_bnd, iec_bnd, jec_bnd',isc_bnd, jsc_bnd, iec_bnd, jec_bnd
!i0 = is - isc_bnd ; j0 = js - jsc_bnd
- i0 = 0; j0 = 0 ! TODO: is this right?
+ i0 = 0; j0 = 0 ! TODO: is this right?
Irho0 = 1.0/CS%Rho0
diff --git a/config_src/mct_driver/ocn_cap_methods.F90 b/config_src/mct_driver/ocn_cap_methods.F90
index a8ac5310c7..7723f51a6c 100644
--- a/config_src/mct_driver/ocn_cap_methods.F90
+++ b/config_src/mct_driver/ocn_cap_methods.F90
@@ -50,11 +50,11 @@ subroutine ocn_import(x2o, ind, grid, ice_ocean_boundary, ocean_public, logunit,
! rotate taux and tauy from true zonal/meridional to local coordinates
! taux
ice_ocean_boundary%u_flux(i,j) = GRID%cos_rot(i,j) * x2o(ind%x2o_Foxx_taux,k) &
- + GRID%sin_rot(i,j) * x2o(ind%x2o_Foxx_tauy,k)
+ - GRID%sin_rot(i,j) * x2o(ind%x2o_Foxx_tauy,k)
! tauy
ice_ocean_boundary%v_flux(i,j) = GRID%cos_rot(i,j) * x2o(ind%x2o_Foxx_tauy,k) &
- - GRID%sin_rot(i,j) * x2o(ind%x2o_Foxx_taux,k)
+ + GRID%sin_rot(i,j) * x2o(ind%x2o_Foxx_taux,k)
! liquid precipitation (rain)
ice_ocean_boundary%lprec(i,j) = x2o(ind%x2o_Faxa_rain,k)
@@ -188,9 +188,9 @@ subroutine ocn_export(ind, ocn_public, grid, o2x, dt_int, ncouple_per_day)
o2x(ind%o2x_So_t, n) = ocn_public%t_surf(ig,jg) * grid%mask2dT(i,j)
o2x(ind%o2x_So_s, n) = ocn_public%s_surf(ig,jg) * grid%mask2dT(i,j)
! rotate ocn current from local tripolar grid to true zonal/meridional (inverse transformation)
- o2x(ind%o2x_So_u, n) = (grid%cos_rot(i,j) * ocn_public%u_surf(ig,jg) - &
+ o2x(ind%o2x_So_u, n) = (grid%cos_rot(i,j) * ocn_public%u_surf(ig,jg) + &
grid%sin_rot(i,j) * ocn_public%v_surf(ig,jg)) * grid%mask2dT(i,j)
- o2x(ind%o2x_So_v, n) = (grid%cos_rot(i,j) * ocn_public%v_surf(ig,jg) + &
+ o2x(ind%o2x_So_v, n) = (grid%cos_rot(i,j) * ocn_public%v_surf(ig,jg) - &
grid%sin_rot(i,j) * ocn_public%u_surf(ig,jg)) * grid%mask2dT(i,j)
! boundary layer depth (m)
@@ -270,8 +270,8 @@ subroutine ocn_export(ind, ocn_public, grid, o2x, dt_int, ncouple_per_day)
n = 0
do j=grid%jsc, grid%jec ; do i=grid%isc,grid%iec
n = n+1
- o2x(ind%o2x_So_dhdx, n) = grid%cos_rot(i,j) * sshx(i,j) - grid%sin_rot(i,j) * sshy(i,j)
- o2x(ind%o2x_So_dhdy, n) = grid%cos_rot(i,j) * sshy(i,j) + grid%sin_rot(i,j) * sshx(i,j)
+ o2x(ind%o2x_So_dhdx, n) = grid%cos_rot(i,j) * sshx(i,j) + grid%sin_rot(i,j) * sshy(i,j)
+ o2x(ind%o2x_So_dhdy, n) = grid%cos_rot(i,j) * sshy(i,j) - grid%sin_rot(i,j) * sshx(i,j)
enddo; enddo
end subroutine ocn_export
diff --git a/config_src/nuopc_driver/MOM_surface_forcing.F90 b/config_src/nuopc_driver/MOM_surface_forcing.F90
index 8348088b8a..69dda6b6d3 100644
--- a/config_src/nuopc_driver/MOM_surface_forcing.F90
+++ b/config_src/nuopc_driver/MOM_surface_forcing.F90
@@ -163,7 +163,7 @@ module MOM_surface_forcing
real, pointer, dimension(:,:) :: q_flux =>NULL() !< specific humidity flux [kg/m2/s]
real, pointer, dimension(:,:) :: salt_flux =>NULL() !< salt flux [kg/m2/s]
real, pointer, dimension(:,:) :: seaice_melt_heat =>NULL() !< sea ice and snow melt heat flux [W/m2]
- real, pointer, dimension(:,:) :: seaice_melt_water =>NULL() !< water flux due to sea ice and snow melting [kg/m2/s]
+ real, pointer, dimension(:,:) :: seaice_melt =>NULL() !< water flux due to sea ice and snow melting [kg/m2/s]
real, pointer, dimension(:,:) :: lw_flux =>NULL() !< long wave radiation [W/m2]
real, pointer, dimension(:,:) :: sw_flux_vis_dir =>NULL() !< direct visible sw radiation [W/m2]
real, pointer, dimension(:,:) :: sw_flux_vis_dif =>NULL() !< diffuse visible sw radiation [W/m2]
@@ -457,18 +457,18 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, US, CS, &
fluxes%heat_content_frunoff(i,j) = IOB%calving_hflx(i-i0,j-j0) * G%mask2dT(i,j)
if (associated(IOB%lw_flux)) &
- fluxes%LW(i,j) = IOB%lw_flux(i-i0,j-j0) * G%mask2dT(i,j)
+ fluxes%LW(i,j) = IOB%lw_flux(i-i0,j-j0) * G%mask2dT(i,j)
if (associated(IOB%t_flux)) &
- fluxes%sens(i,j) = IOB%t_flux(i-i0,j-j0) * G%mask2dT(i,j)
+ fluxes%sens(i,j) = IOB%t_flux(i-i0,j-j0) * G%mask2dT(i,j)
- ! ! sea ice and snow melt heat flux [W/m2]
- ! if (associated(fluxes%seaice_melt_heat)) &
- ! fluxes%seaice_melt_heat(i,j) = G%mask2dT(i,j) * IOB%seaice_melt_heat(i-i0,j-j0)
+ ! sea ice and snow melt heat flux [W/m2]
+ if (associated(IOB%seaice_melt_heat)) &
+ fluxes%seaice_melt_heat(i,j) = G%mask2dT(i,j) * IOB%seaice_melt_heat(i-i0,j-j0)
- ! ! water flux due to sea ice and snow melt [kg/m2/s]
- ! if (associated(fluxes%seaice_melt)) &
- ! fluxes%seaice_melt(i,j) = G%mask2dT(i,j) * IOB%seaice_melt_water(i-i0,j-j0)
+ ! water flux due to sea ice and snow melt [kg/m2/s]
+ if (associated(IOB%seaice_melt)) &
+ fluxes%seaice_melt(i,j) = G%mask2dT(i,j) * IOB%seaice_melt(i-i0,j-j0)
fluxes%latent(i,j) = 0.0
if (associated(IOB%fprec)) then
@@ -540,10 +540,9 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, US, CS, &
sign_for_net_FW_bug = 1.
if (CS%use_net_FW_adjustment_sign_bug) sign_for_net_FW_bug = -1.
do j=js,je ; do i=is,ie
- net_FW(i,j) = (((fluxes%lprec(i,j) + fluxes%fprec(i,j)) + &
+ net_FW(i,j) = (((fluxes%lprec(i,j) + fluxes%fprec(i,j) + fluxes%seaice_melt(i,j)) + &
(fluxes%lrunoff(i,j) + fluxes%frunoff(i,j))) + &
(fluxes%evap(i,j) + fluxes%vprec(i,j)) ) * G%areaT(i,j)
- ! net_FW(i,j) = netFW(i,j) + fluxes%seaice_melt(i,j) * G%areaT(i,j)
! The following contribution appears to be calculating the volume flux of sea-ice
! melt. This calculation is clearly WRONG if either sea-ice has variable
@@ -1087,7 +1086,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt,
call get_param(param_file, mdl, "USE_NET_FW_ADJUSTMENT_SIGN_BUG", &
CS%use_net_FW_adjustment_sign_bug, &
"If true, use the wrong sign for the adjustment to "//&
- "the net fresh-water.", default=.true.)
+ "the net fresh-water.", default=.false.)
call get_param(param_file, mdl, "ADJUST_NET_FRESH_WATER_BY_SCALING", &
CS%adjust_net_fresh_water_by_scaling, &
"If true, adjustments to net fresh water to achieve zero net are "//&
@@ -1368,8 +1367,8 @@ subroutine ice_ocn_bnd_type_chksum(id, timestep, iobt)
write(outunit,100) 'iobt%t_flux ' , mpp_chksum( iobt%t_flux )
write(outunit,100) 'iobt%q_flux ' , mpp_chksum( iobt%q_flux )
write(outunit,100) 'iobt%salt_flux ' , mpp_chksum( iobt%salt_flux )
- !write(outunit,100) 'iobt%seaice_melt_heat' , mpp_chksum( iobt%seaice_melt_heat)
- !write(outunit,100) 'iobt%seaice_melt_water' , mpp_chksum( iobt%seaice_melt_water)
+ write(outunit,100) 'iobt%seaice_melt_heat' , mpp_chksum( iobt%seaice_melt_heat)
+ write(outunit,100) 'iobt%seaice_melt ' , mpp_chksum( iobt%seaice_melt )
write(outunit,100) 'iobt%lw_flux ' , mpp_chksum( iobt%lw_flux )
write(outunit,100) 'iobt%sw_flux_vis_dir' , mpp_chksum( iobt%sw_flux_vis_dir)
write(outunit,100) 'iobt%sw_flux_vis_dif' , mpp_chksum( iobt%sw_flux_vis_dif)
diff --git a/config_src/nuopc_driver/mom_cap.F90 b/config_src/nuopc_driver/mom_cap.F90
index 24e60388b4..3992aae530 100644
--- a/config_src/nuopc_driver/mom_cap.F90
+++ b/config_src/nuopc_driver/mom_cap.F90
@@ -204,6 +204,8 @@
!! --------------------------|------------|-----------------|---------------------------------------|-------------------
!! inst_pres_height_surface | Pa | p | pressure of overlying sea ice and atmosphere
!! mass_of_overlying_sea_ice | kg | mi | mass of overlying sea ice | |
+!! seaice_melt_heat | W m-2 | seaice_melt_heat| sea ice and snow melt heat flux | |
+!! seaice_melt | kg m-2 s-1 | seaice_melt | water flux due to sea ice and snow melting | |
!! mean_calving_heat_flx | W m-2 | calving_hflx | heat flux, relative to 0C, of frozen land water into ocean
!! mean_calving_rate | kg m-2 s-1 | calving | mass flux of frozen runoff | |
!! mean_evap_rate | kg m-2 s-1 | q_flux | specific humidity flux |
@@ -377,6 +379,7 @@ module mom_cap_mod
use ESMF, only: ESMF_COORDSYS_SPH_DEG, ESMF_GridCreate, ESMF_INDEX_DELOCAL
use ESMF, only: ESMF_MESHLOC_ELEMENT, ESMF_RC_VAL_OUTOFRANGE, ESMF_StateGet
use ESMF, only: ESMF_TimePrint, ESMF_AlarmSet, ESMF_FieldGet
+use ESMF, only: operator(==), operator(/=), operator(+), operator(-)
! TODO ESMF_GridCompGetInternalState does not have an explicit Fortran interface.
!! Model does not compile with "use ESMF, only: ESMF_GridCompGetInternalState"
@@ -753,6 +756,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
integer :: userRc
character(len=512) :: restartfile ! Path/Name of restart file
character(len=*), parameter :: subname='(mom_cap:InitializeAdvertise)'
+ character(len=32) :: calendar
!--------------------------------
rc = ESMF_SUCCESS
@@ -804,7 +808,35 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
call fms_init(mpi_comm_mom)
call constants_init
call field_manager_init
- call set_calendar_type (JULIAN)
+
+ ! determine the calendar
+ if (cesm_coupled) then
+ call NUOPC_CompAttributeGet(gcomp, name="calendar", value=cvalue, &
+ isPresent=isPresent, isSet=isSet, rc=rc)
+ if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
+ line=__LINE__, &
+ file=__FILE__)) &
+ return ! bail out
+ if (isPresent .and. isSet) then
+ read(cvalue,*) calendar
+ select case (trim(calendar))
+ case ("NO_LEAP")
+ call set_calendar_type (NOLEAP)
+ case ("GREGORIAN")
+ call set_calendar_type (GREGORIAN)
+ case default
+ call ESMF_LogSetError(ESMF_RC_ARG_BAD, &
+ msg=subname//": Calendar not supported in MOM6: "//trim(calendar), &
+ line=__LINE__, file=__FILE__, rcToReturn=rc)
+ end select
+ else
+ call set_calendar_type (NOLEAP)
+ endif
+
+ else
+ call set_calendar_type (JULIAN)
+ endif
+
call diag_manager_init
! this ocean connector will be driven at set interval
@@ -960,6 +992,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
Ice_ocean_boundary% sw_flux_nir_dif (isc:iec,jsc:jec), &
Ice_ocean_boundary% lprec (isc:iec,jsc:jec), &
Ice_ocean_boundary% fprec (isc:iec,jsc:jec), &
+ Ice_ocean_boundary% seaice_melt_heat (isc:iec,jsc:jec),&
+ Ice_ocean_boundary% seaice_melt (isc:iec,jsc:jec), &
Ice_ocean_boundary% mi (isc:iec,jsc:jec), &
Ice_ocean_boundary% p (isc:iec,jsc:jec), &
Ice_ocean_boundary% runoff (isc:iec,jsc:jec), &
@@ -981,6 +1015,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
Ice_ocean_boundary%sw_flux_nir_dif = 0.0
Ice_ocean_boundary%lprec = 0.0
Ice_ocean_boundary%fprec = 0.0
+ Ice_ocean_boundary%seaice_melt = 0.0
+ Ice_ocean_boundary%seaice_melt_heat= 0.0
Ice_ocean_boundary%mi = 0.0
Ice_ocean_boundary%p = 0.0
Ice_ocean_boundary%runoff = 0.0
@@ -1030,8 +1066,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
call fld_list_add(fldsToOcn_num, fldsToOcn, "inst_pres_height_surface" , "will provide")
call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_rofl" , "will provide") !-> liquid runoff
call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_rofi" , "will provide") !-> ice runoff
- !call fld_list_add(fldsToOcn_num, fldsToOcn, "seaice_melt_water" , "will provide")
- !call fld_list_add(fldsToOcn_num, fldsToOcn, "seaice_melt_heat" , "will provide")
+ call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_fresh_water_to_ocean_rate", "will provide")
+ call fld_list_add(fldsToOcn_num, fldsToOcn, "net_heat_flx_to_ocn" , "will provide")
!call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_runoff_rate" , "will provide")
!call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_calving_rate" , "will provide")
@@ -1728,6 +1764,7 @@ subroutine ModelAdvance(gcomp, rc)
! local variables
integer :: userRc
logical :: existflag, isPresent, isSet
+ logical :: do_advance = .true.
type(ESMF_Clock) :: clock!< ESMF Clock class definition
type(ESMF_Alarm) :: alarm
type(ESMF_State) :: importState, exportState
@@ -1810,8 +1847,47 @@ subroutine ModelAdvance(gcomp, rc)
file=__FILE__)) &
return ! bail out
- Time = esmf2fms_time(currTime)
Time_step_coupled = esmf2fms_time(timeStep)
+ Time = esmf2fms_time(currTime)
+
+ !---------------
+ ! Apply ocean lag for startup runs:
+ !---------------
+
+ if (cesm_coupled) then
+ if (trim(runtype) == "initial") then
+
+ ! Do not call MOM6 timestepping routine if the first cpl tstep of a startup run
+ if (currTime == startTime) then
+ call ESMF_LogWrite("MOM6 - Skipping the first coupling timestep", ESMF_LOGMSG_INFO, rc=rc)
+ if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
+ line=__LINE__, &
+ file=__FILE__)) &
+ return ! bail out
+ do_advance = .false.
+ else
+ do_advance = .true.
+ endif
+
+ ! If the second cpl tstep of a startup run, step back a cpl tstep and advance for two cpl tsteps
+ if (currTime == startTime + timeStep) then
+ call ESMF_LogWrite("MOM6 - Stepping back one coupling timestep", ESMF_LOGMSG_INFO, rc=rc)
+ if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
+ line=__LINE__, &
+ file=__FILE__)) &
+ return ! bail out
+ Time = esmf2fms_time(currTime-timeStep) ! i.e., startTime
+
+ call ESMF_LogWrite("MOM6 - doubling the coupling timestep", ESMF_LOGMSG_INFO, rc=rc)
+ if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
+ line=__LINE__, &
+ file=__FILE__)) &
+ return ! bail out
+ Time_step_coupled = 2 * esmf2fms_time(timeStep)
+ endif
+ endif
+ endif
+
!---------------
! Write diagnostics for import
@@ -1854,7 +1930,9 @@ subroutine ModelAdvance(gcomp, rc)
!---------------
if(profile_memory) call ESMF_VMLogMemInfo("Entering MOM update_ocean_model: ")
- call update_ocean_model(Ice_ocean_boundary, ocean_state, ocean_public, Time, Time_step_coupled)
+ if (do_advance) then
+ call update_ocean_model(Ice_ocean_boundary, ocean_state, ocean_public, Time, Time_step_coupled)
+ endif
if(profile_memory) call ESMF_VMLogMemInfo("Leaving MOM update_ocean_model: ")
!---------------
diff --git a/config_src/nuopc_driver/mom_cap_methods.F90 b/config_src/nuopc_driver/mom_cap_methods.F90
index d893685aec..e6bdbea307 100644
--- a/config_src/nuopc_driver/mom_cap_methods.F90
+++ b/config_src/nuopc_driver/mom_cap_methods.F90
@@ -294,6 +294,7 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary,
!----
! salt flux from ice
!----
+ ice_ocean_boundary%salt_flux(:,:) = 0._ESMF_KIND_R8
call state_getimport(importState, 'mean_salt_rate', &
isc, iec, jsc, jec, ice_ocean_boundary%salt_flux,rc=rc)
if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
@@ -304,22 +305,24 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary,
! !----
! ! snow&ice melt heat flux (W/m^2)
! !----
- ! call state_getimport(importState, 'seaice_melt_heat', &
- ! isc, iec, jsc, jec, ice_ocean_boundary%seaice_melt_heat,rc=rc)
- ! if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
- ! line=__LINE__, &
- ! file=__FILE__)) &
- ! return ! bail out
+ ice_ocean_boundary%seaice_melt_heat(:,:) = 0._ESMF_KIND_R8
+ call state_getimport(importState, 'net_heat_flx_to_ocn', &
+ isc, iec, jsc, jec, ice_ocean_boundary%seaice_melt_heat,rc=rc)
+ if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
+ line=__LINE__, &
+ file=__FILE__)) &
+ return ! bail out
- ! !----
- ! ! snow&ice melt water flux (W/m^2)
- ! !----
- ! call state_getimport(importState, 'seaice_melt_water', &
- ! isc, iec, jsc, jec, ice_ocean_boundary%seaice_melt_water,rc=rc)
- ! if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
- ! line=__LINE__, &
- ! file=__FILE__)) &
- ! return ! bail out
+ ! !----
+ ! ! snow&ice melt water flux (W/m^2)
+ ! !----
+ ice_ocean_boundary%seaice_melt(:,:) = 0._ESMF_KIND_R8
+ call state_getimport(importState, 'mean_fresh_water_to_ocean_rate', &
+ isc, iec, jsc, jec, ice_ocean_boundary%seaice_melt,rc=rc)
+ if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
+ line=__LINE__, &
+ file=__FILE__)) &
+ return ! bail out
!----
! mass of overlying ice
@@ -373,7 +376,6 @@ subroutine mom_export(ocean_public, ocean_grid, ocean_state, exportState, clock,
rc = ESMF_SUCCESS
- ! Use Adcroft's rule of reciprocals; it does the right thing here.
call ESMF_ClockGet( clock, timeStep=timeStep, rc=rc)
if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
line=__LINE__, &
@@ -386,6 +388,7 @@ subroutine mom_export(ocean_public, ocean_grid, ocean_state, exportState, clock,
file=__FILE__)) &
return ! bail out
+ ! Use Adcroft's rule of reciprocals; it does the right thing here.
if (real(dt_int) > 0.0) then
inv_dt_int = 1.0 / real(dt_int)
else
diff --git a/src/core/MOM.F90 b/src/core/MOM.F90
index de7f01421d..301969ed50 100644
--- a/src/core/MOM.F90
+++ b/src/core/MOM.F90
@@ -50,6 +50,7 @@ module MOM
use MOM_ALE, only : ALE_init, ALE_end, ALE_main, ALE_CS, adjustGridForIntegrity
use MOM_ALE, only : ALE_getCoordinate, ALE_getCoordinateUnits, ALE_writeCoordinateFile
use MOM_ALE, only : ALE_updateVerticalGridType, ALE_remap_init_conds, ALE_register_diags
+use MOM_barotropic, only : Barotropic_CS
use MOM_boundary_update, only : call_OBC_register, OBC_register_end, update_OBC_CS
use MOM_coord_initialization, only : MOM_initialize_coord
use MOM_diabatic_driver, only : diabatic, diabatic_driver_init, diabatic_CS
@@ -321,7 +322,8 @@ module MOM
!< Pointer to the control structure for the MEKE updates
type(VarMix_CS), pointer :: VarMix => NULL()
!< Pointer to the control structure for the variable mixing module
-
+ type(Barotropic_CS), pointer :: Barotropic_CSp => NULL()
+ !< Pointer to the control structure for the barotropic module
type(tracer_registry_type), pointer :: tracer_Reg => NULL()
!< Pointer to the MOM tracer registry
type(tracer_advect_CS), pointer :: tracer_adv_CSp => NULL()
@@ -963,8 +965,8 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
call step_MOM_dyn_split_RK2(u, v, h, CS%tv, CS%visc, Time_local, dt, forces, &
p_surf_begin, p_surf_end, CS%uh, CS%vh, CS%uhtr, CS%vhtr, &
- CS%eta_av_bc, G, GV, US, CS%dyn_split_RK2_CSp, calc_dtbt, CS%VarMix, CS%MEKE,&
- waves=waves)
+ CS%eta_av_bc, G, GV, US, CS%dyn_split_RK2_CSp, calc_dtbt, CS%VarMix, &
+ CS%MEKE, CS%thickness_diffuse_CSp, waves=waves)
if (showCallTree) call callTree_waypoint("finished step_MOM_dyn_split (step_MOM)")
elseif (CS%do_dynamics) then ! ------------------------------------ not SPLIT
@@ -1157,7 +1159,9 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, &
call enable_averaging(dtdia, Time_end_thermo, CS%diag)
- call apply_oda_tracer_increments(dtdia,G,tv,h,CS%odaCS)
+ if (associated(CS%odaCS)) then
+ call apply_oda_tracer_increments(dtdia,G,tv,h,CS%odaCS)
+ endif
if (update_BBL) then
! Calculate the BBL properties and store them inside visc (u,h).
@@ -2294,11 +2298,13 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
call VarMix_init(Time, G, GV, US, param_file, diag, CS%VarMix)
call set_visc_init(Time, G, GV, US, param_file, diag, CS%visc, CS%set_visc_CSp, restart_CSp, CS%OBC)
+ call thickness_diffuse_init(Time, G, GV, US, param_file, diag, CS%CDp, CS%thickness_diffuse_CSp)
if (CS%split) then
allocate(eta(SZI_(G),SZJ_(G))) ; eta(:,:) = 0.0
call initialize_dyn_split_RK2(CS%u, CS%v, CS%h, CS%uh, CS%vh, eta, Time, &
G, GV, US, param_file, diag, CS%dyn_split_RK2_CSp, restart_CSp, &
CS%dt, CS%ADp, CS%CDp, MOM_internal_state, CS%VarMix, CS%MEKE, &
+ CS%thickness_diffuse_CSp, &
CS%OBC, CS%update_OBC_CSp, CS%ALE_CSp, CS%set_visc_CSp, &
CS%visc, dirs, CS%ntrunc, calc_dtbt=calc_dtbt)
if (CS%dtbt_reset_period > 0.0) then
@@ -2326,7 +2332,6 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
endif
call callTree_waypoint("dynamics initialized (initialize_MOM)")
- call thickness_diffuse_init(Time, G, GV, US, param_file, diag, CS%CDp, CS%thickness_diffuse_CSp)
CS%mixedlayer_restrat = mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, &
CS%mixedlayer_restrat_CSp, restart_CSp)
if (CS%mixedlayer_restrat) then
diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90
index 2f1eb68961..33450e8a3d 100644
--- a/src/core/MOM_barotropic.F90
+++ b/src/core/MOM_barotropic.F90
@@ -56,7 +56,7 @@ module MOM_barotropic
#endif
public btcalc, bt_mass_source, btstep, barotropic_init, barotropic_end
-public register_barotropic_restarts, set_dtbt
+public register_barotropic_restarts, set_dtbt, barotropic_get_tav
! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
@@ -4363,6 +4363,29 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS,
end subroutine barotropic_init
+!> Copies ubtav and vbtav from private type into arrays
+subroutine barotropic_get_tav(CS, ubtav, vbtav, G)
+ type(barotropic_CS), pointer :: CS !< Control structure for
+ !! this module
+ type(ocean_grid_type), intent(in) :: G !< Grid structure
+ real, dimension(SZIB_(G),SZJ_(G)), intent(inout) :: ubtav!< zonal barotropic vel.
+ !! ave. over baroclinic time-step (m s-1)
+ real, dimension(SZI_(G),SZJB_(G)), intent(inout) :: vbtav!< meridional barotropic vel.
+ !! ave. over baroclinic time-step (m s-1)
+ ! Local variables
+ integer :: i,j
+
+ do j = G%jsc, G%jec ; do I = G%isc-1, G%iec
+ ubtav(I,j) = CS%ubtav(I,j)
+ enddo ; enddo
+
+ do J = G%jsc-1, G%jec ; do i = G%isc, G%iec
+ vbtav(i,J) = CS%vbtav(i,J)
+ enddo ; enddo
+
+end subroutine barotropic_get_tav
+
+
!> Clean up the barotropic control structure.
subroutine barotropic_end(CS)
type(barotropic_CS), pointer :: CS !< Control structure to clear out.
diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90
index d862fae71d..c6bc7b5c6a 100644
--- a/src/core/MOM_dynamics_split_RK2.F90
+++ b/src/core/MOM_dynamics_split_RK2.F90
@@ -52,6 +52,7 @@ module MOM_dynamics_split_RK2
use MOM_open_boundary, only : open_boundary_test_extern_h
use MOM_PressureForce, only : PressureForce, PressureForce_init, PressureForce_CS
use MOM_set_visc, only : set_viscous_ML, set_visc_CS
+use MOM_thickness_diffuse, only : thickness_diffuse_CS
use MOM_tidal_forcing, only : tidal_forcing_init, tidal_forcing_CS
use MOM_unit_scaling, only : unit_scale_type
use MOM_vert_friction, only : vertvisc, vertvisc_coef, vertvisc_remnant
@@ -184,6 +185,8 @@ module MOM_dynamics_split_RK2
type(PressureForce_CS), pointer :: PressureForce_CSp => NULL()
!> A pointer to the barotropic stepping control structure
type(barotropic_CS), pointer :: barotropic_CSp => NULL()
+ !> A pointer to a structure containing interface height diffusivities
+ type(thickness_diffuse_CS), pointer :: thickness_diffuse_CSp => NULL()
!> A pointer to the vertical viscosity control structure
type(vertvisc_CS), pointer :: vertvisc_CSp => NULL()
!> A pointer to the set_visc control structure
@@ -230,7 +233,7 @@ module MOM_dynamics_split_RK2
subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, &
Time_local, dt, forces, p_surf_begin, p_surf_end, &
uh, vh, uhtr, vhtr, eta_av, &
- G, GV, US, CS, calc_dtbt, VarMix, MEKE, Waves)
+ G, GV, US, CS, calc_dtbt, VarMix, MEKE, thickness_diffuse_CSp, Waves)
type(ocean_grid_type), intent(inout) :: G !< ocean grid structure
type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
@@ -267,8 +270,12 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, &
logical, intent(in) :: calc_dtbt !< if true, recalculate barotropic time step
type(VarMix_CS), pointer :: VarMix !< specify the spatially varying viscosities
type(MEKE_type), pointer :: MEKE !< related to mesoscale eddy kinetic energy param
+ type(thickness_diffuse_CS), pointer :: thickness_diffuse_CSp!< Pointer to a structure containing
+ !! interface height diffusivities
type(wave_parameters_CS), optional, pointer :: Waves !< A pointer to a structure containing
!! fields related to the surface wave conditions
+
+ ! local variables
real :: dt_pred ! The time step for the predictor part of the baroclinic time stepping.
real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: up ! Predicted zonal velocity [m s-1].
@@ -682,7 +689,8 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, &
! diffu = horizontal viscosity terms (u_av)
call cpu_clock_begin(id_clock_horvisc)
call horizontal_viscosity(u_av, v_av, h_av, CS%diffu, CS%diffv, &
- MEKE, Varmix, G, GV, US, CS%hor_visc_CSp, OBC=CS%OBC)
+ MEKE, Varmix, G, GV, US, CS%hor_visc_CSp, &
+ OBC=CS%OBC, BT=CS%barotropic_CSp)
call cpu_clock_end(id_clock_horvisc)
if (showCallTree) call callTree_wayPoint("done with horizontal_viscosity (step_MOM_dyn_split_RK2)")
@@ -953,7 +961,8 @@ end subroutine register_restarts_dyn_split_RK2
!! dynamic core, including diagnostics and the cpu clocks.
subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param_file, &
diag, CS, restart_CS, dt, Accel_diag, Cont_diag, MIS, &
- VarMix, MEKE, OBC, update_OBC_CSp, ALE_CSp, setVisc_CSp, &
+ VarMix, MEKE, thickness_diffuse_CSp, &
+ OBC, update_OBC_CSp, ALE_CSp, setVisc_CSp, &
visc, dirs, ntrunc, calc_dtbt)
type(ocean_grid_type), intent(inout) :: G !< ocean grid structure
type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
@@ -981,6 +990,10 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
!! diagnostic pointers
type(VarMix_CS), pointer :: VarMix !< points to spatially variable viscosities
type(MEKE_type), pointer :: MEKE !< points to mesoscale eddy kinetic energy fields
+! type(Barotropic_CS), pointer :: Barotropic_CSp !< Pointer to the control structure for
+! !! the barotropic module
+ type(thickness_diffuse_CS), pointer :: thickness_diffuse_CSp !< Pointer to the control structure
+ !! used for the isopycnal height diffusive transport.
type(ocean_OBC_type), pointer :: OBC !< points to OBC related fields
type(update_OBC_CS), pointer :: update_OBC_CSp !< points to OBC update related fields
type(ALE_CS), pointer :: ALE_CSp !< points to ALE control structure
@@ -992,6 +1005,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
!! truncated (this should be 0).
logical, intent(out) :: calc_dtbt !< If true, recalculate the barotropic time step
+ ! local variables
real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_tmp
character(len=40) :: mdl = "MOM_dynamics_split_RK2" ! This module's name.
character(len=48) :: thickness_units, flux_units, eta_rest_name
@@ -1136,7 +1150,8 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
if (.not. query_initialized(CS%diffu,"diffu",restart_CS) .or. &
.not. query_initialized(CS%diffv,"diffv",restart_CS)) &
call horizontal_viscosity(u, v, h, CS%diffu, CS%diffv, MEKE, VarMix, &
- G, GV, US, CS%hor_visc_CSp, OBC=CS%OBC)
+ G, GV, US, CS%hor_visc_CSp, &
+ OBC=CS%OBC, BT=CS%barotropic_CSp)
if (.not. query_initialized(CS%u_av,"u2", restart_CS) .or. &
.not. query_initialized(CS%u_av,"v2", restart_CS)) then
CS%u_av(:,:,:) = u(:,:,:)
diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90
index 0995725536..dd03e11f42 100644
--- a/src/core/MOM_dynamics_unsplit.F90
+++ b/src/core/MOM_dynamics_unsplit.F90
@@ -50,7 +50,6 @@ module MOM_dynamics_unsplit
!* *
!********+*********+*********+*********+*********+*********+*********+**
-
use MOM_variables, only : vertvisc_type, thermo_var_ptrs
use MOM_variables, only : accel_diag_ptrs, ocean_internal_state, cont_diag_ptrs
use MOM_forcing_type, only : mech_forcing
@@ -76,6 +75,7 @@ module MOM_dynamics_unsplit
use MOM_time_manager, only : operator(-), operator(>), operator(*), operator(/)
use MOM_ALE, only : ALE_CS
+use MOM_barotropic, only : barotropic_CS
use MOM_boundary_update, only : update_OBC_data, update_OBC_CS
use MOM_continuity, only : continuity, continuity_init, continuity_CS
use MOM_CoriolisAdv, only : CorAdCalc, CoriolisAdv_init, CoriolisAdv_CS
@@ -91,6 +91,7 @@ module MOM_dynamics_unsplit
use MOM_open_boundary, only : open_boundary_zero_normal_flow
use MOM_PressureForce, only : PressureForce, PressureForce_init, PressureForce_CS
use MOM_set_visc, only : set_viscous_ML, set_visc_CS
+use MOM_thickness_diffuse, only : thickness_diffuse_CS
use MOM_tidal_forcing, only : tidal_forcing_init, tidal_forcing_CS
use MOM_unit_scaling, only : unit_scale_type
use MOM_vert_friction, only : vertvisc, vertvisc_coef
diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90
index be81a8b25e..b5b547b362 100644
--- a/src/core/MOM_dynamics_unsplit_RK2.F90
+++ b/src/core/MOM_dynamics_unsplit_RK2.F90
@@ -75,6 +75,7 @@ module MOM_dynamics_unsplit_RK2
use MOM_ALE, only : ALE_CS
use MOM_boundary_update, only : update_OBC_data, update_OBC_CS
+use MOM_barotropic, only : barotropic_CS
use MOM_continuity, only : continuity, continuity_init, continuity_CS
use MOM_CoriolisAdv, only : CorAdCalc, CoriolisAdv_init, CoriolisAdv_CS
use MOM_debugging, only : check_redundant
@@ -88,6 +89,7 @@ module MOM_dynamics_unsplit_RK2
use MOM_open_boundary, only : open_boundary_zero_normal_flow
use MOM_PressureForce, only : PressureForce, PressureForce_init, PressureForce_CS
use MOM_set_visc, only : set_viscous_ML, set_visc_CS
+use MOM_thickness_diffuse, only : thickness_diffuse_CS
use MOM_tidal_forcing, only : tidal_forcing_init, tidal_forcing_CS
use MOM_unit_scaling, only : unit_scale_type
use MOM_vert_friction, only : vertvisc, vertvisc_coef
@@ -230,7 +232,6 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt,
type(MEKE_type), pointer :: MEKE !< A pointer to a structure containing
!! fields related to the Mesoscale
!! Eddy Kinetic Energy.
-
! Local variables
real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_av, hp
real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: up
diff --git a/src/diagnostics/MOM_obsolete_params.F90 b/src/diagnostics/MOM_obsolete_params.F90
index 797db75240..0e563648f5 100644
--- a/src/diagnostics/MOM_obsolete_params.F90
+++ b/src/diagnostics/MOM_obsolete_params.F90
@@ -69,6 +69,9 @@ subroutine find_obsolete_params(param_file)
hint="Instead use OBC_SEGMENT_XXX_DATA.")
call obsolete_char(param_file, "EXTEND_OBC_SEGMENTS", &
hint="This option is no longer needed, nor supported.")
+ call obsolete_char(param_file, "MEKE_VISCOSITY_COEFF", &
+ hint="This option has been replaced by MEKE_VISCOSITY_COEFF_KU and \n" //&
+ " MEKE_VISCOSITY_COEFF_AU. Please set these parameters instead.")
nseg = 0
call read_param(param_file, "OBC_NUMBER_OF_SEGMENTS", nseg)
do l_seg = 1,nseg
diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90
index 5020a4cbe7..e8f1fecf60 100644
--- a/src/ice_shelf/MOM_ice_shelf.F90
+++ b/src/ice_shelf/MOM_ice_shelf.F90
@@ -321,7 +321,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces)
if (CS%mass_from_file) call update_shelf_mass(G, CS, ISS, Time)
endif
- if (CS%DEBUG) then
+ if (CS%debug) then
call hchksum(fluxes%frac_shelf_h, "frac_shelf_h before apply melting", G%HI, haloshift=0)
call hchksum(state%sst, "sst before apply melting", G%HI, haloshift=0)
call hchksum(state%sss, "sss before apply melting", G%HI, haloshift=0)
@@ -650,7 +650,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces)
endif
endif
- if (CS%DEBUG) call MOM_forcing_chksum("Before add shelf flux", fluxes, G, US, haloshift=0)
+ if (CS%debug) call MOM_forcing_chksum("Before add shelf flux", fluxes, G, CS%US, haloshift=0)
call add_shelf_flux(G, CS, state, fluxes)
@@ -692,7 +692,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces)
call cpu_clock_end(id_clock_shelf)
- if (CS%DEBUG) call MOM_forcing_chksum("End of shelf calc flux", fluxes, G, US, haloshift=0)
+ if (CS%debug) call MOM_forcing_chksum("End of shelf calc flux", fluxes, G, CS%US, haloshift=0)
end subroutine shelf_calc_flux
@@ -1055,8 +1055,8 @@ subroutine add_shelf_flux(G, CS, state, fluxes)
endif
enddo ; enddo
- if (CS%DEBUG) then
- write(mesg,*) 'Mean melt flux [kg m-2 s-1], dt = ', mean_melt_flux, CS%time_step
+ if (CS%debug) then
+ write(mesg,*) 'Mean melt flux (kg/(m^2 s)), dt = ', mean_melt_flux, CS%time_step
call MOM_mesg(mesg)
call MOM_forcing_chksum("After constant sea level", fluxes, G, CS%US, haloshift=0)
endif
@@ -1514,7 +1514,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl
if (G%areaT(i,j) > 0.0) fluxes%frac_shelf_h(i,j) = ISS%area_shelf_h(i,j) / G%areaT(i,j)
enddo ; enddo ; endif
- if (CS%DEBUG) then
+ if (CS%debug) then
call hchksum(fluxes%frac_shelf_h, "IS init: frac_shelf_h", G%HI, haloshift=0)
endif
diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90
index b1c970871b..415ae3d813 100644
--- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90
+++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90
@@ -928,7 +928,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u, v, iters, time)
call ice_shelf_solve_inner(CS, ISS, G, u, v, TAUDX, TAUDY, H_node, float_cond, &
ISS%hmask, conv_flag, iters, time, Phi, Phisub)
- if (CS%DEBUG) then
+ if (CS%debug) then
call qchksum(u, "u shelf", G%HI, haloshift=2)
call qchksum(v, "v shelf", G%HI, haloshift=2)
endif
@@ -3600,7 +3600,7 @@ subroutine ice_shelf_temp(CS, ISS, G, US, time_step, melt_rate, Time)
call pass_var(CS%t_shelf, G%domain)
call pass_var(CS%tmask, G%domain)
- if (CS%DEBUG) then
+ if (CS%debug) then
call hchksum(CS%t_shelf, "temp after front", G%HI, haloshift=3)
endif
diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90
index 447f2eefc0..a17bfc6aa9 100644
--- a/src/parameterizations/lateral/MOM_MEKE.F90
+++ b/src/parameterizations/lateral/MOM_MEKE.F90
@@ -11,6 +11,7 @@ module MOM_MEKE
use MOM_diag_mediator, only : diag_ctrl, time_type
use MOM_domains, only : create_group_pass, do_group_pass
use MOM_domains, only : group_pass_type
+use MOM_domains, only : pass_var, pass_vector
use MOM_error_handler, only : MOM_error, FATAL, WARNING, NOTE, MOM_mesg
use MOM_file_parser, only : read_param, get_param, log_version, param_file_type
use MOM_grid, only : ocean_grid_type
@@ -33,6 +34,7 @@ module MOM_MEKE
! Parameters
real :: MEKE_FrCoeff !< Efficiency of conversion of ME into MEKE [nondim]
real :: MEKE_GMcoeff !< Efficiency of conversion of PE into MEKE [nondim]
+ real :: MEKE_GMECoeff !< Efficiency of conversion of MEKE into ME by GME [nondim]
real :: MEKE_damping !< Local depth-independent MEKE dissipation rate [s-1].
real :: MEKE_Cd_scale !< The ratio of the bottom eddy velocity to the column mean
!! eddy velocity, i.e. sqrt(2*MEKE). This should be less than 1
@@ -41,6 +43,11 @@ module MOM_MEKE
real :: MEKE_min_gamma!< Minimum value of gamma_b^2 allowed [nondim]
real :: MEKE_Ct !< Coefficient in the \f$\gamma_{bt}\f$ expression [nondim]
logical :: visc_drag !< If true use the vertvisc_type to calculate bottom drag.
+ logical :: Jansen15_drag !< If true use the bottom drag formulation from Jansen et al. (2015)
+ logical :: MEKE_GEOMETRIC !< If true, uses the GM coefficient formulation from the GEOMETRIC
+ !! framework (Marshall et al., 2012)
+ logical :: GM_src_alt !< If true, use the GM energy conversion form S^2*N^2*kappa rather
+ !! than the streamfunction for the MEKE GM source term.
logical :: Rd_as_max_scale !< If true the length scale can not exceed the
!! first baroclinic deformation radius.
logical :: use_old_lscale !< Use the old formula for mixing length scale.
@@ -54,8 +61,11 @@ module MOM_MEKE
real :: MEKE_K4 !< Background bi-harmonic diffusivity (of MEKE) [m4 s-1]
real :: KhMEKE_Fac !< A factor relating MEKE%Kh to the diffusivity used for
!! MEKE itself [nondim].
- real :: viscosity_coeff !< The scaling coefficient in the expression for
- !! viscosity used to parameterize lateral momentum mixing
+ real :: viscosity_coeff_Ku !< The scaling coefficient in the expression for
+ !! viscosity used to parameterize lateral harmonic momentum mixing
+ !! by unresolved eddies represented by MEKE.
+ real :: viscosity_coeff_Au !< The scaling coefficient in the expression for
+ !! viscosity used to parameterize lateral biharmonic momentum mixing
!! by unresolved eddies represented by MEKE.
real :: Lfixed !< Fixed mixing length scale [m].
real :: aDeform !< Weighting towards deformation scale of mixing length [nondim]
@@ -77,8 +87,8 @@ module MOM_MEKE
!>@{ Diagnostic handles
integer :: id_MEKE = -1, id_Ue = -1, id_Kh = -1, id_src = -1
integer :: id_Ub = -1, id_Ut = -1
- integer :: id_GM_src = -1, id_mom_src = -1, id_decay = -1
- integer :: id_KhMEKE_u = -1, id_KhMEKE_v = -1, id_Ku = -1
+ integer :: id_GM_src = -1, id_mom_src = -1, id_GME_snk = -1, id_decay = -1
+ integer :: id_KhMEKE_u = -1, id_KhMEKE_v = -1, id_Ku = -1, id_Au = -1
integer :: id_Le = -1, id_gamma_b = -1, id_gamma_t = -1
integer :: id_Lrhines = -1, id_Leady = -1
!!@}
@@ -87,7 +97,9 @@ module MOM_MEKE
integer :: id_clock_pass !< Clock for group pass calls
type(group_pass_type) :: pass_MEKE !< Type for group halo pass calls
type(group_pass_type) :: pass_Kh !< Type for group halo pass calls
+ type(group_pass_type) :: pass_Kh_diff !< Type for group halo pass calls
type(group_pass_type) :: pass_Ku !< Type for group halo pass calls
+ type(group_pass_type) :: pass_Au !< Type for group halo pass calls
type(group_pass_type) :: pass_del2MEKE !< Type for group halo pass calls
end type MEKE_CS
@@ -101,8 +113,8 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure.
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2].
- real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: SN_u !< Eady growth rate at u-points [s-1].
- real, dimension(SZI_(G),SZJB_(G)), intent(in) :: SN_v !< Eady growth rate at v-points [s-1].
+ real, dimension(SZIB_(G),SZJ_(G)), intent(inout) :: SN_u !< Eady growth rate at u-points [s-1].
+ real, dimension(SZI_(G),SZJB_(G)), intent(inout) :: SN_v !< Eady growth rate at v-points [s-1].
type(vertvisc_type), intent(in) :: visc !< The vertical viscosity type.
real, intent(in) :: dt !< Model(baroclinic) time-step [s].
type(MEKE_CS), pointer :: CS !< MEKE control structure.
@@ -117,11 +129,13 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
MEKE_decay, & ! The MEKE decay timescale [s-1].
MEKE_GM_src, & ! The MEKE source from thickness mixing [m2 s-3].
MEKE_mom_src, & ! The MEKE source from momentum [m2 s-3].
+ MEKE_GME_snk, & ! The MEKE sink from GME backscatter [m2 s-3].
drag_rate_visc, &
drag_rate, & ! The MEKE spindown timescale due to bottom drag [s-1].
LmixScale, & ! Square of eddy mixing length [m2].
barotrFac2, & ! Ratio of EKE_barotropic / EKE [nondim]
bottomFac2 ! Ratio of EKE_bottom / EKE [nondim]
+
real, dimension(SZIB_(G),SZJ_(G)) :: &
MEKE_uflux, & ! The zonal diffusive flux of MEKE [kg m2 s-3].
Kh_u, & ! The zonal diffusivity that is actually used [m2 s-1].
@@ -168,6 +182,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
if (CS%debug) then
if (associated(MEKE%mom_src)) call hchksum(MEKE%mom_src, 'MEKE mom_src',G%HI)
+ if (associated(MEKE%GME_snk)) call hchksum(MEKE%GME_snk, 'MEKE GME_snk',G%HI)
if (associated(MEKE%GM_src)) call hchksum(MEKE%GM_src, 'MEKE GM_src',G%HI)
if (associated(MEKE%MEKE)) call hchksum(MEKE%MEKE, 'MEKE MEKE',G%HI)
call uvchksum("MEKE SN_[uv]", SN_u, SN_v, G%HI)
@@ -282,13 +297,26 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
enddo ; enddo
endif
- if (associated(MEKE%GM_src)) then
- !$OMP parallel do default(shared)
+ if (associated(MEKE%GME_snk)) then
+!$OMP do
do j=js,je ; do i=is,ie
- src(i,j) = src(i,j) - CS%MEKE_GMcoeff*I_mass(i,j)*MEKE%GM_src(i,j)
+ src(i,j) = src(i,j) - CS%MEKE_GMECoeff*I_mass(i,j)*MEKE%GME_snk(i,j)
enddo ; enddo
endif
+ if (associated(MEKE%GM_src)) then
+!$OMP do
+ if (CS%GM_src_alt) then
+ do j=js,je ; do i=is,ie
+ src(i,j) = src(i,j) - CS%MEKE_GMcoeff*MEKE%GM_src(i,j) / MAX(1.0,G%bathyT(i,j))
+ enddo ; enddo
+ else
+ do j=js,je ; do i=is,ie
+ src(i,j) = src(i,j) - CS%MEKE_GMcoeff*I_mass(i,j)*MEKE%GM_src(i,j)
+ enddo ; enddo
+ endif
+ endif
+
! Increase EKE by a full time-steps worth of source
!$OMP parallel do default(shared)
do j=js,je ; do i=is,ie
@@ -297,24 +325,39 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
if (use_drag_rate) then
! Calculate a viscous drag rate (includes BBL contributions from mean flow and eddies)
- !$OMP parallel do default(shared)
- do j=js,je ; do i=is,ie
- drag_rate(i,j) = (Rho0 * I_mass(i,j)) * sqrt( drag_rate_visc(i,j)**2 &
+!$OMP do
+ if (CS%Jansen15_drag) then
+ do j=js,je ; do i=is,ie
+ drag_rate(i,j) = (cdrag2/MAX(1.0,G%bathyT(i,j))) * sqrt(CS%MEKE_Uscale**2 + drag_rate_visc(i,j)**2 + &
+ 2.0*bottomFac2(i,j)*MEKE%MEKE(i,j)) * 2.0 * bottomFac2(i,j)*MEKE%MEKE(i,j)
+ enddo ; enddo
+ else
+ do j=js,je ; do i=is,ie
+ drag_rate(i,j) = (Rho0 * I_mass(i,j)) * sqrt( drag_rate_visc(i,j)**2 &
+ cdrag2 * ( max(0.0, 2.0*bottomFac2(i,j)*MEKE%MEKE(i,j)) + CS%MEKE_Uscale**2 ) )
- enddo ; enddo
+ enddo ; enddo
+ endif
endif
! First stage of Strang splitting
- !$OMP parallel do default(shared) private(ldamping)
- do j=js,je ; do i=is,ie
- ldamping = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j)
- if (MEKE%MEKE(i,j)<0.) ldamping = 0.
- ! notice that the above line ensures a damping only if MEKE is positive,
- ! while leaving MEKE unchanged if it is negative
- MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1.0 + sdt_damp*ldamping)
- MEKE_decay(i,j) = ldamping*G%mask2dT(i,j)
- enddo ; enddo
-
+!$OMP do
+ if (CS%Jansen15_drag) then
+ do j=js,je ; do i=is,ie
+ ldamping = CS%MEKE_damping + drag_rate(i,j)
+ MEKE%MEKE(i,j) = MEKE%MEKE(i,j) - MIN(MEKE%MEKE(i,j),sdt_damp*drag_rate(i,j))
+ MEKE_decay(i,j) = ldamping*G%mask2dT(i,j)
+ enddo ; enddo
+ else
+ do j=js,je ; do i=is,ie
+ ldamping = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j)
+ if (MEKE%MEKE(i,j)<0.) ldamping = 0.
+ ! notice that the above line ensures a damping only if MEKE is positive,
+ ! while leaving MEKE unchanged if it is negative
+ MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1.0 + sdt_damp*ldamping)
+ MEKE_decay(i,j) = ldamping*G%mask2dT(i,j)
+ enddo ; enddo
+ endif
+!$OMP end parallel
if (CS%MEKE_KH >= 0.0 .or. CS%KhMEKE_FAC > 0.0 .or. CS%MEKE_K4 >= 0.0) then
! Update halos for lateral or bi-harmonic diffusion
call cpu_clock_begin(CS%id_clock_pass)
@@ -392,6 +435,8 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
! Limit Kh to avoid CFL violations.
if (associated(MEKE%Kh)) &
Kh_here = max(0.,CS%MEKE_Kh) + CS%KhMEKE_Fac*0.5*(MEKE%Kh(i,j)+MEKE%Kh(i+1,j))
+ if (associated(MEKE%Kh_diff)) &
+ Kh_here = max(0.,CS%MEKE_Kh) + CS%KhMEKE_Fac*0.5*(MEKE%Kh_diff(i,j)+MEKE%Kh_diff(i+1,j))
Inv_Kh_max = 2.0*sdt * ((G%dy_Cu(I,j)*G%IdxCu(I,j)) * &
max(G%IareaT(i,j),G%IareaT(i+1,j)))
if (Kh_here*Inv_Kh_max > 0.25) Kh_here = 0.25 / Inv_Kh_max
@@ -405,6 +450,8 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
do J=js-1,je ; do i=is,ie
if (associated(MEKE%Kh)) &
Kh_here = max(0.,CS%MEKE_Kh) + CS%KhMEKE_Fac*0.5*(MEKE%Kh(i,j)+MEKE%Kh(i,j+1))
+ if (associated(MEKE%Kh_diff)) &
+ Kh_here = max(0.,CS%MEKE_Kh) + CS%KhMEKE_Fac*0.5*(MEKE%Kh_diff(i,j)+MEKE%Kh_diff(i,j+1))
Inv_Kh_max = 2.0*sdt * ((G%dx_Cv(i,J)*G%IdyCv(i,J)) * &
max(G%IareaT(i,j),G%IareaT(i,j+1)))
if (Kh_here*Inv_Kh_max > 0.25) Kh_here = 0.25 / Inv_Kh_max
@@ -454,65 +501,90 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
if (sdt>sdt_damp) then
! Recalculate the drag rate, since MEKE has changed.
if (use_drag_rate) then
- !$OMP parallel do default(shared)
- do j=js,je ; do i=is,ie
- drag_rate(i,j) = (Rho0 * I_mass(i,j)) * sqrt( drag_rate_visc(i,j)**2 &
- + cdrag2 * ( max(0.0, 2.0*bottomFac2(i,j)*MEKE%MEKE(i,j)) + CS%MEKE_Uscale**2 ) )
- enddo ; enddo
+!$OMP do
+ if (CS%Jansen15_drag) then
+ do j=js,je ; do i=is,ie
+ ldamping = CS%MEKE_damping + drag_rate(i,j)
+ MEKE%MEKE(i,j) = MEKE%MEKE(i,j) -sdt_damp*drag_rate(i,j)
+ MEKE_decay(i,j) = ldamping*G%mask2dT(i,j)
+ enddo ; enddo
+ else
+ do j=js,je ; do i=is,ie
+ drag_rate(i,j) = (Rho0 * I_mass(i,j)) * sqrt( drag_rate_visc(i,j)**2 &
+ + cdrag2 * ( max(0.0, 2.0*bottomFac2(i,j)*MEKE%MEKE(i,j)) + CS%MEKE_Uscale**2 ) )
+ enddo ; enddo
+ do j=js,je ; do i=is,ie
+ ldamping = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j)
+ if (MEKE%MEKE(i,j)<0.) ldamping = 0.
+ ! notice that the above line ensures a damping only if MEKE is positive,
+ ! while leaving MEKE unchanged if it is negative
+ MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1.0 + sdt_damp*ldamping)
+ MEKE_decay(i,j) = ldamping*G%mask2dT(i,j)
+ enddo ; enddo
+ endif
endif
- !$OMP parallel do default(shared) private(ldamping)
- do j=js,je ; do i=is,ie
- ldamping = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j)
- if (MEKE%MEKE(i,j)<0.) ldamping = 0.
- ! notice that the above line ensures a damping only if MEKE is positive,
- ! while leaving MEKE unchanged if it is negative
- MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1.0 + sdt_damp*ldamping)
- MEKE_decay(i,j) = 0.5 * G%mask2dT(i,j) * (MEKE_decay(i,j) + ldamping)
- enddo ; enddo
+!$OMP do
endif
endif ! MEKE_KH>=0
+ ! do j=js,je ; do i=is,ie
+ ! MEKE%MEKE(i,j) = MAX(MEKE%MEKE(i,j),0.0)
+ ! enddo ; enddo
+
+!$OMP end parallel
call cpu_clock_begin(CS%id_clock_pass)
call do_group_pass(CS%pass_MEKE, G%Domain)
call cpu_clock_end(CS%id_clock_pass)
! Calculate diffusivity for main model to use
if (CS%MEKE_KhCoeff>0.) then
- if (CS%use_old_lscale) then
- if (CS%Rd_as_max_scale) then
- !$OMP parallel do default(shared)
- do j=js,je ; do i=is,ie
- MEKE%Kh(i,j) = (CS%MEKE_KhCoeff &
+ if (.not.CS%MEKE_GEOMETRIC) then
+ if (CS%use_old_lscale) then
+ if (CS%Rd_as_max_scale) then
+ !$OMP parallel do default(shared)
+ do j=js,je ; do i=is,ie
+ MEKE%Kh(i,j) = (CS%MEKE_KhCoeff &
* sqrt(2.*max(0.,barotrFac2(i,j)*MEKE%MEKE(i,j))*G%areaT(i,j))) &
* min(MEKE%Rd_dx_h(i,j), 1.0)
- enddo ; enddo
+ enddo ; enddo
+ else
+ !$OMP parallel do default(shared)
+ do j=js,je ; do i=is,ie
+ MEKE%Kh(i,j) = CS%MEKE_KhCoeff*sqrt(2.*max(0.,barotrFac2(i,j)*MEKE%MEKE(i,j))*G%areaT(i,j))
+ enddo ; enddo
+ endif
else
!$OMP parallel do default(shared)
do j=js,je ; do i=is,ie
- MEKE%Kh(i,j) = CS%MEKE_KhCoeff*sqrt(2.*max(0.,barotrFac2(i,j)*MEKE%MEKE(i,j))*G%areaT(i,j))
+ MEKE%Kh(i,j) = (CS%MEKE_KhCoeff*sqrt(2.*max(0.,barotrFac2(i,j)*MEKE%MEKE(i,j)))*LmixScale(i,j))
enddo ; enddo
endif
- else
- !$OMP parallel do default(shared)
- do j=js,je ; do i=is,ie
- MEKE%Kh(i,j) = (CS%MEKE_KhCoeff*sqrt(2.*max(0.,barotrFac2(i,j)*MEKE%MEKE(i,j)))*LmixScale(i,j))
- enddo ; enddo
- endif
- call cpu_clock_begin(CS%id_clock_pass)
- call do_group_pass(CS%pass_Kh, G%Domain)
- call cpu_clock_end(CS%id_clock_pass)
+ call cpu_clock_begin(CS%id_clock_pass)
+ call do_group_pass(CS%pass_Kh, G%Domain)
+ call cpu_clock_end(CS%id_clock_pass)
+ endif
endif
! Calculate viscosity for the main model to use
- if (CS%viscosity_coeff/=0.) then
+ if (CS%viscosity_coeff_Ku /=0.) then
do j=js,je ; do i=is,ie
- MEKE%Ku(i,j) = CS%viscosity_coeff*sqrt(2.*max(0.,MEKE%MEKE(i,j)))*LmixScale(i,j)
+ MEKE%Ku(i,j) = CS%viscosity_coeff_Ku*sqrt(2.*max(0.,MEKE%MEKE(i,j)))*LmixScale(i,j)
enddo ; enddo
call cpu_clock_begin(CS%id_clock_pass)
call do_group_pass(CS%pass_Ku, G%Domain)
call cpu_clock_end(CS%id_clock_pass)
endif
+ if (CS%viscosity_coeff_Au /=0.) then
+ do j=js,je ; do i=is,ie
+ MEKE%Au(i,j) = CS%viscosity_coeff_Au*sqrt(2.*max(0.,MEKE%MEKE(i,j)))*LmixScale(i,j)**3
+ enddo ; enddo
+ call cpu_clock_begin(CS%id_clock_pass)
+ call do_group_pass(CS%pass_Au, G%Domain)
+ call cpu_clock_end(CS%id_clock_pass)
+ endif
+
+
! Offer fields for averaging.
if (CS%id_MEKE>0) call post_data(CS%id_MEKE, MEKE%MEKE, CS%diag)
if (CS%id_Ue>0) call post_data(CS%id_Ue, sqrt(max(0.,2.0*MEKE%MEKE)), CS%diag)
@@ -520,12 +592,14 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
if (CS%id_Ut>0) call post_data(CS%id_Ut, sqrt(max(0.,2.0*MEKE%MEKE*barotrFac2)), CS%diag)
if (CS%id_Kh>0) call post_data(CS%id_Kh, MEKE%Kh, CS%diag)
if (CS%id_Ku>0) call post_data(CS%id_Ku, MEKE%Ku, CS%diag)
+ if (CS%id_Au>0) call post_data(CS%id_Au, MEKE%Au, CS%diag)
if (CS%id_KhMEKE_u>0) call post_data(CS%id_KhMEKE_u, Kh_u, CS%diag)
if (CS%id_KhMEKE_v>0) call post_data(CS%id_KhMEKE_v, Kh_v, CS%diag)
if (CS%id_src>0) call post_data(CS%id_src, src, CS%diag)
if (CS%id_decay>0) call post_data(CS%id_decay, MEKE_decay, CS%diag)
if (CS%id_GM_src>0) call post_data(CS%id_GM_src, MEKE%GM_src, CS%diag)
if (CS%id_mom_src>0) call post_data(CS%id_mom_src, MEKE%mom_src, CS%diag)
+ if (CS%id_GME_snk>0) call post_data(CS%id_GME_snk, MEKE%GME_snk, CS%diag)
if (CS%id_Le>0) call post_data(CS%id_Le, LmixScale, CS%diag)
if (CS%id_gamma_b>0) then
do j=js,je ; do i=is,ie
@@ -695,6 +769,7 @@ subroutine MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_m
EKE = 0.
endif
MEKE%MEKE(i,j) = EKE
+! MEKE%MEKE(i,j) = (US%Z_to_m*G%bathyT(i,j)*SN / (8*CS%cdrag))**2
enddo ; enddo
end subroutine MEKE_equilibrium
@@ -864,7 +939,7 @@ logical function MEKE_init(Time, G, param_file, diag, CS, MEKE, restart_CS)
type(MOM_restart_CS), pointer :: restart_CS !< Restart control structure for MOM_MEKE.
! Local variables
integer :: is, ie, js, je, isd, ied, jsd, jed, nz
- logical :: laplacian, useVarMix, coldStart
+ logical :: laplacian, biharmonic, useVarMix, coldStart
! This include declares and sets the variable "version".
#include "version_variable.h"
character(len=40) :: mdl = "MOM_MEKE" ! This module's name.
@@ -919,10 +994,17 @@ logical function MEKE_init(Time, G, param_file, diag, CS, MEKE, restart_CS)
"into MEKE by the thickness mixing parameterization. "//&
"If MEKE_GMCOEFF is negative, this conversion is not "//&
"used or calculated.", units="nondim", default=-1.0)
+ call get_param(param_file, mdl, "MEKE_GEOMETRIC", CS%MEKE_GEOMETRIC, &
+ "If MEKE_GEOMETRIC is true, uses the GM coefficient formulation "//&
+ "from the GEOMETRIC framework (Marshall et al., 2012).", default=.false.)
call get_param(param_file, mdl, "MEKE_FRCOEFF", CS%MEKE_FrCoeff, &
"The efficiency of the conversion of mean energy into "//&
"MEKE. If MEKE_FRCOEFF is negative, this conversion "//&
"is not used or calculated.", units="nondim", default=-1.0)
+ call get_param(param_file, mdl, "MEKE_GMECOEFF", CS%MEKE_GMECoeff, &
+ "The efficiency of the conversion of MEKE into mean energy "//&
+ "by GME. If MEKE_GMECOEFF is negative, this conversion "//&
+ "is not used or calculated.", units="nondim", default=-1.0)
call get_param(param_file, mdl, "MEKE_BGSRC", CS%MEKE_BGsrc, &
"A background energy source for MEKE.", units="W kg-1", &
default=0.0)
@@ -947,6 +1029,12 @@ logical function MEKE_init(Time, G, param_file, diag, CS, MEKE, restart_CS)
call get_param(param_file, mdl, "MEKE_USCALE", CS%MEKE_Uscale, &
"The background velocity that is combined with MEKE to "//&
"calculate the bottom drag.", units="m s-1", default=0.0)
+ call get_param(param_file, mdl, "MEKE_JANSEN15_DRAG", CS%Jansen15_drag, &
+ "If true, use the bottom drag formulation from Jansen et al. (2015) "//&
+ "to calculate the drag acting on MEKE.", default=.false.)
+ call get_param(param_file, mdl, "MEKE_GM_SRC_ALT", CS%GM_src_alt, &
+ "If true, use the GM energy conversion form S^2*N^2*kappa rather "//&
+ "than the streamfunction for the MEKE GM source term.", default=.false.)
call get_param(param_file, mdl, "MEKE_VISC_DRAG", CS%visc_drag, &
"If true, use the vertvisc_type to calculate the bottom "//&
"drag acting on MEKE.", default=.true.)
@@ -971,10 +1059,16 @@ logical function MEKE_init(Time, G, param_file, diag, CS, MEKE, restart_CS)
"If true, the length scale used by MEKE is the minimum of "//&
"the deformation radius or grid-spacing. Only used if "//&
"MEKE_OLD_LSCALE=True", units="nondim", default=.false.)
- call get_param(param_file, mdl, "MEKE_VISCOSITY_COEFF", CS%viscosity_coeff, &
- "If non-zero, is the scaling coefficient in the expression for "//&
- "viscosity used to parameterize lateral momentum mixing by "//&
- "unresolved eddies represented by MEKE. Can be negative to "//&
+ call get_param(param_file, mdl, "MEKE_VISCOSITY_COEFF_KU", CS%viscosity_coeff_Ku, &
+ "If non-zero, is the scaling coefficient in the expression for"//&
+ "viscosity used to parameterize harmonic lateral momentum mixing by"//&
+ "unresolved eddies represented by MEKE. Can be negative to"//&
+ "represent backscatter from the unresolved eddies.", &
+ units="nondim", default=0.0)
+ call get_param(param_file, mdl, "MEKE_VISCOSITY_COEFF_AU", CS%viscosity_coeff_Au, &
+ "If non-zero, is the scaling coefficient in the expression for"//&
+ "viscosity used to parameterize biharmonic lateral momentum mixing by"//&
+ "unresolved eddies represented by MEKE. Can be negative to"//&
"represent backscatter from the unresolved eddies.", &
units="nondim", default=0.0)
call get_param(param_file, mdl, "MEKE_FIXED_MIXING_LENGTH", CS%Lfixed, &
@@ -1028,8 +1122,13 @@ logical function MEKE_init(Time, G, param_file, diag, CS, MEKE, restart_CS)
"the velocity field to the bottom stress.", units="nondim", &
default=0.003)
call get_param(param_file, mdl, "LAPLACIAN", laplacian, default=.false., do_not_log=.true.)
- if (CS%viscosity_coeff/=0. .and. .not. laplacian) call MOM_error(FATAL, &
- "LAPLACIAN must be true if MEKE_VISCOSITY_COEFF is true.")
+ call get_param(param_file, mdl, "BIHARMONIC", biharmonic, default=.false., do_not_log=.true.)
+
+ if (CS%viscosity_coeff_Ku/=0. .and. .not. laplacian) call MOM_error(FATAL, &
+ "LAPLACIAN must be true if MEKE_VISCOSITY_COEFF_KU is true.")
+
+ if (CS%viscosity_coeff_Au/=0. .and. .not. biharmonic) call MOM_error(FATAL, &
+ "BIHARMONIC must be true if MEKE_VISCOSITY_COEFF_AU is true.")
call get_param(param_file, mdl, "DEBUG", CS%debug, default=.false., do_not_log=.true.)
@@ -1054,10 +1153,18 @@ logical function MEKE_init(Time, G, param_file, diag, CS, MEKE, restart_CS)
call create_group_pass(CS%pass_Kh, MEKE%Kh, G%Domain)
call do_group_pass(CS%pass_Kh, G%Domain)
endif
+ if (associated(MEKE%Kh_diff)) then
+ call create_group_pass(CS%pass_Kh_diff, MEKE%Kh_diff, G%Domain)
+ call do_group_pass(CS%pass_Kh_diff, G%Domain)
+ endif
if (associated(MEKE%Ku)) then
call create_group_pass(CS%pass_Ku, MEKE%Ku, G%Domain)
call do_group_pass(CS%pass_Ku, G%Domain)
endif
+ if (associated(MEKE%Au)) then
+ call create_group_pass(CS%pass_Au, MEKE%Au, G%Domain)
+ call do_group_pass(CS%pass_Au, G%Domain)
+ endif
if (allocated(CS%del2MEKE)) then
call create_group_pass(CS%pass_del2MEKE, CS%del2MEKE, G%Domain)
call do_group_pass(CS%pass_del2MEKE, G%Domain)
@@ -1074,6 +1181,9 @@ logical function MEKE_init(Time, G, param_file, diag, CS, MEKE, restart_CS)
CS%id_Ku = register_diag_field('ocean_model', 'MEKE_KU', diag%axesT1, Time, &
'MEKE derived lateral viscosity', 'm2 s-1')
if (.not. associated(MEKE%Ku)) CS%id_Ku = -1
+ CS%id_Au = register_diag_field('ocean_model', 'MEKE_AU', diag%axesT1, Time, &
+ 'MEKE derived lateral biharmonic viscosity', 'm4 s-1')
+ if (.not. associated(MEKE%Au)) CS%id_Au = -1
CS%id_Ue = register_diag_field('ocean_model', 'MEKE_Ue', diag%axesT1, Time, &
'MEKE derived eddy-velocity scale', 'm s-1')
if (.not. associated(MEKE%MEKE)) CS%id_Ue = -1
@@ -1093,6 +1203,9 @@ logical function MEKE_init(Time, G, param_file, diag, CS, MEKE, restart_CS)
CS%id_mom_src = register_diag_field('ocean_model', 'MEKE_mom_src',diag%axesT1, Time, &
'MEKE energy available from momentum', 'W m-2')
if (.not. associated(MEKE%mom_src)) CS%id_mom_src = -1
+ CS%id_GME_snk = register_diag_field('ocean_model', 'MEKE_GME_snk',diag%axesT1, Time, &
+ 'MEKE energy lost to GME backscatter', 'W m-2')
+ if (.not. associated(MEKE%GME_snk)) CS%id_GME_snk = -1
CS%id_Le = register_diag_field('ocean_model', 'MEKE_Le', diag%axesT1, Time, &
'Eddy mixing length used in the MEKE derived eddy diffusivity', 'm')
CS%id_Lrhines = register_diag_field('ocean_model', 'MEKE_Lrhines', diag%axesT1, Time, &
@@ -1132,7 +1245,8 @@ subroutine MEKE_alloc_register_restart(HI, param_file, MEKE, restart_CS)
type(MOM_restart_CS), pointer :: restart_CS !< Restart control structure for MOM_MEKE.
! Local variables
type(vardesc) :: vd
- real :: MEKE_GMcoeff, MEKE_FrCoeff, MEKE_KHCoeff, MEKE_viscCoeff
+ real :: MEKE_GMcoeff, MEKE_FrCoeff, MEKE_GMECoeff, MEKE_KHCoeff, MEKE_viscCoeff_Ku, MEKE_viscCoeff_Au
+ logical :: Use_KH_in_MEKE
logical :: useMEKE
integer :: isd, ied, jsd, jed
@@ -1142,9 +1256,11 @@ subroutine MEKE_alloc_register_restart(HI, param_file, MEKE, restart_CS)
! Read these parameters to determine what should be in the restarts
MEKE_GMcoeff =-1.; call read_param(param_file,"MEKE_GMCOEFF",MEKE_GMcoeff)
MEKE_FrCoeff =-1.; call read_param(param_file,"MEKE_FRCOEFF",MEKE_FrCoeff)
+ MEKE_GMEcoeff =-1.; call read_param(param_file,"MEKE_GMECOEFF",MEKE_GMEcoeff)
MEKE_KhCoeff =1.; call read_param(param_file,"MEKE_KHCOEFF",MEKE_KhCoeff)
- MEKE_viscCoeff =0.; call read_param(param_file,"MEKE_VISCOSITY_COEFF",MEKE_viscCoeff)
-
+ MEKE_viscCoeff_Ku =0.; call read_param(param_file,"MEKE_VISCOSITY_COEFF_KU",MEKE_viscCoeff_Ku)
+ MEKE_viscCoeff_Au =0.; call read_param(param_file,"MEKE_VISCOSITY_COEFF_AU",MEKE_viscCoeff_Au)
+ Use_KH_in_MEKE = .false.; call read_param(param_file,"USE_KH_IN_MEKE", Use_KH_in_MEKE)
! Allocate control structure
if (associated(MEKE)) then
call MOM_error(WARNING, "MEKE_alloc_register_restart called with an associated "// &
@@ -1164,9 +1280,12 @@ subroutine MEKE_alloc_register_restart(HI, param_file, MEKE, restart_CS)
if (MEKE_GMcoeff>=0.) then
allocate(MEKE%GM_src(isd:ied,jsd:jed)) ; MEKE%GM_src(:,:) = 0.0
endif
- if (MEKE_FrCoeff>=0.) then
+ if (MEKE_FrCoeff>=0. .or. MEKE_GMECoeff>=0.) then
allocate(MEKE%mom_src(isd:ied,jsd:jed)) ; MEKE%mom_src(:,:) = 0.0
endif
+ if (MEKE_GMECoeff>=0.) then
+ allocate(MEKE%GME_snk(isd:ied,jsd:jed)) ; MEKE%GME_snk(:,:) = 0.0
+ endif
if (MEKE_KhCoeff>=0.) then
allocate(MEKE%Kh(isd:ied,jsd:jed)) ; MEKE%Kh(:,:) = 0.0
vd = var_desc("MEKE_Kh", "m2 s-1",hor_grid='h',z_grid='1', &
@@ -1174,12 +1293,25 @@ subroutine MEKE_alloc_register_restart(HI, param_file, MEKE, restart_CS)
call register_restart_field(MEKE%Kh, vd, .false., restart_CS)
endif
allocate(MEKE%Rd_dx_h(isd:ied,jsd:jed)) ; MEKE%Rd_dx_h(:,:) = 0.0
- if (MEKE_viscCoeff/=0.) then
+ if (MEKE_viscCoeff_Ku/=0.) then
allocate(MEKE%Ku(isd:ied,jsd:jed)) ; MEKE%Ku(:,:) = 0.0
- vd = var_desc("MEKE_Ah", "m2 s-1", hor_grid='h', z_grid='1', &
+ vd = var_desc("MEKE_Ku", "m2 s-1", hor_grid='h', z_grid='1', &
longname="Lateral viscosity from Mesoscale Eddy Kinetic Energy")
call register_restart_field(MEKE%Ku, vd, .false., restart_CS)
endif
+ if (Use_Kh_in_MEKE) then
+ allocate(MEKE%Kh_diff(isd:ied,jsd:jed)) ; MEKE%Kh_diff(:,:) = 0.0
+ vd = var_desc("MEKE_Kh_diff", "m2 s-1",hor_grid='h',z_grid='1', &
+ longname="Copy of thickness diffusivity for diffusing MEKE")
+ call register_restart_field(MEKE%Kh_diff, vd, .false., restart_CS)
+ endif
+
+ if (MEKE_viscCoeff_Au/=0.) then
+ allocate(MEKE%Au(isd:ied,jsd:jed)) ; MEKE%Au(:,:) = 0.0
+ vd = var_desc("MEKE_Au", "m4 s-1", hor_grid='h', z_grid='1', &
+ longname="Lateral biharmonic viscosity from Mesoscale Eddy Kinetic Energy")
+ call register_restart_field(MEKE%Au, vd, .false., restart_CS)
+ endif
end subroutine MEKE_alloc_register_restart
@@ -1196,8 +1328,11 @@ subroutine MEKE_end(MEKE, CS)
if (associated(MEKE%MEKE)) deallocate(MEKE%MEKE)
if (associated(MEKE%GM_src)) deallocate(MEKE%GM_src)
if (associated(MEKE%mom_src)) deallocate(MEKE%mom_src)
+ if (associated(MEKE%GME_snk)) deallocate(MEKE%GME_snk)
if (associated(MEKE%Kh)) deallocate(MEKE%Kh)
+ if (associated(MEKE%Kh_diff)) deallocate(MEKE%Kh_diff)
if (associated(MEKE%Ku)) deallocate(MEKE%Ku)
+ if (associated(MEKE%Au)) deallocate(MEKE%Au)
if (allocated(CS%del2MEKE)) deallocate(CS%del2MEKE)
deallocate(MEKE)
@@ -1365,9 +1500,13 @@ end subroutine MEKE_end
!! \subsection section_MEKE_viscosity Viscosity derived from MEKE
!!
!! As for \f$ \kappa_M \f$, the predicted eddy velocity scale can be
-!! used to form an eddy viscosity:
+!! used to form a harmonic eddy viscosity,
+!!
+!! \f[ \kappa_u = \gamma_u \sqrt{ U_e^2 A_\Delta } \f]
+!!
+!! as well as a biharmonic eddy viscosity,
!!
-!! \f[ \kappa_u = \gamma_u \sqrt{ U_e^2 A_\Delta } . \f]
+!! \f[ \kappa_4 = \gamma_4 \sqrt{ U_e^2 A_\Delta^3 } \f]
!!
!! \subsection section_MEKE_limit_case Limit cases for local source-dissipative balance
!!
@@ -1404,7 +1543,8 @@ end subroutine MEKE_end
!! | \f$ \kappa_4 \f$ | MEKE_K4 |
!! | \f$ \gamma_\kappa \f$ | MEKE_KHCOEFF |
!! | \f$ \gamma_M \f$ | MEKE_KHMEKE_FAC |
-!! | \f$ \gamma_u \f$ | MEKE_VISCOSITY_COEFF |
+!! | \f$ \gamma_u \f$ | MEKE_VISCOSITY_COEFF_KU |
+!! | \f$ \gamma_4 \f$ | MEKE_VISCOSITY_COEFF_AU |
!! | \f$ \gamma_{min}^2 \f$| MEKE_MIN_GAMMA2 |
!! | \f$ \alpha_d \f$ | MEKE_ALPHA_DEFORM |
!! | \f$ \alpha_f \f$ | MEKE_ALPHA_FRICT |
diff --git a/src/parameterizations/lateral/MOM_MEKE_types.F90 b/src/parameterizations/lateral/MOM_MEKE_types.F90
index 22ed34c6c2..95106f1fdb 100644
--- a/src/parameterizations/lateral/MOM_MEKE_types.F90
+++ b/src/parameterizations/lateral/MOM_MEKE_types.F90
@@ -11,17 +11,21 @@ module MOM_MEKE_types
MEKE => NULL(), & !< Vertically averaged eddy kinetic energy [m2 s-2].
GM_src => NULL(), & !< MEKE source due to thickness mixing (GM) [W m-2].
mom_src => NULL(),& !< MEKE source from lateral friction in the momentum equations [W m-2].
- Kh => NULL(), & !< The MEKE-derived lateral mixing coefficient [m2 s-1.
+ GME_snk => NULL(),& !< MEKE sink from GME backscatter in the momentum equations [W m-2].
+ Kh => NULL(), & !< The MEKE-derived lateral mixing coefficient [m2 s-1].
+ Kh_diff => NULL(), & !< Uses the non-MEKE-derived thickness diffusion coefficient to diffuse MEKE [m2 s-1].
Rd_dx_h => NULL() !< The deformation radius compared with the grid spacing [nondim].
!! Rd_dx_h is copied from VarMix_CS.
real, dimension(:,:), pointer :: Ku => NULL() !< The MEKE-derived lateral viscosity coefficient [m2 s-1].
!! This viscosity can be negative when representing backscatter
!! from unresolved eddies (see Jansen and Held, 2014).
+ real, dimension(:,:), pointer :: Au => NULL() !< The MEKE-derived lateral biharmonic viscosity coefficient [m4 s-1].
! Parameters
real :: KhTh_fac = 1.0 !< Multiplier to map Kh(MEKE) to KhTh [nondim]
real :: KhTr_fac = 1.0 !< Multiplier to map Kh(MEKE) to KhTr [nondim].
real :: backscatter_Ro_pow = 0.0 !< Power in Rossby number function for backscatter.
real :: backscatter_Ro_c = 0.0 !< Coefficient in Rossby number function for backscatter.
+
end type MEKE_type
end module MOM_MEKE_types
diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90
index c43f6744d0..919ad02820 100644
--- a/src/parameterizations/lateral/MOM_hor_visc.F90
+++ b/src/parameterizations/lateral/MOM_hor_visc.F90
@@ -5,12 +5,14 @@ module MOM_hor_visc
use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
use MOM_diag_mediator, only : diag_ctrl, time_type
-use MOM_domains, only : pass_var
-use MOM_error_handler, only : MOM_error, FATAL, WARNING
+use MOM_domains, only : pass_var, CORNER, pass_vector
+use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe
use MOM_file_parser, only : get_param, log_version, param_file_type
use MOM_grid, only : ocean_grid_type
+use MOM_lateral_mixing_coeffs, only : VarMix_CS, calc_QG_Leith_viscosity
+use MOM_barotropic, only : barotropic_CS, barotropic_get_tav
+use MOM_thickness_diffuse, only : thickness_diffuse_CS, thickness_diffuse_get_KH
use MOM_io, only : MOM_read_data, slasher
-use MOM_lateral_mixing_coeffs, only : VarMix_CS
use MOM_MEKE_types, only : MEKE_type
use MOM_open_boundary, only : ocean_OBC_type, OBC_DIRECTION_E, OBC_DIRECTION_W
use MOM_open_boundary, only : OBC_DIRECTION_N, OBC_DIRECTION_S, OBC_NONE
@@ -54,8 +56,11 @@ module MOM_hor_visc
!! viscosity. KH is the background value.
logical :: Modified_Leith !< If true, use extra component of Leith viscosity
!! to damp divergent flow. To use, still set Leith_Kh=.TRUE.
+ logical :: use_beta_in_Leith !< If true, includes the beta term in the Leith viscosity
logical :: Leith_Ah !< If true, use a biharmonic form of 2D Leith
!! nonlinear eddy viscosity. AH is the background.
+ logical :: use_QG_Leith_visc !< If true, use QG Leith nonlinear eddy viscosity.
+ !! KH is the background value.
logical :: bound_Coriolis !< If true & SMAGORINSKY_AH is used, the biharmonic
!! viscosity is modified to include a term that
!! scales quadratically with the velocity shears.
@@ -73,6 +78,8 @@ module MOM_hor_visc
!! of state. This is set depending on ANISOTROPIC_MODE.
logical :: res_scale_MEKE !< If true, the viscosity contribution from MEKE is scaled by
!! the resolution function.
+ logical :: use_GME !< If true, use GME backscatter scheme.
+
real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: Kh_bg_xx
!< The background Laplacian viscosity at h points [m2 s-1].
!! The actual viscosity may be the larger of this
@@ -85,7 +92,7 @@ module MOM_hor_visc
!< The background biharmonic viscosity at h points [m4 s-1].
!! The actual viscosity may be the larger of this
!! viscosity and the Smagorinsky and Leith viscosities.
- real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: Biharm_Const2_xx
+ real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: Biharm5_const2_xx
!< A constant relating the biharmonic viscosity to the
!! square of the velocity shear [m4 s]. This value is
!! set to be the magnitude of the Coriolis terms once the
@@ -107,7 +114,7 @@ module MOM_hor_visc
!< The background biharmonic viscosity at q points [m4 s-1].
!! The actual viscosity may be the larger of this
!! viscosity and the Smagorinsky and Leith viscosities.
- real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: Biharm_Const2_xy
+ real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: Biharm5_const2_xy
!< A constant relating the biharmonic viscosity to the
!! square of the velocity shear [m4 s]. This value is
!! set to be the magnitude of the Coriolis terms once the
@@ -141,16 +148,18 @@ module MOM_hor_visc
! The following variables are precalculated time-invariant combinations of
! parameters and metric terms.
real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: &
- Laplac_Const_xx, & !< Laplacian metric-dependent constants [nondim]
- Biharm_Const_xx, & !< Biharmonic metric-dependent constants [nondim]
- Laplac3_Const_xx, & !< Laplacian metric-dependent constants [nondim]
- Biharm5_Const_xx !< Biharmonic metric-dependent constants [nondim]
+ Laplac2_const_xx, & !< Laplacian metric-dependent constants [nondim]
+ Biharm5_const_xx, & !< Biharmonic metric-dependent constants [nondim]
+ Laplac3_const_xx, & !< Laplacian metric-dependent constants [nondim]
+ Biharm_const_xx, & !< Biharmonic metric-dependent constants [nondim]
+ Biharm_const2_xx !< Biharmonic metric-dependent constants [nondim]
real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: &
- Laplac_Const_xy, & !< Laplacian metric-dependent constants [nondim]
- Biharm_Const_xy, & !< Biharmonic metric-dependent constants [nondim]
- Laplac3_Const_xy, & !< Laplacian metric-dependent constants [nondim]
- Biharm5_Const_xy !< Biharmonic metric-dependent constants [nondim]
+ Laplac2_const_xy, & !< Laplacian metric-dependent constants [nondim]
+ Biharm5_const_xy, & !< Biharmonic metric-dependent constants [nondim]
+ Laplac3_const_xy, & !< Laplacian metric-dependent constants [nondim]
+ Biharm_const_xy, & !< Biharmonic metric-dependent constants [nondim]
+ Biharm_const2_xy !< Biharmonic metric-dependent constants [nondim]
type(diag_ctrl), pointer :: diag => NULL() !< structure to regulate diagnostics
@@ -159,9 +168,14 @@ module MOM_hor_visc
integer :: id_diffu = -1, id_diffv = -1
integer :: id_Ah_h = -1, id_Ah_q = -1
integer :: id_Kh_h = -1, id_Kh_q = -1
+ integer :: id_GME_coeff_h = -1, id_GME_coeff_q = -1
+ integer :: id_vort_xy_q = -1, id_div_xx_h = -1
integer :: id_FrictWork = -1, id_FrictWorkIntz = -1
+ integer :: id_FrictWorkMax = -1
+ integer :: id_FrictWork_diss = -1, id_FrictWork_GME = -1
!!@}
+
end type hor_visc_CS
contains
@@ -178,7 +192,8 @@ module MOM_hor_visc
!! u[is-2:ie+2,js-2:je+2]
!! v[is-2:ie+2,js-2:je+2]
!! h[is-1:ie+1,js-1:je+1]
-subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, CS, OBC)
+subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, &
+ CS, OBC, BT)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
@@ -186,7 +201,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
intent(in) :: v !< The meridional velocity [m s-1].
real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
- intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
+ intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2].
real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
intent(out) :: diffu !< Zonal acceleration due to convergence of
!! along-coordinate stress tensor [m s-2]
@@ -198,48 +213,92 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
type(VarMix_CS), pointer :: VarMix !< Pointer to a structure with fields that
!! specify the spatially variable viscosities
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
- type(hor_visc_CS), pointer :: CS !< Pontrol structure returned by a previous
+ type(hor_visc_CS), pointer :: CS !< Control structure returned by a previous
!! call to hor_visc_init.
type(ocean_OBC_type), optional, pointer :: OBC !< Pointer to an open boundary condition type
+ type(barotropic_CS), optional, pointer :: BT !< Pointer to a structure containing
+ !! barotropic velocities.
+
! Local variables
real, dimension(SZIB_(G),SZJ_(G)) :: &
- u0, & ! Laplacian of u [m-1 s-1]
- h_u ! Thickness interpolated to u points [H ~> m or kg m-2].
+ u0, & ! Laplacian of u [m-1 s-1]
+ h_u, & ! Thickness interpolated to u points [H ~> m or kg m-2].
+ vort_xy_dy, & ! y-derivative of vertical vorticity (d/dy(dv/dx - du/dy)) [m-1 s-1]
+ div_xx_dx, & ! x-derivative of horizontal divergence (d/dx(du/dx + dv/dy)) [m-1 s-1]
+ ubtav ! zonal barotropic vel. ave. over baroclinic time-step [m s-1]
real, dimension(SZI_(G),SZJB_(G)) :: &
- v0, & ! Laplacian of v [m-1 s-1]
- h_v ! Thickness interpolated to v points [H ~> m or kg m-2].
-
+ v0, & ! Laplacian of v [m-1 s-1]
+ h_v, & ! Thickness interpolated to v points [H ~> m or kg m-2].
+ vort_xy_dx, & ! x-derivative of vertical vorticity (d/dx(dv/dx - du/dy)) [m-1 s-1]
+ div_xx_dy, & ! y-derivative of horizontal divergence (d/dy(du/dx + dv/dy)) [m-1 s-1]
+ vbtav ! meridional barotropic vel. ave. over baroclinic time-step [m s-1]
real, dimension(SZI_(G),SZJ_(G)) :: &
+ dudx_bt, dvdy_bt, & ! components in the barotropic horizontal tension [s-1]
+ div_xx, & ! Estimate of horizontal divergence at h-points [s-1]
sh_xx, & ! horizontal tension (du/dx - dv/dy) including metric terms [s-1]
+ sh_xx_bt, & ! barotropic horizontal tension (du/dx - dv/dy) including metric terms [s-1]
str_xx,& ! str_xx is the diagonal term in the stress tensor [H m2 s-2 ~> m3 s-2 or kg s-2]
+ str_xx_GME,& ! smoothed diagonal term in the stress tensor from GME [H m2 s-2]
bhstr_xx,& ! A copy of str_xx that only contains the biharmonic contribution [H m2 s-2 ~> m3 s-2 or kg s-2]
- div_xx, & ! horizontal divergence (du/dx + dv/dy) including metric terms [s-1]
- FrictWorkIntz ! depth integrated energy dissipated by lateral friction [W m-2]
+ FrictWorkIntz, & ! depth integrated energy dissipated by lateral friction [W m-2]
+ Leith_Kh_h, & ! Leith Laplacian viscosity at h-points [m2 s-1]
+ Leith_Ah_h, & ! Leith bi-harmonic viscosity at h-points [m4 s-1]
+ beta_h, & ! Gradient of planetary vorticity at h-points [m-1 s-1]
+ grad_vort_mag_h, & ! Magnitude of vorticity gradient at h-points [m-1 s-1]
+ grad_vort_mag_h_2d, & ! Magnitude of 2d vorticity gradient at h-points [m-1 s-1]
+ grad_div_mag_h, & ! Magnitude of divergence gradient at h-points [m-1 s-1]
+ dudx, dvdy, & ! components in the horizontal tension [s-1]
+ grad_vel_mag_h, & ! Magnitude of the velocity gradient tensor squared at h-points [s-2]
+ grad_vel_mag_bt_h, & ! Magnitude of the barotropic velocity gradient tensor squared at h-points [s-2]
+ grad_d2vel_mag_h, & ! Magnitude of the Laplacian of the velocity vector, squared [m-2 s-2]
+ max_diss_rate_bt, & ! maximum possible energy dissipated by barotropic lateral friction [m2 s-3]
+ boundary_mask ! A mask that zeroes out cells with at least one land edge
real, dimension(SZIB_(G),SZJB_(G)) :: &
dvdx, dudy, & ! components in the shearing strain [s-1]
+ dvdx_bt, dudy_bt, & ! components in the barotropic shearing strain [s-1]
sh_xy, & ! horizontal shearing strain (du/dy + dv/dx) including metric terms [s-1]
+ sh_xy_bt, & ! barotropic horizontal shearing strain (du/dy + dv/dx) inc. metric terms [s-1]
str_xy, & ! str_xy is the cross term in the stress tensor [H m2 s-2 ~> m3 s-2 or kg s-2]
+ str_xy_GME, & ! smoothed cross term in the stress tensor from GME [H m2 s-2]
bhstr_xy, & ! A copy of str_xy that only contains the biharmonic contribution [H m2 s-2 ~> m3 s-2 or kg s-2]
- vort_xy ! vertical vorticity (dv/dx - du/dy) including metric terms [s-1]
-
- real, dimension(SZI_(G),SZJB_(G)) :: &
- vort_xy_dx, & ! x-derivative of vertical vorticity (d/dx(dv/dx - du/dy)) including metric terms [m-1 s-1]
- div_xx_dy ! y-derivative of horizontal divergence (d/dy(du/dx + dv/dy)) including metric terms [m-1 s-1]
-
- real, dimension(SZIB_(G),SZJ_(G)) :: &
- vort_xy_dy, & ! y-derivative of vertical vorticity (d/dy(dv/dx - du/dy)) including metric terms [m-1 s-1]
- div_xx_dx ! x-derivative of horizontal divergence (d/dx(du/dx + dv/dy)) including metric terms [m-1 s-1]
+ vort_xy, & ! Vertical vorticity (dv/dx - du/dy) including metric terms [s-1]
+ Leith_Kh_q, & ! Leith Laplacian viscosity at q-points [m2 s-1]
+ Leith_Ah_q, & ! Leith bi-harmonic viscosity at q-points [m4 s-1]
+ beta_q, & ! Gradient of planetary vorticity at q-points [m-1 s-1]
+ grad_vort_mag_q, & ! Magnitude of vorticity gradient at q-points [m-1 s-1]
+ grad_vort_mag_q_2d, & ! Magnitude of 2d vorticity gradient at q-points [m-1 s-1]
+ grad_div_mag_q, & ! Magnitude of divergence gradient at q-points [m-1 s-1]
+ grad_vel_mag_q, & ! Magnitude of the velocity gradient tensor squared at q-points [s-2]
+ hq, & ! harmonic mean of the harmonic means of the u- & v point thicknesses, in H; This form guarantees that hq/hu < 4.
+ grad_vel_mag_bt_q ! Magnitude of the barotropic velocity gradient tensor squared at q-points [s-2]
real, dimension(SZIB_(G),SZJB_(G),SZK_(G)) :: &
- Ah_q, & ! biharmonic viscosity at corner points [m4 s-1]
- Kh_q ! Laplacian viscosity at corner points [m2 s-1]
-
+ Ah_q, & ! biharmonic viscosity at corner points [m4 s-1]
+ Kh_q, & ! Laplacian viscosity at corner points [m2 s-1]
+ vort_xy_q, & ! vertical vorticity at corner points [s-1]
+ GME_coeff_q !< GME coeff. at q-points [m2 s-1]
+
+ real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1) :: &
+ KH_u_GME !< interface height diffusivities in u-columns [m2 s-1]
+ real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1) :: &
+ KH_v_GME !< interface height diffusivities in v-columns [m2 s-1]
real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
Ah_h, & ! biharmonic viscosity at thickness points [m4 s-1]
Kh_h, & ! Laplacian viscosity at thickness points [m2 s-1]
- FrictWork ! energy dissipated by lateral friction [W m-2]
-
+ diss_rate, & ! MKE dissipated by parameterized shear production [m2 s-3]
+ max_diss_rate, & ! maximum possible energy dissipated by lateral friction [m2 s-3]
+ target_diss_rate_GME, & ! the maximum theoretical dissipation plus the amount spuriously dissipated
+ ! by friction [m2 s-3]
+ FrictWork, & ! work done by MKE dissipation mechanisms [W m-2]
+ FrictWork_diss, & ! negative definite work done by MKE dissipation mechanisms [W m-2]
+ FrictWorkMax, & ! maximum possible work done by MKE dissipation mechanisms [W m-2]
+ FrictWork_GME, & ! work done by GME [W m-2]
+ div_xx_h ! horizontal divergence [s-1]
+ !real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &
+ real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
+ KH_t_GME, & !< interface height diffusivities in t-columns [m2 s-1]
+ GME_coeff_h !< GME coeff. at h-points [m2 s-1]
real :: Ah ! biharmonic viscosity [m4 s-1]
real :: Kh ! Laplacian viscosity [m2 s-1]
real :: AhSm ! Smagorinsky biharmonic viscosity [m4 s-1]
@@ -250,13 +309,13 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
! viscosity. Here set equal to nondimensional Laplacian Leith constant.
! This is set equal to zero if modified Leith is not used.
real :: Shear_mag ! magnitude of the shear [s-1]
- real :: Vort_mag ! magnitude of the vorticity [s-1]
+ real :: vert_vort_mag ! magnitude of the vertical vorticity gradient [m-1 s-1]
real :: h2uq, h2vq ! temporary variables [H2 ~> m2 or kg2 m-4].
real :: hu, hv ! Thicknesses interpolated by arithmetic means to corner
! points; these are first interpolated to u or v velocity
! points where masks are applied [H ~> m or kg m-2].
- real :: hq ! harmonic mean of the harmonic means of the u- & v- poing thicknesses,
- ! [H ~> m or kg m-2]; This form guarantees that hq/hu < 4.
+! real :: hq ! harmonic mean of the harmonic means of the u- & v-
+! ! point thicknesses, in H; This form guarantees that hq/hu < 4.
real :: h_neglect ! thickness so small it can be lost in roundoff and so neglected [H ~> m or kg m-2]
real :: h_neglect3 ! h_neglect^3 [H3 ~> m3 or kg3 m-6]
real :: hrat_min ! minimum thicknesses at the 4 neighboring
@@ -270,19 +329,30 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
real :: FatH ! abs(f) at h-point for MEKE source term [s-1]
real :: local_strain ! Local variable for interpolating computed strain rates [s-1].
real :: meke_res_fn ! A copy of the resolution scaling factor if being applied to MEKE. Otherwise =1.
-
+ real :: GME_coeff ! The GME (negative) viscosity coefficient [m2 s-1]
+ real :: GME_coeff_limiter ! Maximum permitted value of the GME coefficient [m2 s-1]
+ real :: FWfrac ! Fraction of maximum theoretical energy transfer to use when scaling GME coefficient
+ real :: DY_dxBu, DX_dyBu
+ real :: H0 ! Depth used to scale down GME coefficient in shallow areas [m]
logical :: rescale_Kh, legacy_bound
logical :: find_FrictWork
logical :: apply_OBC = .false.
logical :: use_MEKE_Ku
+ logical :: use_MEKE_Au
integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
integer :: i, j, k, n
-
+ real :: inv_PI3, inv_PI2, inv_PI5
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
h_neglect = GV%H_subroundoff
h_neglect3 = h_neglect**3
+ inv_PI3 = 1.0/((4.0*atan(1.0))**3)
+ inv_PI2 = 1.0/((4.0*atan(1.0))**2)
+ inv_PI5 = inv_PI3 * inv_PI2
+
+ Ah_h(:,:,:) = 0.0
+ Kh_h(:,:,:) = 0.0
if (present(OBC)) then ; if (associated(OBC)) then ; if (OBC%OBC_pe) then
apply_OBC = OBC%Flather_u_BCs_exist_globally .or. OBC%Flather_v_BCs_exist_globally
@@ -310,23 +380,99 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
legacy_bound = (CS%Smagorinsky_Kh .or. CS%Leith_Kh) .and. &
(CS%bound_Kh .and. .not.CS%better_bound_Kh)
- ! Coefficient for modified Leith
- if (CS%Modified_Leith) then
- mod_Leith = 1.0
- else
- mod_Leith = 0.0
- endif
-
! Toggle whether to use a Laplacian viscosity derived from MEKE
use_MEKE_Ku = associated(MEKE%Ku)
+ use_MEKE_Au = associated(MEKE%Au)
+
+ do j=js,je ; do i=is,ie
+ boundary_mask(i,j) = (G%mask2dCu(I,j) * G%mask2dCv(i,J) * G%mask2dCu(I-1,j) * G%mask2dCv(i,J-1))
+ enddo ; enddo
+
+ if (CS%use_GME) then
+ ! GME tapers off above this depth
+ H0 = 1000.0
+ FWfrac = 1.0
+ GME_coeff_limiter = 1e7
+
+ ! initialize diag. array with zeros
+ GME_coeff_h(:,:,:) = 0.0
+ GME_coeff_q(:,:,:) = 0.0
+ str_xx_GME(:,:) = 0.0
+ str_xy_GME(:,:) = 0.0
+
+! call pass_var(boundary_mask, G%Domain, complete=.true.)
+
+ ! Get barotropic velocities and their gradients
+ call barotropic_get_tav(BT, ubtav, vbtav, G)
+ call pass_vector(ubtav, vbtav, G%Domain)
+
+ do j=js,je ; do i=is,ie
+ dudx_bt(i,j) = CS%DY_dxT(i,j)*(G%IdyCu(I,j) * ubtav(I,j) - &
+ G%IdyCu(I-1,j) * ubtav(I-1,j))
+ dvdy_bt(i,j) = CS%DX_dyT(i,j)*(G%IdxCv(i,J) * vbtav(i,J) - &
+ G%IdxCv(i,J-1) * vbtav(i,J-1))
+ enddo; enddo
+
+ call pass_var(dudx_bt, G%Domain, complete=.true.)
+ call pass_var(dvdy_bt, G%Domain, complete=.true.)
+
+ do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2
+ sh_xx_bt(i,j) = dudx_bt(i,j) - dvdy_bt(i,j)
+ enddo ; enddo
+
+ ! Components for the barotropic shearing strain
+ do J=js-2,Jeq+1 ; do I=is-2,Ieq+1
+ dvdx_bt(I,J) = CS%DY_dxBu(I,J)*(vbtav(i+1,J)*G%IdyCv(i+1,J) &
+ - vbtav(i,J)*G%IdyCv(i,J))
+ dudy_bt(I,J) = CS%DX_dyBu(I,J)*(ubtav(I,j+1)*G%IdxCu(I,j+1) &
+ - ubtav(I,j)*G%IdxCu(I,j))
+ enddo ; enddo
+
+ call pass_var(dvdx_bt, G%Domain, position=CORNER, complete=.true.)
+ call pass_var(dudy_bt, G%Domain, position=CORNER, complete=.true.)
+
+ if (CS%no_slip) then
+ do J=js-2,Jeq+1 ; do I=is-2,Ieq+1
+ sh_xy_bt(I,J) = (2.0-G%mask2dBu(I,J)) * ( dvdx_bt(I,J) + dudy_bt(I,J) )
+ enddo ; enddo
+ else
+ do J=js-2,Jeq+1 ; do I=is-2,Ieq+1
+ sh_xy_bt(I,J) = G%mask2dBu(I,J) * ( dvdx_bt(I,J) + dudy_bt(I,J) )
+ enddo ; enddo
+ endif
+
+ ! Get thickness diffusivity for use in GME
+! call thickness_diffuse_get_KH(thickness_diffuse, KH_u_GME, KH_v_GME, G)
+
+ do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2
+ grad_vel_mag_bt_h(i,j) = boundary_mask(i,j) * (dudx_bt(i,j)**2 + dvdy_bt(i,j)**2 + &
+ (0.25*(dvdx_bt(I,J)+dvdx_bt(I-1,J)+dvdx_bt(I,J-1)+dvdx_bt(I-1,J-1)) )**2 + &
+ (0.25*(dudy_bt(I,J)+dudy_bt(I-1,J)+dudy_bt(I,J-1)+dudy_bt(I-1,J-1)) )**2)
+ enddo ; enddo
+
+ if (associated(MEKE)) then ; if (associated(MEKE%mom_src)) then
+ do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2
+ max_diss_rate_bt(i,j) = 2.0*MEKE%MEKE(i,j) * grad_vel_mag_bt_h(i,j)
+ enddo ; enddo
+ endif ; endif
+
+ do J=js-2,Jeq+1 ; do I=is-2,Ieq+1
+ grad_vel_mag_bt_q(I,J) = boundary_mask(i,j) * (dvdx_bt(i,j)**2 + dudy_bt(i,j)**2 + &
+ (0.25*(dudx_bt(i,j)+dudx_bt(i+1,j)+dudx_bt(i,j+1)+dudx_bt(i+1,j+1)))**2 + &
+ (0.25*(dvdy_bt(i,j)+dvdy_bt(i+1,j)+dvdy_bt(i,j+1)+dvdy_bt(i+1,j+1)) )**2)
+ enddo ; enddo
+
+ endif ! use_GME
!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,CS,G,GV,u,v,is,js,ie,je,h, &
!$OMP rescale_Kh,VarMix,h_neglect,h_neglect3, &
!$OMP Kh_h,Ah_h,Kh_q,Ah_q,diffu,apply_OBC,OBC,diffv, &
- !$OMP find_FrictWork,FrictWork,use_MEKE_Ku,MEKE, &
- !$OMP mod_Leith, legacy_bound) &
- !$OMP private(u0, v0, sh_xx, str_xx, visc_bound_rem, &
+ !$OMP find_FrictWork,FrictWork,use_MEKE_Ku, &
+ !$OMP use_MEKE_Au, MEKE, hq, &
+ !$OMP mod_Leith, legacy_bound, div_xx_h, vort_xy_q) &
+ !$OMP private(u0, v0, sh_xx, str_xx, visc_bound_rem, &
!$OMP sh_xy, str_xy, Ah, Kh, AhSm, KhSm, dvdx, dudy, &
+ !$OMP sh_xx_bt, sh_xy_bt, dvdx_bt, dudy_bt, &
!$OMP bhstr_xx, bhstr_xy,FatH,RoScl, hu, hv, h_u, h_v, &
!$OMP vort_xy,vort_xy_dx,vort_xy_dy,Vort_mag,AhLth,KhLth, &
!$OMP div_xx, div_xx_dx, div_xx_dy, local_strain, &
@@ -340,10 +486,11 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
! Calculate horizontal tension
do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2
- sh_xx(i,j) = (CS%DY_dxT(i,j)*(G%IdyCu(I,j) * u(I,j,k) - &
- G%IdyCu(I-1,j) * u(I-1,j,k)) - &
- CS%DX_dyT(i,j)*(G%IdxCv(i,J) * v(i,J,k) - &
- G%IdxCv(i,J-1)*v(i,J-1,k)))
+ dudx(i,j) = CS%DY_dxT(i,j)*(G%IdyCu(I,j) * u(I,j,k) - &
+ G%IdyCu(I-1,j) * u(I-1,j,k))
+ dvdy(i,j) = CS%DX_dyT(i,j)*(G%IdxCv(i,J) * v(i,J,k) - &
+ G%IdxCv(i,J-1) * v(i,J-1,k))
+ sh_xx(i,j) = dudx(i,j) - dvdy(i,j)
enddo ; enddo
! Components for the shearing strain
@@ -480,18 +627,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
endif
enddo ; endif
- ! Calculate horizontal divergence (not from continuity) if needed.
- ! h_u and h_v include modifications at OBCs from above.
- if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then
- do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2
- div_xx(i,j) = ((G%dyCu(I ,j) * u(I ,j,k) * h_u(I ,j) - &
- G%dyCu(I-1,j) * u(I-1,j,k) * h_u(I-1,j) ) + &
- (G%dxCv(i,J ) * v(i,J ,k) * h_v(i,J ) - &
- G%dxCv(i,J-1) * v(i,J-1,k) * h_v(i,J-1) ) )*G%IareaT(i,j)/ &
- (h(i,j,k) + h_neglect)
- enddo ; enddo
- endif
-
! Shearing strain (including no-slip boundary conditions at the 2-D land-sea mask).
! dudy and dvdx include modifications at OBCs from above.
if (CS%no_slip) then
@@ -504,38 +639,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
enddo ; enddo
endif
- if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then
- ! Calculate relative vorticity (including no-slip boundary conditions at the 2-D land-sea mask).
- ! dudy and dvdx include modifications at OBCs from above.
- if (CS%no_slip) then
- do J=js-2,Jeq+1 ; do I=is-2,Ieq+1
- vort_xy(I,J) = (2.0-G%mask2dBu(I,J)) * ( dvdx(I,J) - dudy(I,J) )
- enddo ; enddo
- else
- do J=js-2,Jeq+1 ; do I=is-2,Ieq+1
- vort_xy(I,J) = G%mask2dBu(I,J) * ( dvdx(I,J) - dudy(I,J) )
- enddo ; enddo
- endif
-
- ! Vorticity gradient
- do J=js-2,Jeq+1 ; do I=is-1,Ieq+1
- 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))
- enddo ; enddo
-
- do J=js-1,Jeq+1 ; do I=is-2,Ieq+1
- 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))
- enddo ; enddo
-
- ! Divergence gradient
- do j=js-1,Jeq+1 ; do I=Isq-1,Ieq+1
- div_xx_dx(I,j) = G%IdxCu(I,j)*(div_xx(i+1,j) - div_xx(i,j))
- enddo ; enddo
-
- do J=Jsq-1,Jeq+1 ; do i=is-1,Ieq+1
- div_xx_dy(i,J) = G%IdyCv(i,J)*(div_xx(i,j+1) - div_xx(i,j))
- enddo ; enddo
- endif
-
! Evaluate u0 = x.Div(Grad u) and v0 = y.Div( Grad u)
if (CS%biharmonic) then
do j=js-1,Jeq+1 ; do I=Isq-1,Ieq+1
@@ -562,20 +665,157 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
endif; endif
endif
+ if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then
+
+ ! Components for the vertical vorticity
+ ! Note this a simple re-calculation of shearing components using the same discretization.
+ ! We will consider using a circulation based calculation of vorticity later.
+ ! Also note this will need OBC boundary conditions re-applied...
+ do J=js-2,Jeq+1 ; do I=is-2,Ieq+1
+ DY_dxBu = G%dyBu(I,J) * G%IdxBu(I,J)
+ dvdx(I,J) = DY_dxBu * (v(i+1,J,k) * G%IdyCv(i+1,J) - v(i,J,k) * G%IdyCv(i,J))
+ DX_dyBu = G%dxBu(I,J) * G%IdyBu(I,J)
+ dudy(I,J) = DX_dyBu * (u(I,j+1,k) * G%IdxCu(I,j+1) - u(I,j,k) * G%IdxCu(I,j))
+ enddo ; enddo
+
+ ! Vorticity
+ if (CS%no_slip) then
+ do J=js-2,Jeq+1 ; do I=is-2,Ieq+1
+ vort_xy(I,J) = (2.0-G%mask2dBu(I,J)) * ( dvdx(I,J) - dudy(I,J) )
+ dudy(I,J) = (2.0-G%mask2dBu(I,J)) * dudy(I,J)
+ dvdx(I,J) = (2.0-G%mask2dBu(I,J)) * dvdx(I,J)
+ enddo ; enddo
+ else
+ do J=js-2,Jeq+1 ; do I=is-2,Ieq+1
+ vort_xy(I,J) = G%mask2dBu(I,J) * ( dvdx(I,J) - dudy(I,J) )
+ dudy(I,J) = G%mask2dBu(I,J) * dudy(I,J)
+ dvdx(I,J) = G%mask2dBu(I,J) * dvdx(I,J)
+ enddo ; enddo
+ endif
+
+ call pass_var(vort_xy, G%Domain, position=CORNER, complete=.true.)
+
+ ! Vorticity gradient
+ do J=js-2,Jeq+1 ; do i=is-1,Ieq+1
+ DY_dxBu = G%dyBu(I,J) * G%IdxBu(I,J)
+ vort_xy_dx(i,J) = DY_dxBu * (vort_xy(I,J) * G%IdyCu(I,j) - vort_xy(I-1,J) * G%IdyCu(I-1,j))
+ enddo ; enddo
+
+ do j=js-1,Jeq+1 ; do I=is-2,Ieq+1
+ DX_dyBu = G%dxBu(I,J) * G%IdyBu(I,J)
+ vort_xy_dy(I,j) = DX_dyBu * (vort_xy(I,J) * G%IdxCv(i,J) - vort_xy(I,J-1) * G%IdxCv(i,J-1))
+ enddo ; enddo
+
+ call pass_vector(vort_xy_dy, vort_xy_dx, G%Domain)
+
+ if (CS%modified_Leith) then
+ ! Divergence
+ do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2
+ div_xx(i,j) = 0.5*((G%dyCu(I,j) * u(I,j,k) * (h(i+1,j,k)+h(i,j,k)) - &
+ G%dyCu(I-1,j) * u(I-1,j,k) * (h(i-1,j,k)+h(i,j,k)) ) + &
+ (G%dxCv(i,J) * v(i,J,k) * (h(i,j,k)+h(i,j+1,k)) - &
+ G%dxCv(i,J-1)*v(i,J-1,k)*(h(i,j,k)+h(i,j-1,k))))*G%IareaT(i,j)/ &
+ (h(i,j,k) + GV%H_subroundoff)
+ enddo ; enddo
+
+ call pass_var(div_xx, G%Domain, complete=.true.)
+
+ ! Divergence gradient
+ do j=Jsq-1,Jeq+2 ; do I=is-2,Ieq+1
+ div_xx_dx(I,j) = G%IdxCu(I,j)*(div_xx(i+1,j) - div_xx(i,j))
+ enddo ; enddo
+ do J=js-2,Jeq+1 ; do i=Isq-1,Ieq+2
+ div_xx_dy(i,J) = G%IdyCv(i,J)*(div_xx(i,j+1) - div_xx(i,j))
+ enddo ; enddo
+
+ call pass_vector(div_xx_dx, div_xx_dy, G%Domain)
+
+ ! Magnitude of divergence gradient
+ do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2
+ grad_div_mag_h(i,j) =sqrt((0.5*(div_xx_dx(I,j) + div_xx_dx(I-1,j)))**2 + &
+ (0.5 * (div_xx_dy(i,J) + div_xx_dy(i,J-1)))**2)
+ enddo ; enddo
+ do J=js-2,Jeq+1 ; do I=is-2,Ieq+1
+ grad_div_mag_q(I,J) =sqrt((0.5*(div_xx_dx(I,j) + div_xx_dx(I,j+1)))**2 + &
+ (0.5 * (div_xx_dy(i,J) + div_xx_dy(i+1,J)))**2)
+ enddo ; enddo
+
+ else
+
+ do j=Jsq-1,Jeq+2 ; do I=is-2,Ieq+1
+ div_xx_dx(I,j) = 0.0
+ enddo ; enddo
+ do J=js-2,Jeq+1 ; do i=Isq-1,Ieq+2
+ div_xx_dy(i,J) = 0.0
+ enddo ; enddo
+ do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2
+ grad_div_mag_h(i,j) = 0.0
+ enddo ; enddo
+ do J=js-2,Jeq+1 ; do I=is-2,Ieq+1
+ grad_div_mag_q(I,J) = 0.0
+ enddo ; enddo
+
+ endif ! CS%modified_Leith
+
+ ! Add in beta for the Leith viscosity
+ if (CS%use_beta_in_Leith) then
+ do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2
+ beta_h(i,j) = sqrt( G%dF_dx(i,j)**2 + G%dF_dy(i,j)**2 )
+ enddo; enddo
+ do J=js-2,Jeq+1 ; do I=is-2,Ieq+1
+ beta_q(I,J) = sqrt( (0.25*(G%dF_dx(i,j)+G%dF_dx(i+1,j)+G%dF_dx(i,j+1)+G%dF_dx(i+1,j+1))**2) + &
+ (0.25*(G%dF_dy(i,j)+G%dF_dy(i+1,j)+G%dF_dy(i,j+1)+G%dF_dy(i+1,j+1))**2) )
+ enddo ; enddo
+
+ do J=js-2,Jeq+1 ; do i=is-1,Ieq+1
+ vort_xy_dx(i,J) = vort_xy_dx(i,J) + 0.5 * ( G%dF_dx(i,j) + G%dF_dx(i,j+1))
+ enddo ; enddo
+ do j=js-1,Jeq+1 ; do I=is-2,Ieq+1
+ vort_xy_dy(I,j) = vort_xy_dy(I,j) + 0.5 * ( G%dF_dy(i,j) + G%dF_dy(i+1,j))
+ enddo ; enddo
+ endif ! CS%use_beta_in_Leith
+
+ if (CS%use_QG_Leith_visc) then
+
+ do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2
+ grad_vort_mag_h_2d(i,j) = SQRT((0.5*(vort_xy_dx(i,J) + vort_xy_dx(i,J-1)))**2 + (0.5*(vort_xy_dy(I,j) + &
+ vort_xy_dy(I-1,j)))**2 )
+ enddo; enddo
+ do J=js-2,Jeq+1 ; do I=is-2,Ieq+1
+ grad_vort_mag_q_2d(I,J) = SQRT((0.5*(vort_xy_dx(i,J) + vort_xy_dx(i+1,J)))**2 + (0.5*(vort_xy_dy(I,j) + &
+ vort_xy_dy(I,j+1)))**2 )
+ enddo ; enddo
+
+ call calc_QG_Leith_viscosity(VarMix, G, GV, h, k, div_xx_dx, div_xx_dy, &
+ vort_xy_dx, vort_xy_dy)
+
+ endif
+
+ do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2
+ grad_vort_mag_h(i,j) = SQRT((0.5*(vort_xy_dx(i,J) + vort_xy_dx(i,J-1)))**2 + (0.5*(vort_xy_dy(I,j) + &
+ vort_xy_dy(I-1,j)))**2 )
+ enddo; enddo
+ do J=js-2,Jeq+1 ; do I=is-2,Ieq+1
+ grad_vort_mag_q(I,J) = SQRT((0.5*(vort_xy_dx(i,J) + vort_xy_dx(i+1,J)))**2 + (0.5*(vort_xy_dy(I,j) + &
+ vort_xy_dy(I,j+1)))**2 )
+ enddo ; enddo
+
+ endif ! CS%Leith_Kh
+
meke_res_fn = 1.
- do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+ do J=Jsq,Jeq+1 ; do i=Isq,Ieq+1
if ((CS%Smagorinsky_Kh) .or. (CS%Smagorinsky_Ah)) then
Shear_mag = sqrt(sh_xx(i,j)*sh_xx(i,j) + &
0.25*((sh_xy(I-1,J-1)*sh_xy(I-1,J-1) + sh_xy(I,J)*sh_xy(I,J)) + &
(sh_xy(I-1,J)*sh_xy(I-1,J) + sh_xy(I,J-1)*sh_xy(I,J-1))))
endif
if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then
- Vort_mag = sqrt( &
- 0.5*((vort_xy_dx(i,J-1)*vort_xy_dx(i,J-1) + vort_xy_dx(i,J)*vort_xy_dx(i,J)) + &
- (vort_xy_dy(I-1,j)*vort_xy_dy(I-1,j) + vort_xy_dy(I,j)*vort_xy_dy(I,j))) + &
- 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)) + &
- (div_xx_dy(i,J)*div_xx_dy(i,J) + div_xx_dy(i,J-1)*div_xx_dy(i,J-1))))
+ if (CS%use_QG_Leith_visc) then
+ vert_vort_mag = MIN(grad_vort_mag_h(i,j) + grad_div_mag_h(i,j),3*grad_vort_mag_h_2d(i,j))
+ else
+ vert_vort_mag = grad_vort_mag_h(i,j) + grad_div_mag_h(i,j)
+ endif
endif
if (CS%better_bound_Ah .or. CS%better_bound_Kh) then
hrat_min = min(1.0, min(h_u(I,j), h_u(I-1,j), h_v(i,J), h_v(i,J-1)) / &
@@ -587,8 +827,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
! Determine the Laplacian viscosity at h points, using the
! largest value from several parameterizations.
Kh = CS%Kh_bg_xx(i,j) ! Static (pre-computed) background viscosity
- if (CS%Smagorinsky_Kh) Kh = max( Kh, CS%LAPLAC_CONST_xx(i,j) * Shear_mag )
- if (CS%Leith_Kh) Kh = max( Kh, CS%LAPLAC3_CONST_xx(i,j) * Vort_mag )
+ if (CS%Smagorinsky_Kh) Kh = max( Kh, CS%Laplac2_const_xx(i,j) * Shear_mag )
+ if (CS%Leith_Kh) Kh = max( Kh, CS%Laplac3_const_xx(i,j) * vert_vort_mag*inv_PI3)
! All viscosity contributions above are subject to resolution scaling
if (rescale_Kh) Kh = VarMix%Res_fn_h(i,j) * Kh
if (CS%res_scale_MEKE) meke_res_fn = VarMix%Res_fn_h(i,j)
@@ -610,7 +850,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
endif
endif
- if (CS%id_Kh_h>0) Kh_h(i,j,k) = Kh
+ if ((CS%id_Kh_h>0) .or. find_FrictWork) Kh_h(i,j,k) = Kh
+ if (CS%id_div_xx_h>0) div_xx_h(i,j,k) = div_xx(i,j)
str_xx(i,j) = -Kh * sh_xx(i,j)
else ! not Laplacian
@@ -631,14 +872,13 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
if ((CS%Smagorinsky_Ah) .or. (CS%Leith_Ah)) then
if (CS%Smagorinsky_Ah) then
if (CS%bound_Coriolis) then
- AhSm = Shear_mag * (CS%BIHARM_CONST_xx(i,j) + &
- CS%Biharm_Const2_xx(i,j)*Shear_mag)
+ AhSm = Shear_mag * (CS%Biharm_const_xx(i,j) + &
+ CS%Biharm_const2_xx(i,j)*Shear_mag)
else
- AhSm = CS%BIHARM_CONST_xx(i,j) * Shear_mag
+ AhSm = CS%Biharm_const_xx(i,j) * Shear_mag
endif
endif
- if (CS%Leith_Ah) &
- AhLth = Vort_mag * (CS%BIHARM5_CONST_xx(i,j))
+ if (CS%Leith_Ah) AhLth = CS%biharm5_const_xx(i,j) * vert_vort_mag * inv_PI5
Ah = MAX(MAX(CS%Ah_bg_xx(i,j), AhSm),AhLth)
if (CS%bound_Ah .and. .not.CS%better_bound_Ah) &
Ah = MIN(Ah, CS%Ah_Max_xx(i,j))
@@ -646,11 +886,13 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
Ah = CS%Ah_bg_xx(i,j)
endif ! Smagorinsky_Ah or Leith_Ah
+ if (use_MEKE_Au) Ah = Ah + MEKE%Au(i,j) ! *Add* the MEKE contribution
+
if (CS%better_bound_Ah) then
Ah = MIN(Ah, visc_bound_rem*hrat_min*CS%Ah_Max_xx(i,j))
endif
- if (CS%id_Ah_h>0) Ah_h(i,j,k) = Ah
+ if ((CS%id_Ah_h>0) .or. find_FrictWork) Ah_h(i,j,k) = Ah
str_xx(i,j) = str_xx(i,j) + Ah * &
(CS%DY_dxT(i,j)*(G%IdyCu(I,j)*u0(I,j) - G%IdyCu(I-1,j)*u0(I-1,j)) - &
@@ -664,7 +906,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
endif ! biharmonic
- str_xx(i,j) = str_xx(i,j) * (h(i,j,k) * CS%reduction_xx(i,j))
enddo ; enddo
if (CS%biharmonic) then
@@ -706,23 +947,21 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
0.25*((sh_xx(i,j)*sh_xx(i,j) + sh_xx(i+1,j+1)*sh_xx(i+1,j+1)) + &
(sh_xx(i,j+1)*sh_xx(i,j+1) + sh_xx(i+1,j)*sh_xx(i+1,j))))
endif
- if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) &
- Vort_mag = sqrt( &
- 0.5*((vort_xy_dx(i,J)*vort_xy_dx(i,J) + vort_xy_dx(i+1,J)*vort_xy_dx(i+1,J)) + &
- (vort_xy_dy(I,j)*vort_xy_dy(I,j) + vort_xy_dy(I,j+1)*vort_xy_dy(I,j+1))) + &
- 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)) + &
- (div_xx_dy(i,J)*div_xx_dy(i,J) + div_xx_dy(i+1,J)*div_xx_dy(i+1,J))))
-
+ if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then
+ if (CS%use_QG_Leith_visc) then
+ vert_vort_mag = MIN(grad_vort_mag_q(I,J) + grad_div_mag_q(I,J), 3*grad_vort_mag_q_2d(I,J))
+ else
+ vert_vort_mag = grad_vort_mag_q(I,J) + grad_div_mag_q(I,J)
+ endif
+ endif
h2uq = 4.0 * h_u(I,j) * h_u(I,j+1)
h2vq = 4.0 * h_v(i,J) * h_v(i+1,J)
- !hq = 2.0 * h2uq * h2vq / (h_neglect3 + (h2uq + h2vq) * &
- ! ((h(i,j,k) + h(i+1,j+1,k)) + (h(i,j+1,k) + h(i+1,j,k))))
- hq = 2.0 * h2uq * h2vq / (h_neglect3 + (h2uq + h2vq) * &
- ((h_u(I,j) + h_u(I,j+1)) + (h_v(i,J) + h_v(i+1,J))))
+ hq(I,J) = 2.0 * h2uq * h2vq / (h_neglect3 + (h2uq + h2vq) * &
+ ((h_u(I,j) + h_u(I,j+1)) + (h_v(i,J) + h_v(i+1,J))))
if (CS%better_bound_Ah .or. CS%better_bound_Kh) then
hrat_min = min(1.0, min(h_u(I,j), h_u(I,j+1), h_v(i,J), h_v(i+1,J)) / &
- (hq + h_neglect) )
+ (hq(I,J) + h_neglect) )
visc_bound_rem = 1.0
endif
@@ -736,12 +975,12 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
if ((G%mask2dCu(I,j) + G%mask2dCu(I,j+1)) * &
(G%mask2dCv(i,J) + G%mask2dCv(i+1,J)) == 0.0) then
! Only one of hu and hv is nonzero, so just add them.
- hq = hu + hv
+ hq(I,J) = hu + hv
hrat_min = 1.0
else
! Both hu and hv are nonzero, so take the harmonic mean.
- hq = 2.0 * (hu * hv) / ((hu + hv) + h_neglect)
- hrat_min = min(1.0, min(hu, hv) / (hq + h_neglect) )
+ hq(I,J) = 2.0 * (hu * hv) / ((hu + hv) + h_neglect)
+ hrat_min = min(1.0, min(hu, hv) / (hq(I,J) + h_neglect) )
endif
endif
endif
@@ -750,8 +989,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
! Determine the Laplacian viscosity at q points, using the
! largest value from several parameterizations.
Kh = CS%Kh_bg_xy(i,j) ! Static (pre-computed) background viscosity
- if (CS%Smagorinsky_Kh) Kh = max( Kh, CS%LAPLAC_CONST_xy(I,J) * Shear_mag )
- if (CS%Leith_Kh) Kh = max( Kh, CS%LAPLAC3_CONST_xy(I,J) * Vort_mag)
+ if (CS%Smagorinsky_Kh) Kh = max( Kh, CS%Laplac2_const_xy(I,J) * Shear_mag )
+ if (CS%Leith_Kh) Kh = max( Kh, CS%Laplac3_const_xy(I,J) * vert_vort_mag*inv_PI3)
! All viscosity contributions above are subject to resolution scaling
if (rescale_Kh) Kh = VarMix%Res_fn_q(i,j) * Kh
if (CS%res_scale_MEKE) meke_res_fn = VarMix%Res_fn_q(i,j)
@@ -778,6 +1017,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
endif
if (CS%id_Kh_q>0) Kh_q(I,J,k) = Kh
+ if (CS%id_vort_xy_q>0) vort_xy_q(I,J,k) = vort_xy(I,J)
str_xy(I,J) = -Kh * sh_xy(I,J)
else ! not Laplacian
@@ -798,20 +1038,25 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
if (CS%Smagorinsky_Ah .or. CS%Leith_Ah) then
if (CS%Smagorinsky_Ah) then
if (CS%bound_Coriolis) then
- AhSm = Shear_mag * (CS%BIHARM_CONST_xy(I,J) + &
- CS%Biharm_Const2_xy(I,J)*Shear_mag)
+ AhSm = Shear_mag * (CS%Biharm_const_xy(I,J) + &
+ CS%Biharm_const2_xy(I,J)*Shear_mag)
else
- AhSm = CS%BIHARM_CONST_xy(I,J) * Shear_mag
+ AhSm = CS%Biharm_const_xy(I,J) * Shear_mag
endif
endif
- if (CS%Leith_Ah) &
- AhLth = Vort_mag * (CS%BIHARM5_CONST_xy(I,J))
+ if (CS%Leith_Ah) AhLth = CS%Biharm5_const_xy(I,J) * vert_vort_mag * inv_PI5
Ah = MAX(MAX(CS%Ah_bg_xy(I,J), AhSm),AhLth)
if (CS%bound_Ah .and. .not.CS%better_bound_Ah) &
Ah = MIN(Ah, CS%Ah_Max_xy(I,J))
else
Ah = CS%Ah_bg_xy(I,J)
endif ! Smagorinsky_Ah or Leith_Ah
+
+ if (use_MEKE_Au) then ! *Add* the MEKE contribution
+ Ah = Ah + 0.25*( (MEKE%Au(I,J)+MEKE%Au(I+1,J+1)) &
+ +(MEKE%Au(I+1,J)+MEKE%Au(I,J+1)) )
+ endif
+
if (CS%better_bound_Ah) then
Ah = MIN(Ah, visc_bound_rem*hrat_min*CS%Ah_Max_xy(I,J))
endif
@@ -822,16 +1067,155 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
! Keep a copy of the biharmonic contribution for backscatter parameterization
bhstr_xy(I,J) = Ah * ( dvdx(I,J) + dudy(I,J) ) * &
- (hq * G%mask2dBu(I,J) * CS%reduction_xy(I,J))
+ (hq(I,J) * G%mask2dBu(I,J) * CS%reduction_xy(I,J))
endif ! biharmonic
- if (CS%no_slip) then
- str_xy(I,J) = str_xy(I,J) * (hq * CS%reduction_xy(I,J))
+ enddo ; enddo
+
+
+ if (find_FrictWork) then
+ if (CS%biharmonic) call pass_vector(u0, v0, G%Domain)
+ call pass_var(dudx, G%Domain, complete=.true.)
+ call pass_var(dvdy, G%Domain, complete=.true.)
+ call pass_var(dvdx, G%Domain, position=CORNER, complete=.true.)
+ call pass_var(dudy, G%Domain, position=CORNER, complete=.true.)
+
+ if (CS%Laplacian) then
+ do j=js,je ; do i=is,ie
+ grad_vel_mag_h(i,j) = boundary_mask(i,j) * (dudx(i,j)**2 + dvdy(i,j)**2 + &
+ (0.25*(dvdx(I,J)+dvdx(I-1,J)+dvdx(I,J-1)+dvdx(I-1,J-1)) )**2 + &
+ (0.25*(dudy(I,J)+dudy(I-1,J)+dudy(I,J-1)+dudy(I-1,J-1)) )**2)
+ enddo ; enddo
else
- str_xy(I,J) = str_xy(I,J) * (hq * G%mask2dBu(I,J) * CS%reduction_xy(I,J))
+ do j=js,je ; do i=is,ie
+ grad_vel_mag_h(i,j) = 0.0
+ enddo ; enddo
endif
- enddo ; enddo
+
+ if (CS%biharmonic) then
+ do j=js,je ; do i=is,ie
+ grad_d2vel_mag_h(i,j) = boundary_mask(i,j) * ((0.5*(u0(I,j) + u0(I-1,j)))**2 + &
+ (0.5*(v0(i,J) + v0(i,J-1)))**2)
+ enddo ; enddo
+ else
+ do j=js,je ; do i=is,ie
+ grad_d2vel_mag_h(i,j) = 0.0
+ enddo ; enddo
+ endif
+
+ do j=js,je ; do i=is,ie
+ ! Diagnose -Kh * |del u|^2 - Ah * |del^2 u|^2
+ diss_rate(i,j,k) = -Kh_h(i,j,k) * grad_vel_mag_h(i,j) - &
+ Ah_h(i,j,k) * grad_d2vel_mag_h(i,j)
+
+ if (associated(MEKE)) then ; if (associated(MEKE%mom_src)) then
+ ! This is the maximum possible amount of energy that can be converted
+ ! per unit time, according to theory (multiplied by h)
+ max_diss_rate(i,j,k) = 2.0*MEKE%MEKE(i,j) * sqrt(grad_vel_mag_h(i,j))
+ FrictWork_diss(i,j,k) = diss_rate(i,j,k) * h(i,j,k) * GV%H_to_kg_m2
+ FrictWorkMax(i,j,k) = -max_diss_rate(i,j,k) * h(i,j,k) * GV%H_to_kg_m2
+
+ ! Determine how much work GME needs to do to reach the "target" ratio between
+ ! the amount of work actually done and the maximum allowed by theory. Note that
+ ! we need to add the FrictWork done by the dissipation operators, since this work
+ ! is done only for numerical stability and is therefore spurious
+ if (CS%use_GME) then
+ target_diss_rate_GME(i,j,k) = FWfrac * max_diss_rate(i,j,k) - diss_rate(i,j,k)
+ endif
+
+ else
+ FrictWork_diss(i,j,k) = diss_rate(i,j,k) * h(i,j,k) * GV%H_to_kg_m2
+ endif ; endif
+
+ enddo ; enddo
+ endif
+
+
+ if (CS%use_GME) then
+
+ if (.not. (associated(MEKE))) call MOM_error(FATAL, &
+ "MEKE must be enabled for GME to be used.")
+
+ do J=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+
+ if ((max_diss_rate(i,j,k) > 0) .and. (grad_vel_mag_bt_h(i,j)>0) ) then
+ GME_coeff = (MIN(G%bathyT(i,j)/H0,1.0)**2) * FWfrac*max_diss_rate(i,j,k) / grad_vel_mag_bt_h(i,j)
+! GME_coeff = (MIN(G%bathyT(i,j)/H0,1.0)**2) * FWfrac*target_diss_rate_GME(i,j,k) / grad_vel_mag_bt_h(i,j)
+ else
+ GME_coeff = 0.0
+ endif
+
+ ! apply mask
+ GME_coeff = GME_coeff * boundary_mask(i,j)
+
+ GME_coeff = MIN(GME_coeff, GME_coeff_limiter)
+
+ if ((CS%id_GME_coeff_h>0) .or. find_FrictWork) GME_coeff_h(i,j,k) = GME_coeff
+
+ str_xx_GME(i,j) = GME_coeff * sh_xx_bt(i,j)
+
+ enddo ; enddo
+
+
+ do J=js-1,Jeq ; do I=is-1,Ieq
+
+ if ((max_diss_rate(i,j,k) > 0) .and. (grad_vel_mag_bt_q(i,j)>0) ) then
+ GME_coeff = (MIN(G%bathyT(i,j)/H0,1.0)**2) * FWfrac*max_diss_rate(i,j,k) / grad_vel_mag_bt_q(I,J)
+! GME_coeff = (MIN(G%bathyT(i,j)/H0,1.0)**2) * FWfrac*target_diss_rate_GME(i,j,k) / grad_vel_mag_bt_q(I,J)
+ else
+ GME_coeff = 0.0
+ endif
+
+ ! apply mask
+ GME_coeff = GME_coeff * boundary_mask(i,j)
+
+ GME_coeff = MIN(GME_coeff, GME_coeff_limiter)
+
+ if (CS%id_GME_coeff_q>0) GME_coeff_q(I,J,k) = GME_coeff
+ str_xy_GME(I,J) = GME_coeff * sh_xy_bt(I,J)
+
+ enddo ; enddo
+
+ ! applying GME diagonal term
+ call smooth_GME(CS,G,GME_flux_h=str_xx_GME)
+ call smooth_GME(CS,G,GME_flux_q=str_xy_GME)
+
+ do J=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+ str_xx(i,j) = (str_xx(i,j) + str_xx_GME(i,j)) * (h(i,j,k) * CS%reduction_xx(i,j))
+ enddo ; enddo
+
+ do J=js-1,Jeq ; do I=is-1,Ieq
+ ! GME is applied below
+ if (CS%no_slip) then
+ str_xy(I,J) = (str_xy(I,J) + str_xy_GME(I,J)) * (hq(I,J) * CS%reduction_xy(I,J))
+ else
+ str_xy(I,J) = (str_xy(I,J) + str_xy_GME(I,J)) * (hq(I,J) * G%mask2dBu(I,J) * CS%reduction_xy(I,J))
+ endif
+ enddo ; enddo
+
+ if (associated(MEKE%GME_snk)) then
+ do j=js,je ; do i=is,ie
+ FrictWork_GME(i,j,k) = GME_coeff_h(i,j,k) * h(i,j,k) * GV%H_to_kg_m2 * grad_vel_mag_bt_h(i,j)
+ enddo ; enddo
+ endif
+
+ else ! use_GME
+ do J=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+ str_xx(i,j) = str_xx(i,j) * (h(i,j,k) * CS%reduction_xx(i,j))
+ enddo ; enddo
+
+ do J=js-1,Jeq ; do I=is-1,Ieq
+ if (CS%no_slip) then
+ str_xy(I,J) = str_xy(I,J) * (hq(I,J) * CS%reduction_xy(I,J))
+ else
+ str_xy(I,J) = str_xy(I,J) * (hq(I,J) * G%mask2dBu(I,J) * CS%reduction_xy(I,J))
+ endif
+ enddo ; enddo
+
+ endif ! use_GME
+
+
! Evaluate 1/h x.Div(h Grad u) or the biharmonic equivalent.
do j=js,je ; do I=Isq,Ieq
@@ -877,7 +1261,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
endif
if (find_FrictWork) then ; do j=js,je ; do i=is,ie
- ! Diagnose str_xx*d_x u + str_yy*d_y v + str_xy*(d_y u + d_x v)
+ ! Diagnose str_xx*d_x u - str_yy*d_y v + str_xy*(d_y u + d_x v)
+ ! This is the old formulation that includes energy diffusion
FrictWork(i,j,k) = GV%H_to_kg_m2 * ( &
(str_xx(i,j)*(u(I,j,k)-u(I-1,j,k))*G%IdxT(i,j) &
-str_xx(i,j)*(v(i,J,k)-v(i,J-1,k))*G%IdyT(i,j)) &
@@ -902,6 +1287,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
if (k==1) then
do j=js,je ; do i=is,ie
MEKE%mom_src(i,j) = 0.
+ MEKE%GME_snk(i,j) = 0.
enddo ; enddo
endif
if (MEKE%backscatter_Ro_c /= 0.) then
@@ -936,10 +1322,20 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
enddo ; enddo
else
do j=js,je ; do i=is,ie
- MEKE%mom_src(i,j) = MEKE%mom_src(i,j) + FrictWork(i,j,k)
+ ! MEKE%mom_src now is sign definite because it only uses the dissipation
+ MEKE%mom_src(i,j) = MEKE%mom_src(i,j) + MAX(FrictWork_diss(i,j,k), FrictWorkMax(i,j,k))
enddo ; enddo
+ endif ! MEKE%backscatter
+
+ if (CS%use_GME .and. associated(MEKE)) then
+ if (associated(MEKE%GME_snk)) then
+ do j=js,je ; do i=is,ie
+ MEKE%GME_snk(i,j) = MEKE%GME_snk(i,j) + FrictWork_GME(i,j,k)
+ enddo ; enddo
+ endif
endif
- endif ; endif
+
+ endif ; endif ! find_FrictWork and associated(mom_src)
enddo ! end of k loop
@@ -947,10 +1343,17 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
if (CS%id_diffu>0) call post_data(CS%id_diffu, diffu, CS%diag)
if (CS%id_diffv>0) call post_data(CS%id_diffv, diffv, CS%diag)
if (CS%id_FrictWork>0) call post_data(CS%id_FrictWork, FrictWork, CS%diag)
+ if (CS%id_FrictWorkMax>0) call post_data(CS%id_FrictWorkMax, FrictWorkMax, CS%diag)
+ if (CS%id_FrictWork_diss>0) call post_data(CS%id_FrictWork_diss, FrictWork_diss, CS%diag)
+ if (CS%id_FrictWork_GME>0) call post_data(CS%id_FrictWork_GME, FrictWork_GME, CS%diag)
if (CS%id_Ah_h>0) call post_data(CS%id_Ah_h, Ah_h, CS%diag)
+ if (CS%id_div_xx_h>0) call post_data(CS%id_div_xx_h, div_xx_h, CS%diag)
+ if (CS%id_vort_xy_q>0) call post_data(CS%id_vort_xy_q, vort_xy_q, CS%diag)
if (CS%id_Ah_q>0) call post_data(CS%id_Ah_q, Ah_q, CS%diag)
if (CS%id_Kh_h>0) call post_data(CS%id_Kh_h, Kh_h, CS%diag)
if (CS%id_Kh_q>0) call post_data(CS%id_Kh_q, Kh_q, CS%diag)
+ if (CS%id_GME_coeff_h > 0) call post_data(CS%id_GME_coeff_h, GME_coeff_h, CS%diag)
+ if (CS%id_GME_coeff_q > 0) call post_data(CS%id_GME_coeff_q, GME_coeff_q, CS%diag)
if (CS%id_FrictWorkIntz > 0) then
do j=js,je
@@ -1013,6 +1416,8 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS)
logical :: get_all ! If true, read and log all parameters, regardless of
! whether they are used, to enable spell-checking of
! valid parameters.
+ logical :: split ! If true, use the split time stepping scheme.
+ ! If false and USE_GME = True, issue a FATAL error.
character(len=64) :: inputdir, filename
real :: deg2rad ! Converts degrees to radians
real :: slat_fn ! sin(lat)**Kh_pwr_of_sine
@@ -1047,6 +1452,7 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS)
! cases where the corresponding parameters are not read.
CS%bound_Kh = .false. ; CS%better_bound_Kh = .false. ; CS%Smagorinsky_Kh = .false. ; CS%Leith_Kh = .false.
CS%bound_Ah = .false. ; CS%better_bound_Ah = .false. ; CS%Smagorinsky_Ah = .false. ; CS%Leith_Ah = .false.
+ CS%use_QG_Leith_visc = .false.
CS%bound_Coriolis = .false.
CS%Modified_Leith = .false.
CS%anisotropic = .false.
@@ -1105,12 +1511,27 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS)
"If true, the viscosity contribution from MEKE is scaled by "//&
"the resolution function.", default=.false.)
- if (CS%Leith_Kh .or. get_all) &
+ if (CS%Leith_Kh .or. get_all) then
call get_param(param_file, mdl, "LEITH_LAP_CONST", Leith_Lap_const, &
"The nondimensional Laplacian Leith constant, "//&
- "often ??", units="nondim", default=0.0, &
+ "often set to 1.0", units="nondim", default=0.0, &
fail_if_missing = CS%Leith_Kh)
-
+ call get_param(param_file, mdl, "USE_QG_LEITH_VISC", CS%use_QG_Leith_visc, &
+ "If true, use QG Leith nonlinear eddy viscosity.", &
+ default=.false.)
+ if (CS%use_QG_Leith_visc .and. .not. CS%Leith_Kh) call MOM_error(FATAL, &
+ "MOM_lateral_mixing_coeffs.F90, VarMix_init:"//&
+ "LEITH_KH must be True when USE_QG_LEITH_VISC=True.")
+ endif
+ if (CS%Leith_Kh .or. CS%Leith_Ah .or. get_all) then
+ call get_param(param_file, mdl, "USE_BETA_IN_LEITH", CS%use_beta_in_Leith, &
+ "If true, include the beta term in the Leith nonlinear eddy viscosity.", &
+ default=CS%Leith_Kh)
+ call get_param(param_file, mdl, "MODIFIED_LEITH", CS%modified_Leith, &
+ "If true, add a term to Leith viscosity which is \n"//&
+ "proportional to the gradient of divergence.", &
+ default=.false.)
+ endif
call get_param(param_file, mdl, "BOUND_KH", CS%bound_Kh, &
"If true, the Laplacian coefficient is locally limited "//&
"to be stable.", default=.true.)
@@ -1206,12 +1627,11 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS)
endif
endif
- if (CS%Leith_Ah .or. get_all) then
- call get_param(param_file, mdl, "LEITH_BI_CONST",Leith_bi_const, &
+ if (CS%Leith_Ah .or. get_all) &
+ call get_param(param_file, mdl, "LEITH_BI_CONST", Leith_bi_const, &
"The nondimensional biharmonic Leith constant, "//&
- "typical values are thus far undetermined", units="nondim", default=0.0, &
+ "typical values are thus far undetermined.", units="nondim", default=0.0, &
fail_if_missing = CS%Leith_Ah)
- endif
endif
@@ -1242,6 +1662,17 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS)
"viscosities. The final viscosity is the maximum of the other "//&
"terms and this background value.", default=.false.)
+ call get_param(param_file, mdl, "USE_GME", CS%use_GME, &
+ "If true, use the GM+E backscatter scheme in association \n"//&
+ "with the Gent and McWilliams parameterization.", default=.false.)
+
+ if (CS%use_GME) then
+ call get_param(param_file, mdl, "SPLIT", split, &
+ "Use the split time stepping if true.", default=.true., &
+ do_not_log=.true.)
+ if (.not. split) call MOM_error(FATAL,"ERROR: Currently, USE_GME = True "// &
+ "cannot be used with SPLIT=False.")
+ endif
if (CS%bound_Kh .or. CS%bound_Ah .or. CS%better_bound_Kh .or. CS%better_bound_Ah) &
call get_param(param_file, mdl, "DT", dt, &
@@ -1279,14 +1710,13 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS)
ALLOC_(CS%Kh_Max_xy(IsdB:IedB,JsdB:JedB)) ; CS%Kh_Max_xy(:,:) = 0.0
endif
if (CS%Smagorinsky_Kh) then
- ALLOC_(CS%Laplac_Const_xx(isd:ied,jsd:jed)) ; CS%Laplac_Const_xx(:,:) = 0.0
- ALLOC_(CS%Laplac_Const_xy(IsdB:IedB,JsdB:JedB)) ; CS%Laplac_Const_xy(:,:) = 0.0
+ ALLOC_(CS%Laplac2_const_xx(isd:ied,jsd:jed)) ; CS%Laplac2_const_xx(:,:) = 0.0
+ ALLOC_(CS%Laplac2_const_xy(IsdB:IedB,JsdB:JedB)) ; CS%Laplac2_const_xy(:,:) = 0.0
endif
if (CS%Leith_Kh) then
- ALLOC_(CS%Laplac3_Const_xx(isd:ied,jsd:jed)) ; CS%Laplac3_Const_xx(:,:) = 0.0
- ALLOC_(CS%Laplac3_Const_xy(IsdB:IedB,JsdB:JedB)) ; CS%Laplac3_Const_xy(:,:) = 0.0
+ ALLOC_(CS%Laplac3_const_xx(isd:ied,jsd:jed)) ; CS%Laplac3_const_xx(:,:) = 0.0
+ ALLOC_(CS%Laplac3_const_xy(IsdB:IedB,JsdB:JedB)) ; CS%Laplac3_const_xy(:,:) = 0.0
endif
-
endif
ALLOC_(CS%reduction_xx(isd:ied,jsd:jed)) ; CS%reduction_xx(:,:) = 0.0
ALLOC_(CS%reduction_xy(IsdB:IedB,JsdB:JedB)) ; CS%reduction_xy(:,:) = 0.0
@@ -1334,16 +1764,16 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS)
ALLOC_(CS%Ah_Max_xy(IsdB:IedB,JsdB:JedB)) ; CS%Ah_Max_xy(:,:) = 0.0
endif
if (CS%Smagorinsky_Ah) then
- ALLOC_(CS%Biharm_Const_xx(isd:ied,jsd:jed)) ; CS%Biharm_Const_xx(:,:) = 0.0
- ALLOC_(CS%Biharm_Const_xy(IsdB:IedB,JsdB:JedB)) ; CS%Biharm_Const_xy(:,:) = 0.0
+ ALLOC_(CS%Biharm_const_xx(isd:ied,jsd:jed)) ; CS%Biharm_const_xx(:,:) = 0.0
+ ALLOC_(CS%Biharm_const_xy(IsdB:IedB,JsdB:JedB)) ; CS%Biharm_const_xy(:,:) = 0.0
if (CS%bound_Coriolis) then
- ALLOC_(CS%Biharm_Const2_xx(isd:ied,jsd:jed)) ; CS%Biharm_Const2_xx(:,:) = 0.0
- ALLOC_(CS%Biharm_Const2_xy(IsdB:IedB,JsdB:JedB)) ; CS%Biharm_Const2_xy(:,:) = 0.0
+ ALLOC_(CS%Biharm_const2_xx(isd:ied,jsd:jed)) ; CS%Biharm_const2_xx(:,:) = 0.0
+ ALLOC_(CS%Biharm_const2_xy(IsdB:IedB,JsdB:JedB)) ; CS%Biharm_const2_xy(:,:) = 0.0
endif
endif
if (CS%Leith_Ah) then
- ALLOC_(CS%Biharm5_Const_xx(isd:ied,jsd:jed)) ; CS%Biharm5_Const_xx(:,:) = 0.0
- ALLOC_(CS%Biharm5_Const_xy(IsdB:IedB,JsdB:JedB)) ; CS%Biharm5_Const_xy(:,:) = 0.0
+ ALLOC_(CS%biharm5_const_xx(isd:ied,jsd:jed)) ; CS%biharm5_const_xx(:,:) = 0.0
+ ALLOC_(CS%biharm5_const_xy(IsdB:IedB,JsdB:JedB)) ; CS%biharm5_const_xy(:,:) = 0.0
endif
endif
@@ -1398,9 +1828,8 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS)
! Static factors in the Smagorinsky and Leith schemes
grid_sp_h2 = (2.0*CS%DX2h(i,j)*CS%DY2h(i,j)) / (CS%DX2h(i,j) + CS%DY2h(i,j))
grid_sp_h3 = grid_sp_h2*sqrt(grid_sp_h2)
- if (CS%Smagorinsky_Kh) CS%LAPLAC_CONST_xx(i,j) = Smag_Lap_const * grid_sp_h2
- if (CS%Leith_Kh) CS%LAPLAC3_CONST_xx(i,j) = Leith_Lap_const * grid_sp_h3
-
+ if (CS%Smagorinsky_Kh) CS%Laplac2_const_xx(i,j) = Smag_Lap_const * grid_sp_h2
+ if (CS%Leith_Kh) CS%Laplac3_const_xx(i,j) = Leith_Lap_const * grid_sp_h3
! Maximum of constant background and MICOM viscosity
CS%Kh_bg_xx(i,j) = MAX(Kh, Kh_vel_scale * sqrt(grid_sp_h2))
@@ -1425,9 +1854,8 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS)
! Static factors in the Smagorinsky and Leith schemes
grid_sp_q2 = (2.0*CS%DX2q(I,J)*CS%DY2q(I,J)) / (CS%DX2q(I,J) + CS%DY2q(I,J))
grid_sp_q3 = grid_sp_q2*sqrt(grid_sp_q2)
- if (CS%Smagorinsky_Kh) CS%LAPLAC_CONST_xy(I,J) = Smag_Lap_const * grid_sp_q2
- if (CS%Leith_Kh) CS%LAPLAC3_CONST_xy(I,J) = Leith_Lap_const * grid_sp_q3
-
+ if (CS%Smagorinsky_Kh) CS%Laplac2_const_xy(I,J) = Smag_Lap_const * grid_sp_q2
+ if (CS%Leith_Kh) CS%Laplac3_const_xy(I,J) = Leith_Lap_const * grid_sp_q3
! Maximum of constant background and MICOM viscosity
CS%Kh_bg_xy(I,J) = MAX(Kh, Kh_vel_scale * sqrt(grid_sp_q2))
@@ -1462,7 +1890,7 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS)
CS%Ah_bg_xy(:,:) = 0.0
! The 0.3 below was 0.4 in MOM1.10. The change in hq requires
! this to be less than 1/3, rather than 1/2 as before.
- Ah_Limit = 0.3 / (dt*64.0)
+ if (CS%better_bound_Ah .or. CS%bound_Ah) Ah_Limit = 0.3 / (dt*64.0)
if (CS%Smagorinsky_Ah .and. CS%bound_Coriolis) &
BoundCorConst = 1.0 / (5.0*(bound_Cor_vel*bound_Cor_vel))
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
@@ -1470,18 +1898,17 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS)
grid_sp_h3 = grid_sp_h2*sqrt(grid_sp_h2)
if (CS%Smagorinsky_Ah) then
- CS%BIHARM_CONST_xx(i,j) = Smag_bi_const * (grid_sp_h2 * grid_sp_h2)
+ CS%Biharm_const_xx(i,j) = Smag_bi_const * (grid_sp_h2 * grid_sp_h2)
if (CS%bound_Coriolis) then
- fmax = US%s_to_T*MAX(abs(G%CoriolisBu(I-1,J-1)), abs(G%CoriolisBu(I,J-1)), &
- abs(G%CoriolisBu(I-1,J)), abs(G%CoriolisBu(I,J)))
- CS%Biharm_Const2_xx(i,j) = (grid_sp_h2 * grid_sp_h2 * grid_sp_h2) * &
- (fmax * BoundCorConst)
+ fmax = MAX(abs(G%CoriolisBu(I-1,J-1)), abs(G%CoriolisBu(I,J-1)), &
+ abs(G%CoriolisBu(I-1,J)), abs(G%CoriolisBu(I,J)))
+ CS%Biharm_const2_xx(i,j) = (grid_sp_h2 * grid_sp_h2 * grid_sp_h2) * &
+ (fmax * BoundCorConst)
endif
endif
if (CS%Leith_Ah) then
- CS%BIHARM5_CONST_xx(i,j) = Leith_bi_const * (grid_sp_h2 * grid_sp_h3)
+ CS%biharm5_const_xx(i,j) = Leith_bi_const * (grid_sp_h3 * grid_sp_h2)
endif
-
CS%Ah_bg_xx(i,j) = MAX(Ah, Ah_vel_scale * grid_sp_h2 * sqrt(grid_sp_h2))
if (Ah_time_scale>0.) CS%Ah_bg_xx(i,j) = &
MAX(CS%Ah_bg_xx(i,j), (grid_sp_h2 * grid_sp_h2) / Ah_time_scale)
@@ -1495,15 +1922,14 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS)
grid_sp_q3 = grid_sp_q2*sqrt(grid_sp_q2)
if (CS%Smagorinsky_Ah) then
- CS%BIHARM_CONST_xy(I,J) = Smag_bi_const * (grid_sp_q2 * grid_sp_q2)
+ CS%Biharm_const_xy(I,J) = Smag_bi_const * (grid_sp_q2 * grid_sp_q2)
if (CS%bound_Coriolis) then
- CS%Biharm_Const2_xy(I,J) = (grid_sp_q2 * grid_sp_q2 * grid_sp_q2) * &
- (abs(US%s_to_T*G%CoriolisBu(I,J)) * BoundCorConst)
+ CS%Biharm_const2_xy(I,J) = (grid_sp_q2 * grid_sp_q2 * grid_sp_q2) * &
+ (abs(US%s_to_T*G%CoriolisBu(I,J)) * BoundCorConst)
endif
endif
-
if (CS%Leith_Ah) then
- CS%BIHARM5_CONST_xy(I,J) = Leith_bi_const * (grid_sp_q2 * grid_sp_q3)
+ CS%biharm5_const_xy(i,j) = Leith_bi_const * (grid_sp_q3 * grid_sp_q2)
endif
CS%Ah_bg_xy(I,J) = MAX(Ah, Ah_vel_scale * grid_sp_q2 * sqrt(grid_sp_q2))
@@ -1627,11 +2053,37 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS)
CS%id_Kh_q = register_diag_field('ocean_model', 'Khq', diag%axesBL, Time, &
'Laplacian Horizontal Viscosity at q Points', 'm2 s-1')
+
+ if (CS%Leith_Kh) then
+ CS%id_vort_xy_q = register_diag_field('ocean_model', 'vort_xy_q', diag%axesBL, Time, &
+ 'Vertical vorticity at q Points', 's-1')
+
+ CS%id_div_xx_h = register_diag_field('ocean_model', 'div_xx_h', diag%axesTL, Time, &
+ 'Horizontal divergence at h Points', 's-1')
+ endif
+
+ endif
+
+ if (CS%use_GME) then
+ CS%id_GME_coeff_h = register_diag_field('ocean_model', 'GME_coeff_h', diag%axesTL, Time, &
+ 'GME coefficient at h Points', 'm^2 s-1')
+
+ CS%id_GME_coeff_q = register_diag_field('ocean_model', 'GME_coeff_q', diag%axesBL, Time, &
+ 'GME coefficient at q Points', 'm^2 s-1')
+
+ CS%id_FrictWork_GME = register_diag_field('ocean_model','FrictWork_GME',diag%axesTL,Time,&
+ 'Integral work done by lateral friction terms in GME (excluding diffusion of energy)', 'W m-2')
endif
CS%id_FrictWork = register_diag_field('ocean_model','FrictWork',diag%axesTL,Time,&
'Integral work done by lateral friction terms', 'W m-2')
+ CS%id_FrictWork_diss = register_diag_field('ocean_model','FrictWork_diss',diag%axesTL,Time,&
+ 'Integral work done by lateral friction terms (excluding diffusion of energy)', 'W m-2')
+
+ CS%id_FrictWorkMax = register_diag_field('ocean_model','FrictWorkMax',diag%axesTL,Time,&
+ 'Maximum possible integral work done by lateral friction terms', 'W m-2')
+
CS%id_FrictWorkIntz = register_diag_field('ocean_model','FrictWorkIntz',diag%axesT1,Time, &
'Depth integrated work done by lateral friction', 'W m-2', &
cmor_field_name='dispkexyfo', &
@@ -1663,6 +2115,80 @@ subroutine align_aniso_tensor_to_grid(CS, n1, n2)
end subroutine align_aniso_tensor_to_grid
+!> Apply a 1-1-4-1-1 Laplacian filter one time on GME diffusive flux to reduce any
+!! horizontal two-grid-point noise
+subroutine smooth_GME(CS,G,GME_flux_h,GME_flux_q)
+ ! Arguments
+ type(hor_visc_CS), pointer :: CS !< Control structure
+ type(ocean_grid_type), intent(in) :: G !< Ocean grid
+ real, dimension(SZI_(G),SZJ_(G)), optional, intent(inout) :: GME_flux_h!< GME diffusive flux
+ !! at h points
+ real, dimension(SZIB_(G),SZJB_(G)), optional, intent(inout) :: GME_flux_q!< GME diffusive flux
+ !! at q points
+
+ ! local variables
+ real, dimension(SZI_(G),SZJ_(G)) :: GME_flux_h_original
+ real, dimension(SZIB_(G),SZJB_(G)) :: GME_flux_q_original
+ real :: wc, ww, we, wn, ws ! averaging weights for smoothing
+ integer :: i, j, k, s
+
+ !do s=1,CS%n_smooth
+ do s=1,1
+
+ ! Update halos
+ if (present(GME_flux_h)) then
+ call pass_var(GME_flux_h, G%Domain)
+ GME_flux_h_original = GME_flux_h
+ ! apply smoothing on GME
+ do j = G%jsc, G%jec
+ do i = G%isc, G%iec
+ ! skip land points
+ if (G%mask2dT(i,j)==0.) cycle
+
+ ! compute weights
+ ww = 0.125 * G%mask2dT(i-1,j)
+ we = 0.125 * G%mask2dT(i+1,j)
+ ws = 0.125 * G%mask2dT(i,j-1)
+ wn = 0.125 * G%mask2dT(i,j+1)
+ wc = 1.0 - (ww+we+wn+ws)
+
+ GME_flux_h(i,j) = wc * GME_flux_h_original(i,j) &
+ + ww * GME_flux_h_original(i-1,j) &
+ + we * GME_flux_h_original(i+1,j) &
+ + ws * GME_flux_h_original(i,j-1) &
+ + wn * GME_flux_h_original(i,j+1)
+ enddo; enddo
+ endif
+
+ ! Update halos
+ if (present(GME_flux_q)) then
+ call pass_var(GME_flux_q, G%Domain, position=CORNER, complete=.true.)
+ GME_flux_q_original = GME_flux_q
+ ! apply smoothing on GME
+ do J = G%JscB, G%JecB
+ do I = G%IscB, G%IecB
+ ! skip land points
+ if (G%mask2dBu(I,J)==0.) cycle
+
+ ! compute weights
+ ww = 0.125 * G%mask2dBu(I-1,J)
+ we = 0.125 * G%mask2dBu(I+1,J)
+ ws = 0.125 * G%mask2dBu(I,J-1)
+ wn = 0.125 * G%mask2dBu(I,J+1)
+ wc = 1.0 - (ww+we+wn+ws)
+
+ GME_flux_q(I,J) = wc * GME_flux_q_original(I,J) &
+ + ww * GME_flux_q_original(I-1,J) &
+ + we * GME_flux_q_original(I+1,J) &
+ + ws * GME_flux_q_original(I,J-1) &
+ + wn * GME_flux_q_original(I,J+1)
+ enddo; enddo
+ endif
+
+ enddo ! s-loop
+
+end subroutine smooth_GME
+
!> Deallocates any variables allocated in hor_visc_init.
subroutine hor_visc_end(CS)
type(hor_visc_CS), pointer :: CS !< The control structure returned by a
@@ -1680,10 +2206,10 @@ subroutine hor_visc_end(CS)
DEALLOC_(CS%Kh_Max_xx) ; DEALLOC_(CS%Kh_Max_xy)
endif
if (CS%Smagorinsky_Kh) then
- DEALLOC_(CS%Laplac_Const_xx) ; DEALLOC_(CS%Laplac_Const_xy)
+ DEALLOC_(CS%Laplac2_const_xx) ; DEALLOC_(CS%Laplac2_const_xy)
endif
if (CS%Leith_Kh) then
- DEALLOC_(CS%Laplac3_Const_xx) ; DEALLOC_(CS%Laplac3_Const_xy)
+ DEALLOC_(CS%Laplac3_const_xx) ; DEALLOC_(CS%Laplac3_const_xy)
endif
endif
@@ -1695,13 +2221,13 @@ subroutine hor_visc_end(CS)
DEALLOC_(CS%Ah_Max_xx) ; DEALLOC_(CS%Ah_Max_xy)
endif
if (CS%Smagorinsky_Ah) then
- DEALLOC_(CS%Biharm_Const_xx) ; DEALLOC_(CS%Biharm_Const_xy)
+ DEALLOC_(CS%Biharm5_const_xx) ; DEALLOC_(CS%Biharm5_const_xy)
if (CS%bound_Coriolis) then
- DEALLOC_(CS%Biharm_Const2_xx) ; DEALLOC_(CS%Biharm_Const2_xy)
+ DEALLOC_(CS%Biharm5_const2_xx) ; DEALLOC_(CS%Biharm5_const2_xy)
endif
endif
if (CS%Leith_Ah) then
- DEALLOC_(CS%Biharm5_Const_xx) ; DEALLOC_(CS%Biharm5_Const_xy)
+ DEALLOC_(CS%Biharm_const_xx) ; DEALLOC_(CS%Biharm_const_xy)
endif
endif
if (CS%anisotropic) then
diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90
index 2a855f4416..d2a3d6d1f6 100644
--- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90
+++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90
@@ -8,7 +8,7 @@ module MOM_lateral_mixing_coeffs
use MOM_diag_mediator, only : register_diag_field, safe_alloc_ptr, post_data
use MOM_diag_mediator, only : diag_ctrl, time_type, query_averaging_enabled
use MOM_domains, only : create_group_pass, do_group_pass
-use MOM_domains, only : group_pass_type, pass_var
+use MOM_domains, only : group_pass_type, pass_var, pass_vector
use MOM_file_parser, only : get_param, log_version, param_file_type
use MOM_interface_heights, only : find_eta
use MOM_isopycnal_slopes, only : calc_isoneutral_slopes
@@ -90,9 +90,23 @@ module MOM_lateral_mixing_coeffs
real, dimension(:,:,:), pointer :: &
slope_x => NULL(), & !< Zonal isopycnal slope [nondim]
slope_y => NULL(), & !< Meridional isopycnal slope [nondim]
+ N2_u => NULL(), & !< Brunt-Vaisala frequency at u-points [s-2]
+ N2_v => NULL(), & !< Brunt-Vaisala frequency at v-points [s-2]
ebt_struct => NULL() !< Vertical structure function to scale diffusivities with [nondim]
+ real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_) :: &
+ Laplac3_const_u !< Laplacian metric-dependent constants [nondim]
+
+ real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_) :: &
+ Laplac3_const_v !< Laplacian metric-dependent constants [nondim]
+
+ real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: &
+ KH_u_QG !< QG Leith GM coefficient at u-points [m2 s-1]
+
+ real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: &
+ KH_v_QG !< QG Leith GM coefficient at v-points [m2 s-1]
! Parameters
+ logical :: use_Visbeck !< Use Visbeck formulation for thickness diffusivity
integer :: VarMix_Ktop !< Top layer to start downward integrals
real :: Visbeck_L_scale !< Fixed length scale in Visbeck formula
real :: Res_coef_khth !< A non-dimensional number that determines the function
@@ -110,12 +124,16 @@ module MOM_lateral_mixing_coeffs
!! and especially 2 are coded to be more efficient.
real :: Visbeck_S_max !< Upper bound on slope used in Eady growth rate [nondim].
+ ! Leith parameters
+ logical :: use_QG_Leith_GM !< If true, uses the QG Leith viscosity as the GM coefficient
+ logical :: use_beta_in_QG_Leith !< If true, includes the beta term in the QG Leith GM coefficient
+
! Diagnostics
!>@{
!! Diagnostic identifier
integer :: id_SN_u=-1, id_SN_v=-1, id_L2u=-1, id_L2v=-1, id_Res_fn = -1
integer :: id_N2_u=-1, id_N2_v=-1, id_S2_u=-1, id_S2_v=-1
- integer :: id_Rd_dx=-1
+ integer :: id_Rd_dx=-1, id_KH_u_QG = -1, id_KH_v_QG = -1
type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the
!! timing of diagnostic output.
!>@}
@@ -126,6 +144,7 @@ module MOM_lateral_mixing_coeffs
end type VarMix_CS
public VarMix_init, calc_slope_functions, calc_resoln_function
+public calc_QG_Leith_viscosity
contains
@@ -409,12 +428,12 @@ subroutine calc_slope_functions(h, tv, dt, G, GV, US, CS)
endif
if (query_averaging_enabled(CS%diag)) then
- if (CS%id_SN_u > 0) call post_data(CS%id_SN_u, CS%SN_u, CS%diag)
- if (CS%id_SN_v > 0) call post_data(CS%id_SN_v, CS%SN_v, CS%diag)
- if (CS%id_L2u > 0) call post_data(CS%id_L2u, CS%L2u, CS%diag)
- if (CS%id_L2v > 0) call post_data(CS%id_L2v, CS%L2v, CS%diag)
- if (CS%id_N2_u > 0) call post_data(CS%id_N2_u, N2_u, CS%diag)
- if (CS%id_N2_v > 0) call post_data(CS%id_N2_v, N2_v, CS%diag)
+ if (CS%id_SN_u > 0) call post_data(CS%id_SN_u, CS%SN_u, CS%diag)
+ if (CS%id_SN_v > 0) call post_data(CS%id_SN_v, CS%SN_v, CS%diag)
+ if (CS%id_L2u > 0) call post_data(CS%id_L2u, CS%L2u, CS%diag)
+ if (CS%id_L2v > 0) call post_data(CS%id_L2v, CS%L2v, CS%diag)
+ if (CS%id_N2_u > 0) call post_data(CS%id_N2_u, CS%N2_u, CS%diag)
+ if (CS%id_N2_v > 0) call post_data(CS%id_N2_v, CS%N2_v, CS%diag)
endif
end subroutine calc_slope_functions
@@ -707,6 +726,154 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slop
end subroutine calc_slope_functions_using_just_e
+!> Calculates the Leith Laplacian and bi-harmonic viscosity coefficients
+subroutine calc_QG_Leith_viscosity(CS, G, GV, h, k, div_xx_dx, div_xx_dy, vort_xy_dx, vort_xy_dy)
+ type(VarMix_CS), pointer :: CS !< Variable mixing coefficients
+ type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
+ type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
+! real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< Zonal flow (m s-1)
+! real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< Meridional flow (m s-1)
+ real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: h !< Layer thickness (m or kg m-2)
+ integer, intent(in) :: k !< Layer for which to calculate vorticity magnitude
+ real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: div_xx_dx !< x-derivative of horizontal divergence (d/dx(du/dx + dv/dy)) (m-1 s-1)
+ real, dimension(SZI_(G),SZJB_(G)), intent(in) :: div_xx_dy !< y-derivative of horizontal divergence (d/dy(du/dx + dv/dy)) (m-1 s-1)
+ real, dimension(SZI_(G),SZJB_(G)), intent(inout) :: vort_xy_dx !< x-derivative of vertical vorticity (d/dx(dv/dx - du/dy)) (m-1 s-1)
+ real, dimension(SZIB_(G),SZJ_(G)), intent(inout) :: vort_xy_dy !< y-derivative of vertical vorticity (d/dy(dv/dx - du/dy)) (m-1 s-1)
+! real, dimension(SZI_(G),SZJ_(G)), intent(out) :: Leith_Kh_h !< Leith Laplacian viscosity at h-points (m2 s-1)
+! real, dimension(SZIB_(G),SZJB_(G)), intent(out) :: Leith_Kh_q !< Leith Laplacian viscosity at q-points (m2 s-1)
+! real, dimension(SZI_(G),SZJ_(G)), intent(out) :: Leith_Ah_h !< Leith bi-harmonic viscosity at h-points (m4 s-1)
+! real, dimension(SZIB_(G),SZJB_(G)), intent(out) :: Leith_Ah_q !< Leith bi-harmonic viscosity at q-points (m4 s-1)
+
+ ! Local variables
+! real, dimension(SZIB_(G),SZJB_(G)) :: vort_xy, & ! Vertical vorticity (dv/dx - du/dy) (s-1)
+! dudy, & ! Meridional shear of zonal velocity (s-1)
+! dvdx ! Zonal shear of meridional velocity (s-1)
+ real, dimension(SZI_(G),SZJB_(G)) :: &
+! vort_xy_dx, & ! x-derivative of vertical vorticity (d/dx(dv/dx - du/dy)) (m-1 s-1)
+! div_xx_dy, & ! y-derivative of horizontal divergence (d/dy(du/dx + dv/dy)) (m-1 s-1)
+ dslopey_dz, & ! z-derivative of y-slope at v-points (m-1)
+ h_at_v, & ! Thickness at v-points (m or kg m-2)
+ beta_v, & ! Beta at v-points (m-1 s-1)
+ grad_vort_mag_v, & ! mag. of vort. grad. at v-points (s-1)
+ grad_div_mag_v ! mag. of div. grad. at v-points (s-1)
+
+ real, dimension(SZIB_(G),SZJ_(G)) :: &
+! vort_xy_dy, & ! y-derivative of vertical vorticity (d/dy(dv/dx - du/dy)) (m-1 s-1)
+! div_xx_dx, & ! x-derivative of horizontal divergence (d/dx(du/dx + dv/dy)) (m-1 s-1)
+ dslopex_dz, & ! z-derivative of x-slope at u-points (m-1)
+ h_at_u, & ! Thickness at u-points (m or kg m-2)
+ beta_u, & ! Beta at u-points (m-1 s-1)
+ grad_vort_mag_u, & ! mag. of vort. grad. at u-points (s-1)
+ grad_div_mag_u ! mag. of div. grad. at u-points (s-1)
+! real, dimension(SZI_(G),SZJ_(G)) :: div_xx ! Estimate of horizontal divergence at h-points (s-1)
+! real :: mod_Leith, DY_dxBu, DX_dyBu, vert_vort_mag
+ real :: h_at_slope_above, h_at_slope_below, Ih, f
+ integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq,nz
+ real :: inv_PI3
+
+ is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
+ Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
+ nz = G%ke
+
+ inv_PI3 = 1.0/((4.0*atan(1.0))**3)
+
+ ! update halos
+ call pass_var(h, G%Domain)
+
+ if ((k > 1) .and. (k < nz)) then
+
+ ! Add in stretching term for the QG Leith vsicosity
+! if (CS%use_QG_Leith) then
+! do j=js-1,Jeq+1 ; do I=is-2,Ieq+1
+ do j=js-2,Jeq+2 ; do I=is-2,Ieq+1
+ h_at_slope_above = 2. * ( h(i,j,k-1) * h(i+1,j,k-1) ) * ( h(i,j,k) * h(i+1,j,k) ) / &
+ ( ( h(i,j,k-1) * h(i+1,j,k-1) ) * ( h(i,j,k) + h(i+1,j,k) ) &
+ + ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k-1) + h(i+1,j,k-1) ) + GV%H_subroundoff )
+ h_at_slope_below = 2. * ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k+1) * h(i+1,j,k+1) ) / &
+ ( ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k+1) + h(i+1,j,k+1) ) &
+ + ( h(i,j,k+1) * h(i+1,j,k+1) ) * ( h(i,j,k) + h(i+1,j,k) ) + GV%H_subroundoff )
+ Ih = 1./ ( ( h_at_slope_above + h_at_slope_below + GV%H_subroundoff ) * GV%H_to_m )
+ dslopex_dz(I,j) = 2. * ( CS%slope_x(i,j,k) - CS%slope_x(i,j,k+1) ) * Ih
+ h_at_u(I,j) = 2. * ( h_at_slope_above * h_at_slope_below ) * Ih
+ enddo ; enddo
+! do J=js-2,Jeq+1 ; do i=is-1,Ieq+1
+ do J=js-2,Jeq+1 ; do i=is-2,Ieq+2
+ h_at_slope_above = 2. * ( h(i,j,k-1) * h(i,j+1,k-1) ) * ( h(i,j,k) * h(i,j+1,k) ) / &
+ ( ( h(i,j,k-1) * h(i,j+1,k-1) ) * ( h(i,j,k) + h(i,j+1,k) ) &
+ + ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k-1) + h(i,j+1,k-1) ) + GV%H_subroundoff )
+ h_at_slope_below = 2. * ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k+1) * h(i,j+1,k+1) ) / &
+ ( ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k+1) + h(i,j+1,k+1) ) &
+ + ( h(i,j,k+1) * h(i,j+1,k+1) ) * ( h(i,j,k) + h(i,j+1,k) ) + GV%H_subroundoff )
+ Ih = 1./ ( ( h_at_slope_above + h_at_slope_below + GV%H_subroundoff ) * GV%H_to_m )
+ dslopey_dz(i,J) = 2. * ( CS%slope_y(i,j,k) - CS%slope_y(i,j,k+1) ) * Ih
+ h_at_v(i,J) = 2. * ( h_at_slope_above * h_at_slope_below ) * Ih
+ enddo ; enddo
+
+ do J=js-2,Jeq+1 ; do i=is-1,Ieq+1
+ f = 0.5 * ( G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J) )
+ vort_xy_dx(i,J) = vort_xy_dx(i,J) - f * &
+ ( ( h_at_u(I,j) * dslopex_dz(I,j) + h_at_u(I-1,j+1) * dslopex_dz(I-1,j+1) ) &
+ + ( h_at_u(I-1,j) * dslopex_dz(I-1,j) + h_at_u(I,j+1) * dslopex_dz(I,j+1) ) ) / &
+ ( ( h_at_u(I,j) + h_at_u(I-1,j+1) ) + ( h_at_u(I-1,j) + h_at_u(I,j+1) ) + GV%H_subroundoff)
+ enddo ; enddo
+
+ do j=js-1,Jeq+1 ; do I=is-2,Ieq+1
+ f = 0.5 * ( G%CoriolisBu(I,J) + G%CoriolisBu(I,J-1) )
+ vort_xy_dy(I,j) = vort_xy_dx(I,j) - f * &
+ ( ( h_at_v(i,J) * dslopey_dz(i,J) + h_at_v(i+1,J-1) * dslopey_dz(i+1,J-1) ) &
+ + ( h_at_v(i,J-1) * dslopey_dz(i,J-1) + h_at_v(i+1,J) * dslopey_dz(i+1,J) ) ) / &
+ ( ( h_at_v(i,J) + h_at_v(i+1,J-1) ) + ( h_at_v(i,J-1) + h_at_v(i+1,J) ) + GV%H_subroundoff)
+ enddo ; enddo
+ endif ! k > 1
+
+ call pass_vector(vort_xy_dy,vort_xy_dx,G%Domain)
+
+ if (CS%use_QG_Leith_GM) then
+ if (CS%use_beta_in_QG_Leith) then
+ do j=Jsq-1,Jeq+2 ; do I=is-2,Ieq+1
+ beta_u(I,j) = sqrt( (0.5*(G%dF_dx(i,j)+G%dF_dx(i+1,j))**2) + &
+ (0.5*(G%dF_dy(i,j)+G%dF_dy(i+1,j))**2) )
+ enddo ; enddo
+ do J=js-2,Jeq+1 ; do i=Isq-1,Ieq+2
+ beta_v(i,J) = sqrt( (0.5*(G%dF_dx(i,j)+G%dF_dx(i,j+1))**2) + &
+ (0.5*(G%dF_dy(i,j)+G%dF_dy(i,j+1))**2) )
+ enddo ; enddo
+ endif
+
+ do j=js-1,Jeq+1 ; do I=is-2,Ieq
+ grad_vort_mag_u(I,j) = SQRT(vort_xy_dy(I,j)**2 + (0.25*(vort_xy_dx(i,J) + vort_xy_dx(i+1,J) &
+ + vort_xy_dx(i,J-1) + vort_xy_dx(i+1,J-1)))**2)
+ grad_div_mag_u(I,j) = SQRT(div_xx_dx(I,j)**2 + (0.25*(div_xx_dy(i,J) + div_xx_dy(i+1,J) &
+ + div_xx_dy(i,J-1) + div_xx_dy(i+1,J-1)))**2)
+ if (CS%use_beta_in_QG_Leith) then
+ CS%KH_u_QG(I,j,k) = MIN(grad_vort_mag_u(I,j) + grad_div_mag_u(I,j), beta_u(I,j)*3) &
+ * CS%Laplac3_const_u(I,j) * inv_PI3
+ else
+ CS%KH_u_QG(I,j,k) = (grad_vort_mag_u(I,j) + grad_div_mag_u(I,j)) &
+ * CS%Laplac3_const_u(I,j) * inv_PI3
+ endif
+ enddo ; enddo
+
+ do J=js-2,Jeq ; do i=is-1,Ieq+1
+ grad_vort_mag_v(i,J) = SQRT(vort_xy_dx(i,J)**2 + (0.25*(vort_xy_dy(I,j) + vort_xy_dy(I-1,j) &
+ + vort_xy_dy(I,j+1) + vort_xy_dy(I-1,j+1)))**2)
+ grad_div_mag_v(i,J) = SQRT(div_xx_dy(i,J)**2 + (0.25*(div_xx_dx(I,j) + div_xx_dx(I-1,j) &
+ + div_xx_dx(I,j+1) + div_xx_dx(I-1,j+1)))**2)
+ if (CS%use_beta_in_QG_Leith) then
+ CS%KH_v_QG(i,J,k) = MIN(grad_vort_mag_v(i,J) + grad_div_mag_v(i,J), beta_v(i,J)*3) &
+ * CS%Laplac3_const_v(i,J) * inv_PI3
+ else
+ CS%KH_v_QG(i,J,k) = (grad_vort_mag_v(i,J) + grad_div_mag_v(i,J)) &
+ * CS%Laplac3_const_v(i,J) * inv_PI3
+ endif
+ enddo ; enddo
+ ! post diagnostics
+ if (CS%id_KH_v_QG > 0) call post_data(CS%id_KH_v_QG, CS%KH_v_QG, CS%diag)
+ if (CS%id_KH_u_QG > 0) call post_data(CS%id_KH_u_QG, CS%KH_u_QG, CS%diag)
+ endif
+
+end subroutine calc_QG_Leith_viscosity
+
!> Initializes the variables mixing coefficients container
subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
type(time_type), intent(in) :: Time !< Current model time
@@ -724,6 +891,9 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
! value is roughly (pi / (the age of the universe) )^2.
logical :: Gill_equatorial_Ld, use_FGNV_streamfn, use_MEKE, in_use
real :: MLE_front_length
+ real :: Leith_Lap_const ! The non-dimensional coefficient in the Leith viscosity
+ real :: grid_sp_u2, grid_sp_u3
+ real :: grid_sp_v2, grid_sp_v3 ! Intermediate quantities for Leith metrics
! This include declares and sets the variable "version".
#include "version_variable.h"
character(len=40) :: mdl = "MOM_lateral_mixing_coeffs" ! This module's name.
@@ -757,6 +927,9 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
"not used. If KHTR_SLOPE_CFF>0 or KhTh_Slope_Cff>0, "//&
"this is set to true regardless of what is in the "//&
"parameter file.", default=.false.)
+ call get_param(param_file, mdl, "USE_VISBECK", CS%use_Visbeck,&
+ "If true, use the Visbeck et al. (1997) formulation for \n"//&
+ "thickness diffusivity.", default=.false.)
call get_param(param_file, mdl, "RESOLN_SCALED_KH", CS%Resoln_scaled_Kh, &
"If true, the Laplacian lateral viscosity is scaled away "//&
"when the first baroclinic deformation radius is well "//&
@@ -828,6 +1001,8 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
in_use = .true.
allocate(CS%slope_x(IsdB:IedB,jsd:jed,G%ke+1)) ; CS%slope_x(:,:,:) = 0.0
allocate(CS%slope_y(isd:ied,JsdB:JedB,G%ke+1)) ; CS%slope_y(:,:,:) = 0.0
+ allocate(CS%N2_u(IsdB:IedB,jsd:jed,G%ke+1)) ; CS%N2_u(:,:,:) = 0.0
+ allocate(CS%N2_v(isd:ied,JsdB:JedB,G%ke+1)) ; CS%N2_v(:,:,:) = 0.0
call get_param(param_file, mdl, "KD_SMOOTH", CS%kappa_smooth, &
"A diapycnal diffusivity that is used to interpolate "//&
"more sensible values of T & S into thin layers.", &
@@ -884,6 +1059,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
'Depth average square of slope magnitude, S^2, at v-points, as used in Visbeck et al.', 's-2')
endif
+ oneOrTwo = 1.0
if (CS%Resoln_scaled_Kh .or. CS%Resoln_scaled_KhTh .or. CS%Resoln_scaled_KhTr) then
CS%calculate_Rd_dx = .true.
CS%calculate_res_fns = .true.
@@ -953,8 +1129,6 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
"is the more appropriate definition.", default=.false.)
if (Gill_equatorial_Ld) then
oneOrTwo = 2.0
- else
- oneOrTwo = 1.0
endif
do J=js-1,Jeq ; do I=is-1,Ieq
@@ -1020,6 +1194,49 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
call wave_speed_init(CS%wave_speed_CSp, use_ebt_mode=CS%Resoln_use_ebt, mono_N2_depth=N2_filter_depth)
endif
+ ! Leith parameters
+ call get_param(param_file, mdl, "USE_QG_LEITH_GM", CS%use_QG_Leith_GM, &
+ "If true, use the QG Leith viscosity as the GM coefficient.", &
+ default=.false.)
+
+ if (CS%Use_QG_Leith_GM) then
+ call get_param(param_file, mdl, "LEITH_LAP_CONST", Leith_Lap_const, &
+ "The nondimensional Laplacian Leith constant, \n"//&
+ "often set to 1.0", units="nondim", default=0.0)
+
+ call get_param(param_file, mdl, "USE_BETA_IN_LEITH", CS%use_beta_in_QG_Leith, &
+ "If true, include the beta term in the Leith nonlinear eddy viscosity.", &
+ default=.true.)
+
+ allocate(CS%Laplac3_const_u(IsdB:IedB,jsd:jed)) ; CS%Laplac3_const_u(:,:) = 0.0
+ allocate(CS%Laplac3_const_v(isd:ied,JsdB:JedB)) ; CS%Laplac3_const_v(:,:) = 0.0
+ allocate(CS%KH_u_QG(IsdB:IedB,jsd:jed,G%ke)) ; CS%KH_u_QG(:,:,:) = 0.0
+ allocate(CS%KH_v_QG(isd:ied,JsdB:JedB,G%ke)) ; CS%KH_v_QG(:,:,:) = 0.0
+ ! register diagnostics
+
+ CS%id_KH_u_QG = register_diag_field('ocean_model', 'KH_u_QG', diag%axesCuL, Time, &
+ 'Horizontal viscosity from Leith QG, at u-points', 'm2 s-1')
+ CS%id_KH_v_QG = register_diag_field('ocean_model', 'KH_v_QG', diag%axesCvL, Time, &
+ 'Horizontal viscosity from Leith QG, at v-points', 'm2 s-1')
+
+ do j=Jsq,Jeq+1 ; do I=is-1,Ieq
+ ! Static factors in the Leith schemes
+ grid_sp_u2 = G%dyCu(I,j)*G%dxCu(I,j)
+ grid_sp_u3 = grid_sp_u2*sqrt(grid_sp_u2)
+ CS%Laplac3_const_u(I,j) = Leith_Lap_const * grid_sp_u3
+ enddo ; enddo
+ do j=js-1,Jeq ; do I=Isq,Ieq+1
+ ! Static factors in the Leith schemes
+ grid_sp_v2 = G%dyCv(i,J)*G%dxCu(i,J)
+ grid_sp_v3 = grid_sp_v2*sqrt(grid_sp_v2)
+ CS%Laplac3_const_v(i,J) = Leith_Lap_const * grid_sp_v3
+ enddo ; enddo
+
+ if (.not. CS%use_stored_slopes) call MOM_error(FATAL, &
+ "MOM_lateral_mixing_coeffs.F90, VarMix_init:"//&
+ "USE_STORED_SLOPES must be True when using QG Leith.")
+ endif
+
! If nothing is being stored in this class then deallocate
if (in_use) then
CS%use_variable_mixing = .true.
diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90
index 4d75494b72..da18462da6 100644
--- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90
+++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90
@@ -7,7 +7,8 @@ module MOM_thickness_diffuse
use MOM_diag_mediator, only : post_data, query_averaging_enabled, diag_ctrl
use MOM_diag_mediator, only : register_diag_field, safe_alloc_ptr, time_type
use MOM_diag_mediator, only : diag_update_remap_grids
-use MOM_error_handler, only : MOM_error, FATAL, WARNING
+use MOM_domains, only : pass_var, CORNER, pass_vector
+use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe
use MOM_EOS, only : calculate_density, calculate_density_derivs
use MOM_file_parser, only : get_param, log_version, param_file_type
use MOM_grid, only : ocean_grid_type
@@ -23,7 +24,7 @@ module MOM_thickness_diffuse
#include
public thickness_diffuse, thickness_diffuse_init, thickness_diffuse_end
-public vert_fill_TS
+public vert_fill_TS, thickness_diffuse_get_KH
! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
@@ -58,11 +59,26 @@ module MOM_thickness_diffuse
!! longer than DT, or 0 (the default) to use DT.
integer :: nkml !< number of layers within mixed layer
logical :: debug !< write verbose checksums for debugging purposes
+ logical :: use_GME_thickness_diffuse !< If true, passes GM coefficients to MOM_hor_visc for use
+ !! with GME closure.
+ logical :: MEKE_GEOMETRIC !< If true, uses the GM coefficient formulation from the GEOMETRIC
+ !! framework (Marshall et al., 2012)
+ real :: MEKE_GEOMETRIC_alpha!< The nondimensional coefficient governing the efficiency of
+ !! the GEOMETRIC thickness difussion [nondim]
+ real :: MEKE_GEOMETRIC_epsilon !< Minimum Eady growth rate for the GEOMETRIC thickness
+ !! diffusivity [s-1].
+ logical :: Use_KH_in_MEKE !< If true, uses the thickness diffusivity calculated here to diffuse MEKE.
+ logical :: GM_src_alt !< If true, use the GM energy conversion form S^2*N^2*kappa rather
+ !! than the streamfunction for the GM source term.
type(diag_ctrl), pointer :: diag => NULL() !< structure used to regulate timing of diagnostics
real, pointer :: GMwork(:,:) => NULL() !< Work by thickness diffusivity [W m-2]
real, pointer :: diagSlopeX(:,:,:) => NULL() !< Diagnostic: zonal neutral slope [nondim]
real, pointer :: diagSlopeY(:,:,:) => NULL() !< Diagnostic: zonal neutral slope [nondim]
+ real, dimension(:,:,:), pointer :: &
+ KH_u_GME => NULL(), & !< interface height diffusivities in u-columns (m2 s-1)
+ KH_v_GME => NULL() !< interface height diffusivities in v-columns (m2 s-1)
+
!>@{
!! Diagnostic identifier
integer :: id_uhGM = -1, id_vhGM = -1, id_GMwork = -1
@@ -123,7 +139,8 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp
real :: h_neglect ! A thickness that is so small it is usually lost
! in roundoff and can be neglected [H ~> m or kg m-2].
real, dimension(:,:), pointer :: cg1 => null() !< Wave speed [m s-1]
- logical :: use_VarMix, Resoln_scaled, use_stored_slopes, khth_use_ebt_struct
+ logical :: use_VarMix, Resoln_scaled, use_stored_slopes, khth_use_ebt_struct, use_Visbeck
+ logical :: use_QG_Leith
integer :: i, j, k, is, ie, js, je, nz
real :: hu(SZI_(G), SZJ_(G)) ! u-thickness [H ~> m or kg m-2]
real :: hv(SZI_(G), SZJ_(G)) ! v-thickness [H ~> m or kg m-2]
@@ -146,17 +163,21 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp
endif
use_VarMix = .false. ; Resoln_scaled = .false. ; use_stored_slopes = .false.
- khth_use_ebt_struct = .false.
+ khth_use_ebt_struct = .false. ; use_Visbeck = .false. ; use_QG_Leith = .false.
+
if (associated(VarMix)) then
use_VarMix = VarMix%use_variable_mixing .and. (CS%KHTH_Slope_Cff > 0.)
Resoln_scaled = VarMix%Resoln_scaled_KhTh
use_stored_slopes = VarMix%use_stored_slopes
khth_use_ebt_struct = VarMix%khth_use_ebt_struct
+ use_Visbeck = VarMix%use_Visbeck
+ use_QG_Leith = VarMix%use_QG_Leith_GM
if (associated(VarMix%cg1)) cg1 => VarMix%cg1
else
cg1 => null()
endif
+
!$OMP parallel do default(none) shared(is,ie,js,je,KH_u_CFL,dt,G,CS)
do j=js,je ; do I=is-1,ie
KH_u_CFL(I,j) = (0.25*CS%max_Khth_CFL) / &
@@ -183,16 +204,25 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp
if (use_VarMix) then
!$OMP do
- do j=js,je ; do I=is-1,ie
- Khth_Loc_u(I,j) = Khth_Loc_u(I,j) + CS%KHTH_Slope_Cff*VarMix%L2u(I,j)*VarMix%SN_u(I,j)
- enddo ; enddo
+ if (use_Visbeck) then
+ do j=js,je ; do I=is-1,ie
+ Khth_Loc_u(I,j) = Khth_Loc_u(I,j) + CS%KHTH_Slope_Cff*VarMix%L2u(I,j)*VarMix%SN_u(I,j)
+ enddo ; enddo
+ endif
endif
if (associated(MEKE)) then ; if (associated(MEKE%Kh)) then
!$OMP do
- do j=js,je ; do I=is-1,ie
- Khth_Loc_u(I,j) = Khth_Loc_u(I,j) + MEKE%KhTh_fac*sqrt(MEKE%Kh(i,j)*MEKE%Kh(i+1,j))
- enddo ; enddo
+ if (CS%MEKE_GEOMETRIC) then
+ do j=js,je ; do I=is-1,ie
+ Khth_Loc_u(I,j) = Khth_Loc_u(I,j) + G%mask2dCu(I,j) * CS%MEKE_GEOMETRIC_alpha * 0.5*(MEKE%MEKE(i,j)+MEKE%MEKE(i+1,j)) / &
+ (VarMix%SN_u(I,j) + CS%MEKE_GEOMETRIC_epsilon)
+ enddo ; enddo
+ else
+ do j=js,je ; do I=is-1,ie
+ Khth_Loc_u(I,j) = Khth_Loc_u(I,j) + MEKE%KhTh_fac*sqrt(MEKE%Kh(i,j)*MEKE%Kh(i+1,j))
+ enddo ; enddo
+ endif
endif ; endif
if (Resoln_scaled) then
@@ -230,6 +260,22 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp
enddo ; enddo ; enddo
endif
+ if (use_VarMix) then
+!$OMP do
+ if (use_QG_Leith) then
+ do k=1,nz ; do j=js,je ; do I=is-1,ie
+ KH_u(I,j,k) = VarMix%KH_u_QG(I,j,k)
+ enddo ; enddo ; enddo
+ endif
+ endif
+
+!$OMP do
+ if (CS%use_GME_thickness_diffuse) then
+ do k=1,nz+1 ; do j=js,je ; do I=is-1,ie
+ CS%KH_u_GME(I,j,k) = KH_u(I,j,k)
+ enddo ; enddo ; enddo
+ endif
+
!$OMP do
do J=js-1,je ; do i=is,ie
Khth_Loc(i,j) = CS%Khth
@@ -237,15 +283,24 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp
if (use_VarMix) then
!$OMP do
- do J=js-1,je ; do i=is,ie
- Khth_Loc(i,j) = Khth_Loc(i,j) + CS%KHTH_Slope_Cff*VarMix%L2v(i,J)*VarMix%SN_v(i,J)
- enddo ; enddo
+ if (use_Visbeck) then
+ do J=js-1,je ; do i=is,ie
+ Khth_Loc(i,j) = Khth_Loc(i,j) + CS%KHTH_Slope_Cff*VarMix%L2v(i,J)*VarMix%SN_v(i,J)
+ enddo ; enddo
+ endif
endif
if (associated(MEKE)) then ; if (associated(MEKE%Kh)) then
!$OMP do
- do J=js-1,je ; do i=is,ie
- Khth_Loc(i,j) = Khth_Loc(i,j) + MEKE%KhTh_fac*sqrt(MEKE%Kh(i,j)*MEKE%Kh(i,j+1))
- enddo ; enddo
+ if (CS%MEKE_GEOMETRIC) then
+ do j=js-1,je ; do I=is,ie
+ Khth_Loc(I,j) = Khth_Loc(I,j) + G%mask2dCv(i,J) * CS%MEKE_GEOMETRIC_alpha * 0.5*(MEKE%MEKE(i,j)+MEKE%MEKE(i,j+1)) / &
+ (VarMix%SN_v(i,J) + CS%MEKE_GEOMETRIC_epsilon)
+ enddo ; enddo
+ else
+ do J=js-1,je ; do i=is,ie
+ Khth_Loc(i,j) = Khth_Loc(i,j) + MEKE%KhTh_fac*sqrt(MEKE%Kh(i,j)*MEKE%Kh(i,j+1))
+ enddo ; enddo
+ endif
endif ; endif
if (Resoln_scaled) then
@@ -273,6 +328,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp
KH_v(i,J,1) = min(KH_v_CFL(i,J), Khth_Loc(i,j))
enddo ; enddo
endif
+
if (khth_use_ebt_struct) then
!$OMP do
do K=2,nz+1 ; do J=js-1,je ; do i=is,ie
@@ -284,6 +340,35 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp
KH_v(i,J,K) = KH_v(i,J,1)
enddo ; enddo ; enddo
endif
+
+ if (use_VarMix) then
+!$OMP do
+ if (use_QG_Leith) then
+ do k=1,nz ; do J=js-1,je ; do i=is,ie
+ KH_v(i,J,k) = VarMix%KH_v_QG(i,J,k)
+ enddo ; enddo ; enddo
+ endif
+ endif
+
+!$OMP do
+ if (CS%use_GME_thickness_diffuse) then
+ do k=1,nz+1 ; do j=js-1,je ; do I=is,ie
+ CS%KH_v_GME(I,j,k) = KH_v(I,j,k)
+ enddo ; enddo ; enddo
+ endif
+
+ if (associated(MEKE)) then ; if (associated(MEKE%Kh)) then
+!$OMP do
+ if (CS%MEKE_GEOMETRIC) then
+ do j=js,je ; do I=is,ie
+ MEKE%Kh(i,j) = CS%MEKE_GEOMETRIC_alpha * MEKE%MEKE(i,j) / &
+ (0.25*(VarMix%SN_u(I,j)+VarMix%SN_u(I-1,j)+VarMix%SN_v(i,J)+VarMix%SN_v(i,J-1)) + &
+ CS%MEKE_GEOMETRIC_epsilon)
+ enddo ; enddo
+ endif
+ endif ; endif
+
+
!$OMP do
do K=1,nz+1 ; do j=js,je ; do I=is-1,ie ; int_slope_u(I,j,K) = 0.0 ; enddo ; enddo ; enddo
!$OMP do
@@ -343,7 +428,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp
! in the case where KH_u and KH_v are depth independent. Otherwise,
! if use thickness weighted average, the variations of thickness with
! depth will place a spurious depth dependence to the diagnosed KH_t.
- if (CS%id_KH_t > 0 .or. CS%id_KH_t1 > 0) then
+ if (CS%id_KH_t > 0 .or. CS%id_KH_t1 > 0 .or. CS%Use_KH_in_MEKE) then
do k=1,nz
! thicknesses across u and v faces, converted to 0/1 mask
! layer average of the interface diffusivities KH_u and KH_v
@@ -364,6 +449,20 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp
/ (hu(I-1,j)+hu(I,j)+hv(i,J-1)+hv(i,J)+h_neglect)
enddo ; enddo
enddo
+
+ if (CS%Use_KH_in_MEKE) then
+ MEKE%Kh_diff(:,:) = 0.0
+ do k=1,nz
+ do j=js,je ; do i=is,ie
+ MEKE%Kh_diff(i,j) = MEKE%Kh_diff(i,j) + KH_t(i,j,k) * h(i,j,k)
+ enddo; enddo
+ enddo
+
+ do j=js,je ; do i=is,ie
+ MEKE%Kh_diff(i,j) = MEKE%Kh_diff(i,j) / MAX(1.0,G%bathyT(i,j))
+ enddo ; enddo
+ endif
+
if (CS%id_KH_t > 0) call post_data(CS%id_KH_t, KH_t, CS%diag)
if (CS%id_KH_t1 > 0) call post_data(CS%id_KH_t1, KH_t(:,:,1), CS%diag)
endif
@@ -435,7 +534,6 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
!! density gradients.
real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), optional, intent(in) :: slope_x !< Isopycnal slope at u-points
real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), optional, intent(in) :: slope_y !< Isopycnal slope at v-points
-
! Local variables
real, dimension(SZI_(G), SZJ_(G), SZK_(G)) :: &
T, & ! The temperature (or density) [degC], with the values in
@@ -448,6 +546,14 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
! by dt [H m2 s-1 ~> m3 s-1 or kg s-1].
h_frac ! The fraction of the mass in the column above the bottom
! interface of a layer that is within a layer [nondim]. 0 m3 s-1 or kg s-1].
@@ -470,6 +576,8 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
real :: Work_u(SZIB_(G), SZJ_(G)) ! The work being done by the thickness
real :: Work_v(SZI_(G), SZJB_(G)) ! diffusion integrated over a cell [W].
real :: Work_h ! The work averaged over an h-cell [W m-2].
+ real :: PE_release_h ! The amount of potential energy released by GM, averaged over an h-cell [m3 s-3].
+ ! The calculation is equal to h * S^2 * N^2 * kappa_GM.
real :: I4dt ! 1 / 4 dt [s-1].
real :: drdiA, drdiB ! Along layer zonal- and meridional- potential density
real :: drdjA, drdjB ! gradients in the layers above (A) and below(B) the
@@ -544,6 +652,11 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
nk_linear = max(GV%nkml, 1)
+ Slope_x_PE(:,:,:) = 0.0
+ Slope_y_PE(:,:,:) = 0.0
+ hN2_x_PE(:,:,:) = 0.0
+ hN2_y_PE(:,:,:) = 0.0
+
find_work = .false.
if (associated(MEKE)) find_work = associated(MEKE%GM_src)
find_work = (associated(CS%GMwork) .or. find_work)
@@ -650,7 +763,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
if (k > nk_linear) then
if (use_EOS) then
- if (CS%use_FGNV_streamfn .or. .not.present_slope_x) then
+ if (CS%use_FGNV_streamfn .or. find_work .or. .not.present_slope_x) then
hg2L = h(i,j,k-1)*h(i,j,k) + h_neglect2
hg2R = h(i+1,j,k-1)*h(i+1,j,k) + h_neglect2
haL = 0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect
@@ -708,6 +821,9 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
int_slope_u(I,j,K) * US%Z_to_m*((e(i+1,j,K)-e(i,j,K)) * G%IdxCu(I,j))
slope2_Ratio_u(I,K) = (1.0 - int_slope_u(I,j,K)) * slope2_Ratio_u(I,K)
endif
+
+ Slope_x_PE(I,j,k) = MIN(Slope,CS%slope_max)
+ hN2_x_PE(I,j,k) = hN2_u(I,K) * US%m_to_Z
if (CS%id_slope_x > 0) CS%diagSlopeX(I,j,k) = Slope
! Estimate the streamfunction at each interface [m3 s-1].
@@ -896,7 +1012,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
if (k > nk_linear) then
if (use_EOS) then
- if (CS%use_FGNV_streamfn .or. .not.present_slope_y) then
+ if (CS%use_FGNV_streamfn .or. find_work .or. .not. present_slope_y) then
hg2L = h(i,j,k-1)*h(i,j,k) + h_neglect2
hg2R = h(i,j+1,k-1)*h(i,j+1,k) + h_neglect2
haL = 0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect
@@ -954,6 +1070,9 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
int_slope_v(i,J,K) * US%Z_to_m*((e(i,j+1,K)-e(i,j,K)) * G%IdyCv(i,J))
slope2_Ratio_v(i,K) = (1.0 - int_slope_v(i,J,K)) * slope2_Ratio_v(i,K)
endif
+
+ Slope_y_PE(i,J,k) = MIN(Slope,CS%slope_max)
+ hN2_y_PE(i,J,k) = hN2_v(i,K) * US%m_to_Z
if (CS%id_slope_y > 0) CS%diagSlopeY(I,j,k) = Slope
! Estimate the streamfunction at each interface [m3 s-1].
@@ -1142,14 +1261,25 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
enddo
endif
+
+ !if (find_work) then ; do j=js,je ; do i=is,ie ; do k=nz,1,-1
if (find_work) then ; do j=js,je ; do i=is,ie
! Note that the units of Work_v and Work_u are W, while Work_h is W m-2.
Work_h = 0.5 * G%IareaT(i,j) * &
((Work_u(I-1,j) + Work_u(I,j)) + (Work_v(i,J-1) + Work_v(i,J)))
+ PE_release_h = -0.25*(Kh_u(I,j,k)*(Slope_x_PE(I,j,k)**2) * hN2_x_PE(I,j,k) + &
+ Kh_u(I-1,j,k)*(Slope_x_PE(I-1,j,k)**2) * hN2_x_PE(I-1,j,k) + &
+ Kh_v(i,J,k)*(Slope_y_PE(i,J,k)**2) * hN2_y_PE(i,J,k) + &
+ Kh_v(i,J-1,k)*(Slope_y_PE(i,J-1,k)**2) * hN2_y_PE(i,J-1,k))
if (associated(CS%GMwork)) CS%GMwork(i,j) = Work_h
if (associated(MEKE)) then ; if (associated(MEKE%GM_src)) then
- MEKE%GM_src(i,j) = MEKE%GM_src(i,j) + Work_h
+ if (CS%GM_src_alt) then
+ MEKE%GM_src(i,j) = MEKE%GM_src(i,j) + PE_release_h
+ else
+ MEKE%GM_src(i,j) = MEKE%GM_src(i,j) + Work_h
+ endif
endif ; endif
+ !enddo ; enddo ; enddo ; endif
enddo ; enddo ; endif
if (CS%id_slope_x > 0) call post_data(CS%id_slope_x, CS%diagSlopeX, CS%diag)
@@ -1739,6 +1869,11 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS)
"marginally unstable value in a pure layered model, but "//&
"much smaller numbers (e.g. 0.1) seem to work better for "//&
"ALE-based models.", units = "nondimensional", default=0.8)
+
+! call get_param(param_file, mdl, "USE_QG_LEITH_GM", CS%QG_Leith_GM, &
+! "If true, use the QG Leith viscosity as the GM coefficient.", &
+! default=.false.)
+
if (CS%max_Khth_CFL < 0.0) CS%max_Khth_CFL = 0.0
call get_param(param_file, mdl, "DETANGLE_INTERFACES", CS%detangle_interfaces, &
"If defined add 3-d structured enhanced interface height "//&
@@ -1783,6 +1918,35 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS)
"If true, write out verbose debugging data.", &
default=.false., debuggingParam=.true.)
+ call get_param(param_file, mdl, "MEKE_GM_SRC_ALT", CS%GM_src_alt, &
+ "If true, use the GM energy conversion form S^2*N^2*kappa rather \n"//&
+ "than the streamfunction for the GM source term.", default=.false.)
+ call get_param(param_file, mdl, "MEKE_GEOMETRIC", CS%MEKE_GEOMETRIC, &
+ "If true, uses the GM coefficient formulation \n"//&
+ "from the GEOMETRIC framework (Marshall et al., 2012).", default=.false.)
+ if (CS%MEKE_GEOMETRIC) then
+
+ call get_param(param_file, mdl, "MEKE_GEOMETRIC_EPSILON", CS%MEKE_GEOMETRIC_epsilon, &
+ "Minimum Eady growth rate used in the calculation of \n"//&
+ "GEOMETRIC thickness diffusivity.", units="s-1", default=1.0e-7)
+
+ call get_param(param_file, mdl, "MEKE_GEOMETRIC_ALPHA", CS%MEKE_GEOMETRIC_alpha, &
+ "The nondimensional coefficient governing the efficiency of the GEOMETRIC \n"//&
+ "thickness diffusion.", units="nondim", default=0.05)
+ endif
+
+ call get_param(param_file, mdl, "USE_KH_IN_MEKE", CS%Use_KH_in_MEKE, &
+ "If true, uses the thickness diffusivity calculated here to diffuse \n"//&
+ "MEKE.", default=.false.)
+
+ call get_param(param_file, mdl, "USE_GME", CS%use_GME_thickness_diffuse, &
+ "If true, use the GM+E backscatter scheme in association \n"//&
+ "with the Gent and McWilliams parameterization.", default=.false.)
+
+ if (CS%use_GME_thickness_diffuse) then
+ call safe_alloc_ptr(CS%KH_u_GME,G%IsdB,G%IedB,G%jsd,G%jed,G%ke+1)
+ call safe_alloc_ptr(CS%KH_v_GME,G%isd,G%ied,G%JsdB,G%JedB,G%ke+1)
+ endif
if (GV%Boussinesq) then ; flux_to_kg_per_s = GV%Rho0
else ; flux_to_kg_per_s = 1. ; endif
@@ -1840,6 +2004,28 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS)
end subroutine thickness_diffuse_init
+!> Copies ubtav and vbtav from private type into arrays
+subroutine thickness_diffuse_get_KH(CS, KH_u_GME, KH_v_GME, G)
+ type(thickness_diffuse_CS), pointer :: CS !< Control structure for
+ !! this module
+ type(ocean_grid_type), intent(in) :: G !< Grid structure
+ real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: KH_u_GME!< interface height
+ !! diffusivities in u-columns [m2 s-1]
+ real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), intent(inout) :: KH_v_GME!< interface height
+ !! diffusivities in v-columns [m2 s-1]
+ ! Local variables
+ integer :: i,j,k
+
+ do k=1,G%ke+1 ; do j = G%jsc, G%jec ; do I = G%isc-1, G%iec
+ KH_u_GME(I,j,k) = CS%KH_u_GME(I,j,k)
+ enddo ; enddo ; enddo
+
+ do k=1,G%ke+1 ; do J = G%jsc-1, G%jec ; do i = G%isc, G%iec
+ KH_v_GME(i,J,k) = CS%KH_v_GME(i,J,k)
+ enddo ; enddo ; enddo
+
+end subroutine thickness_diffuse_get_KH
+
!> Deallocate the thickness diffusion control structure
subroutine thickness_diffuse_end(CS)
type(thickness_diffuse_CS), pointer :: CS !< Control structure for thickness diffusion