From 4dd011148127f7da24c7d4b9ccfc7e87d0d65084 Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Wed, 16 May 2018 18:08:27 -0600 Subject: [PATCH 1/6] add original OBLdepth diag --- src/parameterizations/vertical/MOM_KPP.F90 | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/parameterizations/vertical/MOM_KPP.F90 b/src/parameterizations/vertical/MOM_KPP.F90 index f98185685a..185be4f0b4 100644 --- a/src/parameterizations/vertical/MOM_KPP.F90 +++ b/src/parameterizations/vertical/MOM_KPP.F90 @@ -134,10 +134,12 @@ module MOM_KPP integer :: id_NLT_temp_budget = -1 integer :: id_NLT_saln_budget = -1 integer :: id_EnhK = -1, id_EnhW = -1, id_EnhVt2 = -1 + integer :: id_OBLdepth_original = -1 ! Diagnostics arrays real, allocatable, dimension(:,:) :: OBLdepth !< Depth (positive) of OBL (m) + real, allocatable, dimension(:,:) :: OBLdepth_original !< Depth (positive) of OBL (m) without smoothing real, allocatable, dimension(:,:) :: kOBL !< Level (+fraction) of OBL extent real, allocatable, dimension(:,:) :: OBLdepthprev !< previous Depth (positive) of OBL (m) real, allocatable, dimension(:,:,:) :: dRho !< Bulk difference in density (kg/m3) @@ -437,6 +439,10 @@ logical function KPP_init(paramFile, G, diag, Time, CS, passive, Waves) ! CMOR names are placeholders; must be modified by time period ! for CMOR compliance. Diag manager will be used for omlmax and ! omldamax. + CS%id_OBLdepth_original = register_diag_field('ocean_model', 'KPP_OBLdepth_original', diag%axesT1, Time, & + 'Thickness of the surface Ocean Boundary Layer without smoothing calculated by [CVMix] KPP', 'meter', & + cmor_field_name='oml', cmor_long_name='ocean_mixed_layer_thickness_defined_by_mixing_scheme', & + cmor_units='m', cmor_standard_name='Ocean Mixed Layer Thickness Defined by Mixing Scheme') CS%id_BulkDrho = register_diag_field('ocean_model', 'KPP_BulkDrho', diag%axesTL, Time, & 'Bulk difference in density used in Bulk Richardson number, as used by [CVMix] KPP', 'kg/m3') CS%id_BulkUz2 = register_diag_field('ocean_model', 'KPP_BulkUz2', diag%axesTL, Time, & @@ -498,6 +504,7 @@ logical function KPP_init(paramFile, G, diag, Time, CS, passive, Waves) CS%OBLdepth(:,:) = 0. allocate( CS%kOBL( SZI_(G), SZJ_(G) ) ) CS%kOBL(:,:) = 0. + if (CS%id_OBLdepth_original > 0) allocate( CS%OBLdepth_original( SZI_(G), SZJ_(G) ) ) allocate( CS%OBLdepthprev( SZI_(G), SZJ_(G) ) );CS%OBLdepthprev(:,:)=0.0 if (CS%id_BulkDrho > 0) allocate( CS%dRho( SZI_(G), SZJ_(G), SZK_(G) ) ) @@ -1063,6 +1070,7 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & ! send diagnostics to post_data if (CS%id_OBLdepth > 0) call post_data(CS%id_OBLdepth, CS%OBLdepth, CS%diag) + if (CS%id_OBLdepth_original > 0) call post_data(CS%id_OBLdepth_original,CS%OBLdepth_original,CS%diag) if (CS%id_BulkDrho > 0) call post_data(CS%id_BulkDrho, CS%dRho, CS%diag) if (CS%id_BulkUz2 > 0) call post_data(CS%id_BulkUz2, CS%Uz2, CS%diag) if (CS%id_BulkRi > 0) call post_data(CS%id_BulkRi, CS%BulkRi, CS%diag) @@ -1516,6 +1524,7 @@ subroutine KPP_smooth_BLD(CS,G,GV,h) call pass_var(CS%OBLdepth, G%Domain) OBLdepth_original = CS%OBLdepth + CS%OBLdepth_original = OBLdepth_original ! apply smoothing on OBL depth do j = G%jsc, G%jec From cfc6742025f71fdadef9c9ebc2effdfced6504cc Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Thu, 17 May 2018 15:07:12 -0600 Subject: [PATCH 2/6] prevent OBL depths deeper than the bathymetric depth --- src/parameterizations/vertical/MOM_KPP.F90 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/parameterizations/vertical/MOM_KPP.F90 b/src/parameterizations/vertical/MOM_KPP.F90 index 185be4f0b4..da59d487cc 100644 --- a/src/parameterizations/vertical/MOM_KPP.F90 +++ b/src/parameterizations/vertical/MOM_KPP.F90 @@ -1548,6 +1548,9 @@ subroutine KPP_smooth_BLD(CS,G,GV,h) enddo enddo + ! prevent OBL depths deeper than the bathymetric depth + where (CS%OBLdepth > G%bathyT) CS%OBLdepth = G%bathyT + ! Update kOBL for smoothed OBL depths do j = G%jsc, G%jec do i = G%isc, G%iec From 8c83e21f4679674341cdf7c6224922cdd8425284 Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Thu, 17 May 2018 17:04:53 -0600 Subject: [PATCH 3/6] add viscosities due to tidal mixing --- .../vertical/MOM_set_diffusivity.F90 | 2 +- .../vertical/MOM_tidal_mixing.F90 | 24 +++++++++++++++---- 2 files changed, 20 insertions(+), 6 deletions(-) diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 9906083597..29cd8c7ca7 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -441,7 +441,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & ! Add the Nikurashin and / or tidal bottom-driven mixing call calculate_tidal_mixing(h, N2_bot, j, TKE_to_Kd, maxTKE, G, GV, CS%tm_csp, & - N2_lay, N2_int, Kd, Kd_int, CS%Kd_max) + N2_lay, N2_int, Kd, Kd_int, CS%Kd_max, visc%Kv_slow) ! This adds the diffusion sustained by the energy extracted from the flow ! by the bottom drag. diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90 index 9c60a05451..b659e9149a 100644 --- a/src/parameterizations/vertical/MOM_tidal_mixing.F90 +++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90 @@ -642,7 +642,7 @@ end function tidal_mixing_init !! tidal dissipation and to add the effect of internal-tide-driven mixing to the layer or interface !! diffusivities. subroutine calculate_tidal_mixing(h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, CS, & - N2_lay, N2_int, Kd, Kd_int, Kd_max) + N2_lay, N2_int, Kd, Kd_int, Kd_max, Kv) 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(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses, in H (usually m or kg m-2) @@ -655,10 +655,12 @@ subroutine calculate_tidal_mixing(h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, CS, & real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: Kd real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), optional, intent(inout) :: Kd_int real, intent(inout) :: Kd_max + real, dimension(:,:,:), pointer :: Kv !< The "slow" vertical viscosity at each interface + !! (not layer!) in m2 s-1. if (CS%Int_tide_dissipation .or. CS%Lee_wave_dissipation .or. CS%Lowmode_itidal_dissipation) then if (CS%use_CVMix_tidal) then - call calculate_CVMix_tidal(h, j, G, GV, CS, N2_int, Kd) + call calculate_CVMix_tidal(h, j, G, GV, CS, N2_int, Kd, Kv) else call add_int_tide_diffusivity(h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, CS, & N2_lay, Kd, Kd_int, Kd_max) @@ -669,7 +671,7 @@ subroutine calculate_tidal_mixing(h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, CS, & !> Calls the CVMix routines to compute tidal dissipation and to add the effect of internal-tide-driven !! mixing to the interface diffusivities. -subroutine calculate_CVMix_tidal(h, j, G, GV, CS, N2_int, Kd) +subroutine calculate_CVMix_tidal(h, j, G, GV, CS, N2_int, Kd, Kv) integer, intent(in) :: j type(ocean_grid_type), intent(in) :: G !< Grid structure. type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure @@ -677,6 +679,8 @@ subroutine calculate_CVMix_tidal(h, j, G, GV, CS, N2_int, Kd) real, dimension(SZI_(G),SZK_(G)+1), intent(in) :: N2_int real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses, in H (usually m or kg m-2). real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: Kd + real, dimension(:,:,:), pointer :: Kv !< The "slow" vertical viscosity at each interface + !! (not layer!) in m2 s-1. ! local real, dimension(SZK_(G)+1) :: Kd_tidal !< tidal diffusivity [m2/s] @@ -744,9 +748,14 @@ subroutine calculate_CVMix_tidal(h, j, G, GV, CS, N2_int, Kd) CVMix_params = CS%CVMix_glb_params, & CVMix_tidal_params_user = CS%CVMix_tidal_params) + ! Update diffusivity do k=1,G%ke Kd(i,j,k) = Kd(i,j,k) + 0.5*(Kd_tidal(k) + Kd_tidal(k+1) ) - !TODO: Kv(i,j,k) = ???????????? + enddo + + ! Update viscosity + do k=1,G%ke+1 + Kv(i,j,k) = Kv(i,j,k) + Kv_tidal(k) enddo ! diagnostics @@ -836,9 +845,14 @@ subroutine calculate_CVMix_tidal(h, j, G, GV, CS, N2_int, Kd) CVmix_params = CS%CVMix_glb_params, & CVmix_tidal_params_user = CS%CVMix_tidal_params) + ! Update diffusivity do k=1,G%ke Kd(i,j,k) = Kd(i,j,k) + 0.5*(Kd_tidal(k) + Kd_tidal(k+1) ) - !TODO: Kv(i,j,k) = ???????????? + enddo + + ! Update viscosity + do k=1,G%ke+1 + Kv(i,j,k) = Kv(i,j,k) + Kv_tidal(k) enddo ! diagnostics From e12d9525fff79e46cff9a42118508769fca66f08 Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Thu, 31 May 2018 13:14:14 -0600 Subject: [PATCH 4/6] fix a typo in ocn_comp_mct.F90 --- config_src/mct_driver/ocn_comp_mct.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/config_src/mct_driver/ocn_comp_mct.F90 b/config_src/mct_driver/ocn_comp_mct.F90 index c3967caf6d..aece6abebc 100644 --- a/config_src/mct_driver/ocn_comp_mct.F90 +++ b/config_src/mct_driver/ocn_comp_mct.F90 @@ -1772,7 +1772,7 @@ subroutine update_ocean_model(OS, Ocean_sfc, time_start_update, & call set_net_mass_forcing(OS%fluxes, OS%forces, OS%grid) if (OS%nstep==0) then - call finish_MOM_initialization(OS%Time, OS%dirs, OS%MOM_CSp, S%restart_CSp) + call finish_MOM_initialization(OS%Time, OS%dirs, OS%MOM_CSp, OS%restart_CSp) endif call disable_averaging(OS%diag) From b2fac04eb9924032363ba79876c1e995b77c3987 Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Thu, 31 May 2018 18:43:55 -0600 Subject: [PATCH 5/6] add multi-smoothing and deepening-only smoothing options --- src/parameterizations/vertical/MOM_KPP.F90 | 78 +++++++++++++--------- 1 file changed, 46 insertions(+), 32 deletions(-) diff --git a/src/parameterizations/vertical/MOM_KPP.F90 b/src/parameterizations/vertical/MOM_KPP.F90 index da59d487cc..aee9b28b7e 100644 --- a/src/parameterizations/vertical/MOM_KPP.F90 +++ b/src/parameterizations/vertical/MOM_KPP.F90 @@ -88,7 +88,8 @@ module MOM_KPP character(len=30) :: MatchTechnique !< Method used in CVMix for setting diffusivity and NLT profile functions integer :: NLT_shape !< MOM6 over-ride of CVMix NLT shape function logical :: applyNonLocalTrans !< If True, apply non-local transport to heat and scalars - logical :: smoothBLD !< If True, apply a 1-1-4-1-1 Laplacian filter one time on HBLT. + integer :: n_smooth !< Number of times smoothing operator is applied on OBLdepth. + logical :: deepen_only !< If true, apply OBLdepth smoothing at a cell only if the OBLdepth gets deeper. logical :: KPPzeroDiffusivity !< If True, will set diffusivity and viscosity from KPP to zero !! for testing purposes. logical :: KPPisAdditive !< If True, will add KPP diffusivity to initial diffusivity. @@ -216,10 +217,16 @@ logical function KPP_init(paramFile, G, diag, Time, CS, passive, Waves) 'If False, calculates the non-local transport and tendencies but\n'//& 'purely for diagnostic purposes.', & default=.not. CS%passiveMode) - call get_param(paramFile, mdl, 'SMOOTH_BLD', CS%smoothBLD, & - 'If True, applies a 1-1-4-1-1 Laplacian filter one time on HBLT.\n'// & - 'computed via CVMix to reduce any horizontal two-grid-point noise.', & - default=.false.) + call get_param(paramFile, mdl, 'N_SMOOTH', CS%n_smooth, & + 'The number of times the 1-1-4-1-1 Laplacian filter is applied on\n'// & + 'OBL depth.', & + default=0) + if (CS%n_smooth > 0) then + call get_param(paramFile, mdl, 'DEEPEN_ONLY_VIA_SMOOTHING', CS%deepen_only, & + 'If true, apply OBLdepth smoothing at a cell only if the OBLdepth.\n'// & + 'gets deeper via smoothing.', & + default=.false.) + endif call get_param(paramFile, mdl, 'RI_CRIT', CS%Ri_crit, & 'Critical bulk Richardson number used to define depth of the\n'// & 'surface Ocean Boundary Layer (OBL).', & @@ -1498,7 +1505,7 @@ subroutine KPP_compute_BLD(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, buoyFlux, enddo enddo - if (CS%smoothBLD) call KPP_smooth_BLD(CS,G,GV,h) + if (CS%n_smooth > 0) call KPP_smooth_BLD(CS,G,GV,h) end subroutine KPP_compute_BLD @@ -1518,38 +1525,45 @@ subroutine KPP_smooth_BLD(CS,G,GV,h) real :: wc, ww, we, wn, ws ! averaging weights for smoothing real :: dh ! The local thickness used for calculating interface positions (m) real :: hcorr ! A cumulative correction arising from inflation of vanished layers (m) - integer :: i, j, k + integer :: i, j, k, s - ! Update halos - call pass_var(CS%OBLdepth, G%Domain) + do s=1,CS%n_smooth - OBLdepth_original = CS%OBLdepth - CS%OBLdepth_original = OBLdepth_original + ! Update halos + call pass_var(CS%OBLdepth, G%Domain) - ! apply smoothing on OBL depth - do j = G%jsc, G%jec - do i = G%isc, G%iec + OBLdepth_original = CS%OBLdepth + CS%OBLdepth_original = OBLdepth_original - ! skip land points - if (G%mask2dT(i,j)==0.) cycle + ! apply smoothing on OBL depth + do j = G%jsc, G%jec + do i = G%isc, G%iec - ! 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) - - CS%OBLdepth(i,j) = wc * OBLdepth_original(i,j) & - + ww * OBLdepth_original(i-1,j) & - + we * OBLdepth_original(i+1,j) & - + ws * OBLdepth_original(i,j-1) & - + wn * OBLdepth_original(i,j+1) + ! 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) + + CS%OBLdepth(i,j) = wc * OBLdepth_original(i,j) & + + ww * OBLdepth_original(i-1,j) & + + we * OBLdepth_original(i+1,j) & + + ws * OBLdepth_original(i,j-1) & + + wn * OBLdepth_original(i,j+1) + enddo enddo - enddo - ! prevent OBL depths deeper than the bathymetric depth - where (CS%OBLdepth > G%bathyT) CS%OBLdepth = G%bathyT + ! Apply OBLdepth smoothing at a cell only if the OBLdepth gets deeper via smoothing. + if (CS%deepen_only) CS%OBLdepth = max(CS%OBLdepth,CS%OBLdepth_original) + + ! prevent OBL depths deeper than the bathymetric depth + where (CS%OBLdepth > G%bathyT) CS%OBLdepth = G%bathyT + + enddo ! s-loop ! Update kOBL for smoothed OBL depths do j = G%jsc, G%jec @@ -1571,7 +1585,7 @@ subroutine KPP_smooth_BLD(CS,G,GV,h) iFaceHeight(k+1) = iFaceHeight(k) - dh enddo - CS%kOBL(i,j) = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, CS%OBLdepth(i,j) ) + CS%kOBL(i,j) = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, CS%OBLdepth(i,j) ) enddo enddo From 21d1038b65d14f9337d60f4412fc0fed27d667aa Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Thu, 31 May 2018 18:51:36 -0600 Subject: [PATCH 6/6] Check if Kv is associated before updating it --- .../vertical/MOM_tidal_mixing.F90 | 23 +++++++++---------- 1 file changed, 11 insertions(+), 12 deletions(-) diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90 index b659e9149a..226f7c4918 100644 --- a/src/parameterizations/vertical/MOM_tidal_mixing.F90 +++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90 @@ -754,9 +754,11 @@ subroutine calculate_CVMix_tidal(h, j, G, GV, CS, N2_int, Kd, Kv) enddo ! Update viscosity - do k=1,G%ke+1 - Kv(i,j,k) = Kv(i,j,k) + Kv_tidal(k) - enddo + if (associated(Kv)) then + do k=1,G%ke+1 + Kv(i,j,k) = Kv(i,j,k) + Kv_tidal(k) + enddo + endif ! diagnostics if (associated(dd%Kd_itidal)) then @@ -851,9 +853,11 @@ subroutine calculate_CVMix_tidal(h, j, G, GV, CS, N2_int, Kd, Kv) enddo ! Update viscosity - do k=1,G%ke+1 - Kv(i,j,k) = Kv(i,j,k) + Kv_tidal(k) - enddo + if (associated(Kv)) then + do k=1,G%ke+1 + Kv(i,j,k) = Kv(i,j,k) + Kv_tidal(k) + enddo + endif ! diagnostics if (associated(dd%Kd_itidal)) then @@ -903,12 +907,7 @@ subroutine add_int_tide_diffusivity(h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, CS, real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), optional, intent(inout) :: Kd_int real, intent(inout) :: Kd_max - ! This subroutine adds the effect of internal-tide-driven mixing to the layer diffusivities. - ! The mechanisms considered are (1) local dissipation of internal waves generated by the - ! barotropic flow ("itidal"), (2) local dissipation of internal waves generated by the propagating - ! low modes (rays) of the internal tide ("lowmode"), and (3) local dissipation of internal lee waves. - ! Will eventually need to add diffusivity due to other wave-breaking processes (e.g. Bottom friction, - ! Froude-number-depending breaking, PSI, etc.). + ! local real, dimension(SZI_(G)) :: & htot, & ! total thickness above or below a layer, or the