From 0d2af76aaad719a7c5c1d26df4ab06040e8eafe3 Mon Sep 17 00:00:00 2001 From: Raphael Dussin Date: Mon, 2 Aug 2021 22:10:42 -0400 Subject: [PATCH 1/6] fix refraction code --- .../lateral/MOM_internal_tides.F90 | 42 +++++++++++++++---- 1 file changed, 34 insertions(+), 8 deletions(-) diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index 2bb3c3b0f1..091e131ef3 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -769,6 +769,8 @@ subroutine refract(En, cn, freq, dt, G, US, NAngle, use_PPMang) Flux_E real, dimension(SZI_(G),SZJ_(G),1-stencil:NAngle+stencil) :: & CFL_ang + real, dimension(G%IsdB:G%IedB,G%jsd:G%jed) :: cn_u !< Internal wave group velocity at U-point + real, dimension(G%isd:G%ied,G%JsdB:G%JedB) :: cn_v !< Internal wave group velocity at V-point real :: f2 ! The squared Coriolis parameter [T-2 ~> s-2]. real :: favg ! The average Coriolis parameter at a point [T-1 ~> s-1]. real :: df_dy, df_dx ! The x- and y- gradients of the Coriolis parameter [T-1 L-1 ~> s-1 m-1]. @@ -779,9 +781,36 @@ subroutine refract(En, cn, freq, dt, G, US, NAngle, use_PPMang) real :: cn_subRO ! A tiny wave speed to prevent division by zero [L T-1 ~> m s-1] integer :: is, ie, js, je, asd, aed, na integer :: i, j, a + real :: wgt1, wgt2 + real :: eps is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; na = size(En,3) asd = 1-stencil ; aed = NAngle+stencil + eps=1.0e-20 * US%m_s_to_L_T + + do i=is-1,ie ; do j=js-1,je + wgt1 = 1. + wgt2 = 1. + if (cn(i,j) < eps) wgt1 = 0. + if (cn(i+1,j) < eps) wgt2 = 0. + if (wgt1 + wgt2 >= 1.) then + cn_u(I,j) = (cn(i,j) + cn(i+1,j)) / (wgt1 + wgt2) + else + cn_u(I,j) = 0. + endif + enddo ; enddo + + do i=is-1,ie ; do j=js-1,je + wgt1 = 1. + wgt2 = 1. + if (cn(i,j) < eps) wgt1 = 0. + if (cn(i,j+1) < eps) wgt2 = 0. + if (wgt1 + wgt2 >= 1.) then + cn_v(i,J) = (cn(i,j) + cn(i,j+1)) / (wgt1 + wgt2) + else + cn_v(i,J) = 0. + endif + enddo ; enddo Ifreq = 1.0 / freq cn_subRO = 1e-100*US%m_s_to_L_T ! The hard-coded value here might need to increase. @@ -813,16 +842,13 @@ subroutine refract(En, cn, freq, dt, G, US, NAngle, use_PPMang) (G%CoriolisBu(I,J-1) + G%CoriolisBu(I-1,J))) df_dx = 0.5*((G%CoriolisBu(I,J) + G%CoriolisBu(I,J-1)) - & (G%CoriolisBu(I-1,J) + G%CoriolisBu(I-1,J-1))) * G%IdxT(i,j) - dlnCn_dx = 0.5*( G%IdxCu(I,j) * (cn(i+1,j) - cn(i,j)) / & - (0.5*(cn(i+1,j) + cn(i,j)) + cn_subRO) + & - G%IdxCu(I-1,j) * (cn(i,j) - cn(i-1,j)) / & - (0.5*(cn(i,j) + cn(i-1,j)) + cn_subRO) ) + dlnCn_dx = G%IdxT(i,j) * (cn_u(I,j) - cn_u(I-1,j)) / (cn(i,j) + cn_subRO) + + df_dy = 0.5*((G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J)) - & (G%CoriolisBu(I,J-1) + G%CoriolisBu(I-1,J-1))) * G%IdyT(i,j) - dlnCn_dy = 0.5*( G%IdyCv(i,J) * (cn(i,j+1) - cn(i,j)) / & - (0.5*(cn(i,j+1) + cn(i,j)) + cn_subRO) + & - G%IdyCv(i,J-1) * (cn(i,j) - cn(i,j-1)) / & - (0.5*(cn(i,j) + cn(i,j-1)) + cn_subRO) ) + dlnCn_dy = G%IdyT(i,j) * (cn_v(i,J) - cn_v(i,J-1)) / (cn(i,j) + cn_subRO) + Kmag2 = (freq**2 - f2) / (cn(i,j)**2 + cn_subRO**2) if (Kmag2 > 0.0) then I_Kmag = 1.0 / sqrt(Kmag2) From 268f467a76c8a5680dbd67f051cd98896d37c006 Mon Sep 17 00:00:00 2001 From: Raphael Dussin Date: Mon, 9 Aug 2021 10:55:46 -0400 Subject: [PATCH 2/6] fix indexing and potential division by small number --- .../lateral/MOM_internal_tides.F90 | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index 091e131ef3..dab256f4ab 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -788,7 +788,7 @@ subroutine refract(En, cn, freq, dt, G, US, NAngle, use_PPMang) asd = 1-stencil ; aed = NAngle+stencil eps=1.0e-20 * US%m_s_to_L_T - do i=is-1,ie ; do j=js-1,je + do i=is-1,ie ; do j=js,je wgt1 = 1. wgt2 = 1. if (cn(i,j) < eps) wgt1 = 0. @@ -800,7 +800,7 @@ subroutine refract(En, cn, freq, dt, G, US, NAngle, use_PPMang) endif enddo ; enddo - do i=is-1,ie ; do j=js-1,je + do i=is,ie ; do j=js-1,je wgt1 = 1. wgt2 = 1. if (cn(i,j) < eps) wgt1 = 0. @@ -842,12 +842,16 @@ subroutine refract(En, cn, freq, dt, G, US, NAngle, use_PPMang) (G%CoriolisBu(I,J-1) + G%CoriolisBu(I-1,J))) df_dx = 0.5*((G%CoriolisBu(I,J) + G%CoriolisBu(I,J-1)) - & (G%CoriolisBu(I-1,J) + G%CoriolisBu(I-1,J-1))) * G%IdxT(i,j) - dlnCn_dx = G%IdxT(i,j) * (cn_u(I,j) - cn_u(I-1,j)) / (cn(i,j) + cn_subRO) - - df_dy = 0.5*((G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J)) - & (G%CoriolisBu(I,J-1) + G%CoriolisBu(I-1,J-1))) * G%IdyT(i,j) - dlnCn_dy = G%IdyT(i,j) * (cn_v(i,J) - cn_v(i,J-1)) / (cn(i,j) + cn_subRO) + + if (cn(i,j) < eps) then + dlnCn_dx = 0. + dlnCn_dy = 0. + else + dlnCn_dx = G%IdxT(i,j) * (cn_u(I,j) - cn_u(I-1,j)) / cn(i,j) + dlnCn_dy = G%IdyT(i,j) * (cn_v(i,J) - cn_v(i,J-1)) / cn(i,j) + endif Kmag2 = (freq**2 - f2) / (cn(i,j)**2 + cn_subRO**2) if (Kmag2 > 0.0) then From 9b97495e91adfe5fcd949e40b727892c74fb4832 Mon Sep 17 00:00:00 2001 From: Raphael Dussin Date: Mon, 9 Aug 2021 11:47:52 -0400 Subject: [PATCH 3/6] replace denominator of dlnCn_dxy with a more solid schme --- src/parameterizations/lateral/MOM_internal_tides.F90 | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index dab256f4ab..5932a43fd0 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -813,7 +813,7 @@ subroutine refract(En, cn, freq, dt, G, US, NAngle, use_PPMang) enddo ; enddo Ifreq = 1.0 / freq - cn_subRO = 1e-100*US%m_s_to_L_T ! The hard-coded value here might need to increase. + cn_subRO = 1e-30*US%m_s_to_L_T ! The hard-coded value here might need to increase. Angle_size = (8.0*atan(1.0)) / (real(NAngle)) dt_Angle_size = dt / Angle_size @@ -845,13 +845,8 @@ subroutine refract(En, cn, freq, dt, G, US, NAngle, use_PPMang) df_dy = 0.5*((G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J)) - & (G%CoriolisBu(I,J-1) + G%CoriolisBu(I-1,J-1))) * G%IdyT(i,j) - if (cn(i,j) < eps) then - dlnCn_dx = 0. - dlnCn_dy = 0. - else - dlnCn_dx = G%IdxT(i,j) * (cn_u(I,j) - cn_u(I-1,j)) / cn(i,j) - dlnCn_dy = G%IdyT(i,j) * (cn_v(i,J) - cn_v(i,J-1)) / cn(i,j) - endif + dlnCn_dx = G%IdxT(i,j) * (cn_u(I,j) - cn_u(I-1,j)) / (0.5 * (cn_u(I,j) + cn_u(I-1,j)) + cn_subRO) + dlnCn_dy = G%IdyT(i,j) * (cn_v(i,J) - cn_v(i,J-1)) / (0.5 * (cn_v(i,J) + cn_v(i,J-1)) + cn_subRO) Kmag2 = (freq**2 - f2) / (cn(i,j)**2 + cn_subRO**2) if (Kmag2 > 0.0) then From f6ab01daecaa34cea8752243992866be0621f84a Mon Sep 17 00:00:00 2001 From: Raphael Dussin Date: Mon, 16 Aug 2021 15:13:54 -0400 Subject: [PATCH 4/6] more elegant formulation for cn_u/cn_v --- .../lateral/MOM_internal_tides.F90 | 35 ++++++++----------- 1 file changed, 14 insertions(+), 21 deletions(-) diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index 389243a30c..37fe3ff8b9 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -779,6 +779,7 @@ subroutine refract(En, cn, freq, dt, G, US, NAngle, use_PPMang) CFL_ang real, dimension(G%IsdB:G%IedB,G%jsd:G%jed) :: cn_u !< Internal wave group velocity at U-point real, dimension(G%isd:G%ied,G%JsdB:G%JedB) :: cn_v !< Internal wave group velocity at V-point + real, dimension(G%isd:G%ied,G%jsd:G%jed) :: cnmask !< Local mask for group velocity real :: f2 ! The squared Coriolis parameter [T-2 ~> s-2]. real :: favg ! The average Coriolis parameter at a point [T-1 ~> s-1]. real :: df_dy, df_dx ! The x- and y- gradients of the Coriolis parameter [T-1 L-1 ~> s-1 m-1]. @@ -796,32 +797,24 @@ subroutine refract(En, cn, freq, dt, G, US, NAngle, use_PPMang) asd = 1-stencil ; aed = NAngle+stencil eps=1.0e-20 * US%m_s_to_L_T - do i=is-1,ie ; do j=js,je - wgt1 = 1. - wgt2 = 1. - if (cn(i,j) < eps) wgt1 = 0. - if (cn(i+1,j) < eps) wgt2 = 0. - if (wgt1 + wgt2 >= 1.) then - cn_u(I,j) = (cn(i,j) + cn(i+1,j)) / (wgt1 + wgt2) - else - cn_u(I,j) = 0. - endif + cnmask = merge(1.,0.,cn(:,:) > eps) + + do j=js,je ; do i=is-1,ie + ! wgt = 0 if local cn == 0, wgt = 0.5 if both contiguous values != 0 + ! and wgt = 1 if neighbour cn == 0 + wgt1 = cnmask(i,j)*(1-0.5*cnmask(i+1,j)) + wgt2 = cnmask(i+1,j)*(1-0.5*cnmask(i,j)) + cn_u(I,j) = wgt1 * cn(i,j) + wgt2 * cn(i+1,j) enddo ; enddo - do i=is,ie ; do j=js-1,je - wgt1 = 1. - wgt2 = 1. - if (cn(i,j) < eps) wgt1 = 0. - if (cn(i,j+1) < eps) wgt2 = 0. - if (wgt1 + wgt2 >= 1.) then - cn_v(i,J) = (cn(i,j) + cn(i,j+1)) / (wgt1 + wgt2) - else - cn_v(i,J) = 0. - endif + do j=js-1,je ; do i=is,ie + wgt1 = cnmask(i,j)*(1-0.5*cnmask(i,j+1)) + wgt2 = cnmask(i,j+1)*(1-0.5*cnmask(i,j)) + cn_v(i,J) = wgt1 * cn(i,j) + wgt2 * cn(i,j+1) enddo ; enddo Ifreq = 1.0 / freq - cn_subRO = 1e-30*US%m_s_to_L_T ! The hard-coded value here might need to increase. + cn_subRO = 1e-30*US%m_s_to_L_T Angle_size = (8.0*atan(1.0)) / (real(NAngle)) dt_Angle_size = dt / Angle_size From c703d2241b668b05d05012fc7873a8ee9c8bd307 Mon Sep 17 00:00:00 2001 From: Raphael Dussin Date: Mon, 16 Aug 2021 15:32:30 -0400 Subject: [PATCH 5/6] better definition of cnmask --- src/parameterizations/lateral/MOM_internal_tides.F90 | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index 37fe3ff8b9..dee463b8ff 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -791,13 +791,11 @@ subroutine refract(En, cn, freq, dt, G, US, NAngle, use_PPMang) integer :: is, ie, js, je, asd, aed, na integer :: i, j, a real :: wgt1, wgt2 - real :: eps is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; na = size(En,3) asd = 1-stencil ; aed = NAngle+stencil - eps=1.0e-20 * US%m_s_to_L_T - cnmask = merge(1.,0.,cn(:,:) > eps) + cnmask(:,:) = merge(0.,1.,cn(:,:) == 0.) do j=js,je ; do i=is-1,ie ! wgt = 0 if local cn == 0, wgt = 0.5 if both contiguous values != 0 From bcbf43c7a80d2cbe03282ada79d7a32903675c7a Mon Sep 17 00:00:00 2001 From: Raphael Dussin Date: Mon, 16 Aug 2021 15:39:14 -0400 Subject: [PATCH 6/6] cosmetics --- .../lateral/MOM_internal_tides.F90 | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index dee463b8ff..4d66471408 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -795,20 +795,20 @@ subroutine refract(En, cn, freq, dt, G, US, NAngle, use_PPMang) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; na = size(En,3) asd = 1-stencil ; aed = NAngle+stencil - cnmask(:,:) = merge(0.,1.,cn(:,:) == 0.) + cnmask(:,:) = merge(0., 1., cn(:,:) == 0.) do j=js,je ; do i=is-1,ie ! wgt = 0 if local cn == 0, wgt = 0.5 if both contiguous values != 0 ! and wgt = 1 if neighbour cn == 0 - wgt1 = cnmask(i,j)*(1-0.5*cnmask(i+1,j)) - wgt2 = cnmask(i+1,j)*(1-0.5*cnmask(i,j)) - cn_u(I,j) = wgt1 * cn(i,j) + wgt2 * cn(i+1,j) + wgt1 = cnmask(i,j) - 0.5 * cnmask(i,j) * cnmask(i+1,j) + wgt2 = cnmask(i+1,j) - 0.5 * cnmask(i,j) * cnmask(i+1,j) + cn_u(I,j) = wgt1*cn(i,j) + wgt2*cn(i+1,j) enddo ; enddo do j=js-1,je ; do i=is,ie - wgt1 = cnmask(i,j)*(1-0.5*cnmask(i,j+1)) - wgt2 = cnmask(i,j+1)*(1-0.5*cnmask(i,j)) - cn_v(i,J) = wgt1 * cn(i,j) + wgt2 * cn(i,j+1) + wgt1 = cnmask(i,j) - 0.5 * cnmask(i,j) * cnmask(i,j+1) + wgt2 = cnmask(i,j+1) - 0.5 * cnmask(i,j) * cnmask(i,j+1) + cn_v(i,J) = wgt1*cn(i,j) + wgt2*cn(i,j+1) enddo ; enddo Ifreq = 1.0 / freq