From ae54d2fccb8f16bc70568bf1c331b4be18e84f66 Mon Sep 17 00:00:00 2001 From: William Cooke Date: Mon, 22 Apr 2019 13:56:20 -0400 Subject: [PATCH 1/7] Revert "Implement changes suggested by @Hallberg-NOAA" This reverts commit 8f4af3d9ef927dc4b99d2a44f32a1e0a3ca5c2c3. --- src/tracer/MOM_neutral_diffusion.F90 | 13 ++++++------- src/tracer/MOM_neutral_diffusion_aux.F90 | 2 +- 2 files changed, 7 insertions(+), 8 deletions(-) diff --git a/src/tracer/MOM_neutral_diffusion.F90 b/src/tracer/MOM_neutral_diffusion.F90 index d5a6f45c5f..bca6dd61ca 100644 --- a/src/tracer/MOM_neutral_diffusion.F90 +++ b/src/tracer/MOM_neutral_diffusion.F90 @@ -246,9 +246,6 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, h, T, S, CS) integer :: iMethod real, dimension(SZI_(G)) :: ref_pres ! Reference pressure used to calculate alpha/beta real :: h_neglect, h_neglect_edge - real :: pa_to_H - - pa_to_H = 1. / GV%H_to_pa !### Try replacing both of these with GV%H_subroundoff if (GV%Boussinesq) then @@ -304,13 +301,15 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, h, T, S, CS) else ! Discontinuous reconstruction do k = 1, G%ke if (CS%ref_pres<0) ref_pres(:) = CS%Pint(:,j,k) - call calculate_density_derivs(CS%T_i(:,j,k,1), CS%S_i(:,j,k,1), ref_pres, & - CS%dRdT_i(:,j,k,1), CS%dRdS_i(:,j,k,1), G%isc-1, G%iec-G%isc+3, CS%EOS) + if (CS%stable_cell(i,j,k)) & + call calculate_density_derivs(CS%T_i(:,j,k,1), CS%S_i(:,j,k,1), ref_pres, & + CS%dRdT_i(:,j,k,1), CS%dRdS_i(:,j,k,1), G%isc-1, G%iec-G%isc+3, CS%EOS) if (CS%ref_pres<0) then ref_pres(:) = CS%Pint(:,j,k+1) endif - call calculate_density_derivs(CS%T_i(:,j,k,2), CS%S_i(:,j,k,2), ref_pres, & - CS%dRdT_i(:,j,k,2), CS%dRdS_i(:,j,k,2), G%isc-1, G%iec-G%isc+3, CS%EOS) + if (CS%stable_cell(i,j,k)) & + call calculate_density_derivs(CS%T_i(:,j,k,2), CS%S_i(:,j,k,2), ref_pres, & + CS%dRdT_i(:,j,k,2), CS%dRdS_i(:,j,k,2), G%isc-1, G%iec-G%isc+3, CS%EOS) enddo endif enddo diff --git a/src/tracer/MOM_neutral_diffusion_aux.F90 b/src/tracer/MOM_neutral_diffusion_aux.F90 index 88df1ddbc5..5152a02b88 100644 --- a/src/tracer/MOM_neutral_diffusion_aux.F90 +++ b/src/tracer/MOM_neutral_diffusion_aux.F90 @@ -161,7 +161,7 @@ real function calc_drho(T1, S1, dRdT1, dRdS1, T2, S2, dRdT2, dRdS2) end function calc_drho !> Calculate the difference in neutral density between a reference T, S, alpha, and beta -!! at a point on the polynomial reconstructions of T, S +!! and a poiet on the polynomial reconstructions of T, S subroutine drho_at_pos(CS, T_ref, S_ref, alpha_ref, beta_ref, P_top, P_bot, ppoly_T, ppoly_S, x0, & delta_rho, P_out, T_out, S_out, alpha_avg_out, beta_avg_out, delta_T_out, delta_S_out) type(ndiff_aux_CS_type), intent(in) :: CS !< Control structure with parameters for this module From 1e68a2fb1068a523429f9828c51b8462bc13e8d4 Mon Sep 17 00:00:00 2001 From: William Cooke Date: Mon, 22 Apr 2019 14:04:51 -0400 Subject: [PATCH 2/7] Revert "*Corrected the clock as seen by diabatic processes" This reverts commit bc6c6e65d658f7cdddf0c589ae770feb40287c01. --- src/core/MOM.F90 | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 1b364d7ef0..e345eea4c9 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -40,7 +40,7 @@ module MOM use MOM_restart, only : register_restart_field, query_initialized, save_restart use MOM_restart, only : restart_init, is_new_run, MOM_restart_CS use MOM_spatial_means, only : global_mass_integral -use MOM_time_manager, only : time_type, real_to_time, time_type_to_real, operator(+) +use MOM_time_manager, only : time_type, real_to_time, set_time, time_type_to_real, operator(+) use MOM_time_manager, only : operator(-), operator(>), operator(*), operator(/) use MOM_time_manager, only : operator(>=), increment_date use MOM_unit_tests, only : unit_tests @@ -615,12 +615,15 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS, & rel_time = 0.0 do n=1,n_max - rel_time = rel_time + dt ! The relative time at the end of the step. - ! Set the universally visible time to the middle of the time step. - CS%Time = Time_start + real_to_time(rel_time - 0.5*dt) - ! Set the local time to the end of the time step. - Time_local = Time_start + real_to_time(rel_time) + nt_debug = nt_debug + 1 + + ! Set the universally visible time to the middle of the time step + CS%Time = Time_start + set_time(int(floor(CS%rel_time+0.5*dt+0.5))) + CS%rel_time = CS%rel_time + dt + + ! Set the local time to the end of the time step. + Time_local = Time_start + set_time(int(floor(CS%rel_time+0.5))) if (showCallTree) call callTree_enter("DT cycles (step_MOM) n=",n) !=========================================================================== From deab21f96fbfd19b24bd29c60657c770a0de3217 Mon Sep 17 00:00:00 2001 From: William Cooke Date: Mon, 22 Apr 2019 14:06:56 -0400 Subject: [PATCH 3/7] Revert "*Use rho_ref in finite volume PGF density calcs" This reverts commit 48e90d0b928457e64c80b065c32f1131cc6b6dfb. --- src/equation_of_state/MOM_EOS.F90 | 107 ++++++++++++++++++------------ 1 file changed, 66 insertions(+), 41 deletions(-) diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index b06ffa0a79..9428307bd8 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -968,10 +968,11 @@ subroutine int_density_dz_generic(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, T5(n) = T(i,j) ; S5(n) = S(i,j) p5(n) = -GxRho*(z_t(i,j) - 0.25*real(n-1)*dz) enddo - call calculate_density(T5, S5, p5, r5, 1, 5, EOS, rho_ref) + call calculate_density(T5, S5, p5, r5, 1, 5, EOS) !, rho_ref) ! Use Bode's rule to estimate the pressure anomaly change. - rho_anom = C1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3)) + rho_anom = C1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3)) - & + rho_ref dpa(i-ioff,j-joff) = G_e*dz*rho_anom ! Use a Bode's-rule-like fifth-order accurate estimate of the double integral of ! the pressure anomaly. @@ -1010,10 +1011,11 @@ subroutine int_density_dz_generic(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, do n=2,5 T5(n) = T5(1) ; S5(n) = S5(1) ; p5(n) = p5(n-1) + GxRho*0.25*dz enddo - call calculate_density(T5, S5, p5, r5, 1, 5, EOS, rho_ref) + call calculate_density(T5, S5, p5, r5, 1, 5, EOS) !, rho_ref) ! Use Bode's rule to estimate the pressure anomaly change. - intz(m) = G_e*dz*( C1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3))) + intz(m) = G_e*dz*( C1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + & + 12.0*r5(3)) - rho_ref) enddo ! Use Bode's rule to integrate the bottom pressure anomaly values in x. intx_dpa(i-ioff,j-joff) = C1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + & @@ -1052,10 +1054,11 @@ subroutine int_density_dz_generic(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, T5(n) = T5(1) ; S5(n) = S5(1) p5(n) = p5(n-1) + GxRho*0.25*dz enddo - call calculate_density(T5, S5, p5, r5, 1, 5, EOS, rho_ref) + call calculate_density(T5, S5, p5, r5, 1, 5, EOS) !, rho_ref) ! Use Bode's rule to estimate the pressure anomaly change. - intz(m) = G_e*dz*( C1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3))) + intz(m) = G_e*dz*( C1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + & + 12.0*r5(3)) - rho_ref) enddo ! Use Bode's rule to integrate the values. inty_dpa(i-ioff,j-joff) = C1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + & @@ -1121,34 +1124,54 @@ subroutine int_density_dz_generic_plm (T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, & ! It is assumed that the salinity and temperature profiles are linear in the ! vertical. The top and bottom values within each layer are provided and ! a linear interpolation is used to compute intermediate values. - - ! Local variables - real :: T5((5*HIO%iscB+1):(5*(HIO%iecB+2))) ! Temperatures along a line of subgrid locations [degC]. - real :: S5((5*HIO%iscB+1):(5*(HIO%iecB+2))) ! Salinities along a line of subgrid locations [ppt]. - real :: p5((5*HIO%iscB+1):(5*(HIO%iecB+2))) ! Pressures along a line of subgrid locations [Pa]. - real :: r5((5*HIO%iscB+1):(5*(HIO%iecB+2))) ! Densities along a line of subgrid locations [kg m-3]. - real :: T15((15*HIO%iscB+1):(15*(HIO%iecB+1))) ! Temperatures at an array of subgrid locations [degC]. - real :: S15((15*HIO%iscB+1):(15*(HIO%iecB+1))) ! Salinities at an array of subgrid locations [ppt]. - real :: p15((15*HIO%iscB+1):(15*(HIO%iecB+1))) ! Pressures at an array of subgrid locations [Pa]. - real :: r15((15*HIO%iscB+1):(15*(HIO%iecB+1))) ! Densities at an array of subgrid locations [kg m-3]. - real :: wt_t(5), wt_b(5) ! Top and bottom weights [nondim]. - real :: rho_anom ! A density anomaly [kg m-3]. - real :: w_left, w_right ! Left and right weights [nondim]. - real :: intz(5) ! The gravitational acceleration times the integrals of density - ! with height at the 5 sub-column locations [Pa]. - real, parameter :: C1_90 = 1.0/90.0 ! A rational constant [nondim]. - real :: GxRho ! Gravitational acceleration times density [kg m-1 Z-1 s-2 ~> kg m-2 s-2]. - real :: I_Rho ! The inverse of the reference density [m3 kg-1]. - real :: dz(HIO%iscB:HIO%iecB+1) ! Layer thicknesses at tracer points [Z ~> m]. - real :: dz_x(5,HIO%iscB:HIO%iecB) ! Layer thicknesses along an x-line of subrid locations [Z ~> m]. - real :: dz_y(5,HIO%isc:HIO%iec) ! Layer thicknesses along a y-line of subrid locations [Z ~> m]. - real :: weight_t, weight_b ! Nondimensional wieghts of the top and bottom. - real :: massWeightToggle ! A nondimensional toggle factor (0 or 1). - real :: Ttl, Tbl, Ttr, Tbr ! Temperatures at the velocity cell corners [degC]. - real :: Stl, Sbl, Str, Sbr ! Salinities at the velocity cell corners [ppt]. - real :: hWght ! A topographically limited thicknes weight [Z ~> m]. - real :: hL, hR ! Thicknesses to the left and right [Z ~> m]. - real :: iDenom ! The denominator of the thickness weight expressions [Z-2 ~> m-2]. +! +! Arguments: T - potential temperature relative to the surface in C +! (the 't' and 'b' subscripts refer to the values at +! the top and the bottom of each layer) +! (in) S - salinity in PSU. +! (the 't' and 'b' subscripts refer to the values at +! the top and the bottom of each layer) +! (in) z_t - height at the top of the layer in m. +! (in) z_b - height at the top of the layer in m. +! (in) rho_ref - A mean density, in kg m-3, that is subtracted out to reduce +! the magnitude of each of the integrals. +! (The pressure is calucated as p~=-z*rho_0*G_e.) +! (in) rho_0 - A density, in kg m-3, that is used to calculate the pressure +! (as p~=-z*rho_0*G_e) used in the equation of state. +! (in) G_e - The Earth's gravitational acceleration, in m s-2. +! (in) G - The ocean's grid structure. +! (in) form_of_eos - integer that selects the eqn of state. +! (out) dpa - The change in the pressure anomaly across the layer, +! in Pa. +! (out,opt) intz_dpa - The integral through the thickness of the layer of the +! pressure anomaly relative to the anomaly at the top of +! the layer, in Pa m. +! (out,opt) intx_dpa - The integral in x of the difference between the +! pressure anomaly at the top and bottom of the layer +! divided by the x grid spacing, in Pa. +! (out,opt) inty_dpa - The integral in y of the difference between the +! pressure anomaly at the top and bottom of the layer +! divided by the y grid spacing, in Pa. +! (in,opt) useMassWghtInterp - If true, uses mass weighting to interpolate +! T/S for top and bottom integrals. + + real :: T5((5*HIO%iscB+1):(5*(HIO%iecB+2))) + real :: S5((5*HIO%iscB+1):(5*(HIO%iecB+2))) + real :: p5((5*HIO%iscB+1):(5*(HIO%iecB+2))) + real :: r5((5*HIO%iscB+1):(5*(HIO%iecB+2))) + real :: u5((5*HIO%iscB+1):(5*(HIO%iecB+2))) + real :: T15((15*HIO%iscB+1):(15*(HIO%iecB+1))) + real :: S15((15*HIO%iscB+1):(15*(HIO%iecB+1))) + real :: p15((15*HIO%iscB+1):(15*(HIO%iecB+1))) + real :: r15((15*HIO%iscB+1):(15*(HIO%iecB+1))) + real :: wt_t(5), wt_b(5) + real :: rho_anom + real :: w_left, w_right, intz(5) + real, parameter :: C1_90 = 1.0/90.0 ! Rational constants. + real :: GxRho, I_Rho + real :: dz(HIO%iscB:HIO%iecB+1), dz_x(5,HIO%iscB:HIO%iecB), dz_y(5,HIO%isc:HIO%iec) + real :: weight_t, weight_b, hWght, massWeightToggle + real :: Ttl, Tbl, Ttr, Tbr, Stl, Sbl, Str, Sbr, hL, hR, iDenom integer :: Isq, Ieq, Jsq, Jeq, i, j, m, n integer :: iin, jin, ioff, joff integer :: pos @@ -1184,17 +1207,19 @@ subroutine int_density_dz_generic_plm (T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, & T5(i*5+n) = wt_t(n) * T_t(iin,jin) + wt_b(n) * T_b(iin,jin) enddo enddo - call calculate_density_array(T5, S5, p5, r5, 1, (ieq-isq+2)*5, EOS, rho_ref ) + call calculate_density_array(T5, S5, p5, r5, 1, (ieq-isq+2)*5, EOS) !, rho_ref ) + u5 = r5 - rho_ref do i=isq,ieq+1 ; iin = i+ioff ! Use Bode's rule to estimate the pressure anomaly change. - rho_anom = C1_90*(7.0*(r5(i*5+1)+r5(i*5+5)) + 32.0*(r5(i*5+2)+r5(i*5+4)) + 12.0*r5(i*5+3)) + rho_anom = C1_90*(7.0*(r5(i*5+1)+r5(i*5+5)) + 32.0*(r5(i*5+2)+r5(i*5+4)) + 12.0*r5(i*5+3)) - & + rho_ref dpa(i,j) = G_e*dz(i)*rho_anom if (present(intz_dpa)) then ! Use a Bode's-rule-like fifth-order accurate estimate of ! the double integral of the pressure anomaly. intz_dpa(i,j) = 0.5*G_e*dz(i)**2 * & - (rho_anom - C1_90*(16.0*(r5(i*5+4)-r5(i*5+2)) + 7.0*(r5(i*5+5)-r5(i*5+1))) ) + (rho_anom - C1_90*(16.0*(u5(i*5+4)-u5(i*5+2)) + 7.0*(u5(i*5+5)-u5(i*5+1))) ) endif enddo enddo ! end loops on j @@ -1264,7 +1289,7 @@ subroutine int_density_dz_generic_plm (T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, & enddo enddo - call calculate_density(T15, S15, p15, r15, 1, 15*(ieq-isq+1), EOS, rho_ref) + call calculate_density(T15, S15, p15, r15, 1, 15*(ieq-isq+1), EOS) !, rho_ref) do I=Isq,Ieq ; iin = i+ioff intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j) @@ -1273,7 +1298,7 @@ subroutine int_density_dz_generic_plm (T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, & do m = 2,4 pos = i*15+(m-2)*5 intz(m) = G_e*dz_x(m,i)*( C1_90*(7.0*(r15(pos+1)+r15(pos+5)) + 32.0*(r15(pos+2)+r15(pos+4)) + & - 12.0*r15(pos+3))) + 12.0*r15(pos+3)) - rho_ref) enddo ! Use Bode's rule to integrate the bottom pressure anomaly values in x. intx_dpa(i,j) = C1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + & @@ -1344,7 +1369,7 @@ subroutine int_density_dz_generic_plm (T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, & enddo call calculate_density_array(T15(15*HIO%isc+1:), S15(15*HIO%isc+1:), p15(15*HIO%isc+1:), & - r15(15*HIO%isc+1:), 1, 15*(HIO%iec-HIO%isc+1), EOS, rho_ref) + r15(15*HIO%isc+1:), 1, 15*(HIO%iec-HIO%isc+1), EOS) !, rho_ref) do i=HIO%isc,HIO%iec ; iin = i+ioff intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1) @@ -1353,7 +1378,7 @@ subroutine int_density_dz_generic_plm (T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, & pos = i*15+(m-2)*5 intz(m) = G_e*dz_y(m,i)*( C1_90*(7.0*(r15(pos+1)+r15(pos+5)) + & 32.0*(r15(pos+2)+r15(pos+4)) + & - 12.0*r15(pos+3))) + 12.0*r15(pos+3)) - rho_ref) enddo ! Use Bode's rule to integrate the values. inty_dpa(i,j) = C1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + & From 148d88e5b8fd37b2d1a66d8b828bb288d7148085 Mon Sep 17 00:00:00 2001 From: William Cooke Date: Mon, 22 Apr 2019 16:16:13 -0400 Subject: [PATCH 4/7] Update of MOM6 code to allow SPEAR to reproduce previous answers. Updates MOM code to dev/gfdl as of 2/7/2019 (6dd6f5295b9af1) and reverts 3 answer changing modifications. Update produced by git revert 8f4af3d9ef927dc4b99d2a44f32a1e0a3ca5c2c3 git add src/tracer/MOM_neutral_diffusion.F90 git revert --continue git revert bc6c6e65d658f git add src/core/MOM.F90 git revert --continue git revert 48e90d0b928457 git add src/equation_of_state/MOM_EOS.F90 git revert --continue Some conflict resolution was needed. --- src/core/MOM.F90 | 28 ++++++---------------------- src/tracer/MOM_neutral_diffusion.F90 | 8 ++------ 2 files changed, 8 insertions(+), 28 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index e345eea4c9..aa684f1f26 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -616,14 +616,12 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS, & rel_time = 0.0 do n=1,n_max - nt_debug = nt_debug + 1 - ! Set the universally visible time to the middle of the time step - CS%Time = Time_start + set_time(int(floor(CS%rel_time+0.5*dt+0.5))) - CS%rel_time = CS%rel_time + dt + CS%Time = Time_start + set_time(int(floor(rel_time+0.5*dt+0.5))) + rel_time = rel_time + dt ! Set the local time to the end of the time step. - Time_local = Time_start + set_time(int(floor(CS%rel_time+0.5))) + Time_local = Time_start + set_time(int(floor(rel_time+0.5))) if (showCallTree) call callTree_enter("DT cycles (step_MOM) n=",n) !=========================================================================== @@ -645,15 +643,9 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS, & dtdia = dt*min(ntstep,n_max-(n-1)) endif - end_time_thermo = Time_local - if (dtdia > dt) then - ! If necessary, temporarily reset CS%Time to the center of the period covered - ! by the call to step_MOM_thermo, noting that they begin at the same time. - CS%Time = CS%Time + real_to_time(0.5*(dtdia-dt)) - ! The end-time of the diagnostic interval needs to be set ahead if there - ! are multiple dynamic time steps worth of thermodynamics applied here. - end_time_thermo = Time_local + real_to_time(dtdia-dt) - endif + ! The end-time of the diagnostic interval needs to be set ahead if there + ! are multiple dynamic time steps worth of thermodynamics applied here. + end_time_thermo = Time_local + set_time(int(floor(dtdia-dt+0.5))) ! Apply diabatic forcing, do mixing, and regrid. call step_MOM_thermo(CS, G, GV, US, u, v, h, CS%tv, fluxes, dtdia, & @@ -664,8 +656,6 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS, & CS%t_dyn_rel_thermo = -dtdia if (showCallTree) call callTree_waypoint("finished diabatic_first (step_MOM)") - if (dtdia > dt) & ! Reset CS%Time to its previous value. - CS%Time = Time_start + real_to_time(rel_time - 0.5*dt) endif ! end of block "(CS%diabatic_first .and. (CS%t_dyn_rel_adv==0.0))" if (do_dyn) then @@ -745,18 +735,12 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS, & "before call to diabatic.") endif - ! If necessary, temporarily reset CS%Time to the center of the period covered - ! by the call to step_MOM_thermo, noting that they end at the same time. - if (dtdia > dt) CS%Time = CS%Time - real_to_time(0.5*(dtdia-dt)) - ! Apply diabatic forcing, do mixing, and regrid. call step_MOM_thermo(CS, G, GV, US, u, v, h, CS%tv, fluxes, dtdia, & Time_local, .false., Waves=Waves) CS%time_in_thermo_cycle = CS%time_in_thermo_cycle + dtdia CS%t_dyn_rel_thermo = 0.0 - if (dtdia > dt) & ! Reset CS%Time to its previous value. - CS%Time = Time_start + real_to_time(rel_time - 0.5*dt) endif if (do_dyn) then diff --git a/src/tracer/MOM_neutral_diffusion.F90 b/src/tracer/MOM_neutral_diffusion.F90 index bca6dd61ca..e692157778 100644 --- a/src/tracer/MOM_neutral_diffusion.F90 +++ b/src/tracer/MOM_neutral_diffusion.F90 @@ -380,12 +380,8 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, h, T, S, CS) ! calculates hEff from the fraction of the nondimensional fraction of the layer occupied by ! the... (Please finish this thought. -RWH) if (CS%continuous_reconstruction) then - do k = 1, CS%nsurf-1 ; do j = G%jsc, G%jec ; do I = G%isc-1, G%iec - if (G%mask2dCu(I,j) > 0.) CS%uhEff(I,j,k) = CS%uhEff(I,j,k) * pa_to_H - enddo ; enddo ; enddo - do k = 1, CS%nsurf-1 ; do J = G%jsc-1, G%jec ; do i = G%isc, G%iec - if (G%mask2dCv(i,J) > 0.) CS%vhEff(i,J,k) = CS%vhEff(i,J,k) * pa_to_H - enddo ; enddo ; enddo + CS%uhEff(:,:,:) = CS%uhEff(:,:,:) / GV%H_to_pa + CS%vhEff(:,:,:) = CS%vhEff(:,:,:) / GV%H_to_pa endif if (CS%id_uhEff_2d>0) then From 380803cda68e27d3cb9d301de302a0732253f5ed Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Mon, 1 Jul 2019 13:36:46 -0400 Subject: [PATCH 5/7] Optional use of differing restoring piston velocities for temp and salt --- .../coupled_driver/MOM_surface_forcing.F90 | 26 ++++++++++++++----- 1 file changed, 19 insertions(+), 7 deletions(-) diff --git a/config_src/coupled_driver/MOM_surface_forcing.F90 b/config_src/coupled_driver/MOM_surface_forcing.F90 index 5112a0b64b..0bd9f22afa 100644 --- a/config_src/coupled_driver/MOM_surface_forcing.F90 +++ b/config_src/coupled_driver/MOM_surface_forcing.F90 @@ -113,6 +113,8 @@ module MOM_surface_forcing logical :: restore_temp !< If true, the coupled MOM driver adds a term to restore sea !! surface temperature to a specified value. real :: Flux_const !< Piston velocity for surface restoring [m s-1] + real :: Flux_const_salt !< Piston velocity for surface salt restoring [m s-1] + real :: Flux_const_temp !< Piston velocity for surface temp restoring [m s-1] logical :: salt_restore_as_sflux !< If true, SSS restore as salt flux instead of water flux logical :: adjust_net_srestore_to_zero !< Adjust srestore to zero (for both salt_flux or vprec) logical :: adjust_net_srestore_by_scaling !< Adjust srestore w/o moving zero contour @@ -347,7 +349,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, US, CS, sfc do j=js,je ; do i=is,ie delta_sss = data_restore(i,j)- sfc_state%SSS(i,j) delta_sss = sign(1.0,delta_sss)*min(abs(delta_sss),CS%max_delta_srestore) - fluxes%salt_flux(i,j) = 1.e-3*G%mask2dT(i,j) * (CS%Rho0*CS%Flux_const)* & + fluxes%salt_flux(i,j) = 1.e-3*G%mask2dT(i,j) * (CS%Rho0*CS%Flux_const_salt)* & (CS%basin_mask(i,j)*open_ocn_mask(i,j)*CS%srestore_mask(i,j)) *delta_sss ! kg Salt m-2 s-1 enddo ; enddo if (CS%adjust_net_srestore_to_zero) then @@ -367,7 +369,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, US, CS, sfc delta_sss = sfc_state%SSS(i,j) - data_restore(i,j) delta_sss = sign(1.0,delta_sss)*min(abs(delta_sss),CS%max_delta_srestore) fluxes%vprec(i,j) = (CS%basin_mask(i,j)*open_ocn_mask(i,j)*CS%srestore_mask(i,j))* & - (CS%Rho0*CS%Flux_const) * & + (CS%Rho0*CS%Flux_const_salt) * & delta_sss / (0.5*(sfc_state%SSS(i,j) + data_restore(i,j))) endif enddo ; enddo @@ -393,7 +395,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, US, CS, sfc delta_sst = data_restore(i,j)- sfc_state%SST(i,j) delta_sst = sign(1.0,delta_sst)*min(abs(delta_sst),CS%max_delta_trestore) fluxes%heat_added(i,j) = G%mask2dT(i,j) * CS%trestore_mask(i,j) * & - (CS%Rho0*fluxes%C_p) * delta_sst * CS%Flux_const ! W m-2 + (CS%Rho0*fluxes%C_p) * delta_sst * CS%Flux_const_temp ! W m-2 enddo ; enddo endif @@ -1261,10 +1263,15 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS) if (CS%restore_salt) then call get_param(param_file, mdl, "FLUXCONST", CS%Flux_const, & - "The constant that relates the restoring surface fluxes "//& + "The constant that relates the restoring surface salt fluxes "//& "to the relative surface anomalies (akin to a piston "//& "velocity). Note the non-MKS units.", units="m day-1", & fail_if_missing=.true.) + call get_param(param_file, mdl, "FLUXCONST_SALT", CS%Flux_const_salt, & + "The constant that relates the restoring surface fluxes "//& + "to the relative surface anomalies (akin to a piston "//& + "velocity). Note the non-MKS units.", units="m day-1", & + fail_if_missing=.false., default=CS%Flux_const) call get_param(param_file, mdl, "SALT_RESTORE_FILE", CS%salt_restore_file, & "A file in which to find the surface salinity to use for restoring.", & default="salt_restore.nc") @@ -1273,7 +1280,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS) "SALT_RESTORE_FILE for restoring salinity.", & default="salt") ! Convert CS%Flux_const from m day-1 to m s-1. - CS%Flux_const = CS%Flux_const / 86400.0 + CS%Flux_const_salt = CS%Flux_const_salt / 86400.0 call get_param(param_file, mdl, "SRESTORE_AS_SFLUX", CS%salt_restore_as_sflux, & "If true, the restoring of salinity is applied as a salt "//& @@ -1309,10 +1316,15 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS) if (CS%restore_temp) then call get_param(param_file, mdl, "FLUXCONST", CS%Flux_const, & - "The constant that relates the restoring surface fluxes "//& + "The constant that relates the restoring surface temperature fluxes "//& "to the relative surface anomalies (akin to a piston "//& "velocity). Note the non-MKS units.", units="m day-1", & fail_if_missing=.true.) + call get_param(param_file, mdl, "FLUXCONST_TEMP", CS%Flux_const_temp, & + "The constant that relates the restoring surface temperature fluxes "//& + "to the relative surface anomalies (akin to a piston "//& + "velocity). Note the non-MKS units.", units="m day-1", & + fail_if_missing=.false., default=CS%Flux_const) call get_param(param_file, mdl, "SST_RESTORE_FILE", CS%temp_restore_file, & "A file in which to find the surface temperature to use for restoring.", & default="temp_restore.nc") @@ -1321,7 +1333,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS) "SST_RESTORE_FILE for restoring sst.", & default="temp") ! Convert CS%Flux_const from m day-1 to m s-1. - CS%Flux_const = CS%Flux_const / 86400.0 + CS%Flux_const_temp = CS%Flux_const_temp / 86400.0 call get_param(param_file, mdl, "MAX_DELTA_TRESTORE", CS%max_delta_trestore, & "The maximum sst difference used in restoring terms.", & From dac0c6c78be30ac9263721afab2dfb6311affd67 Mon Sep 17 00:00:00 2001 From: William Cooke Date: Thu, 5 Sep 2019 16:42:15 -0400 Subject: [PATCH 6/7] Correction to ePBL code to mitigate blowup. Also "corrects" some spelling differences in variables (MStar vs mstar) --- .../vertical/MOM_energetic_PBL.F90 | 32 +++++++++++-------- 1 file changed, 18 insertions(+), 14 deletions(-) diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index b486e1e2ca..5527866793 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -1772,16 +1772,16 @@ subroutine find_mstar(CS, US, Buoyancy_Flux, UStar, UStar_Mean,& if (CS%answers_2018) then ! The limit for the balance of rotation and stabilizing is f(L_Ekman,L_Obukhov) - MStar_S = CS%MStar_coef*sqrt(max(0.0,Buoyancy_Flux) / Ustar**2 / & + MStar_S = CS%MStar_coef*sqrt(max(0.0,Buoyancy_Flux) / UStar**2 / & (Abs_Coriolis + 1.e-10*US%T_to_s) ) ! The limit for rotation (Ekman length) limited mixing - MStar_N = CS%C_Ek * log( max( 1., Ustar / (Abs_Coriolis + 1.e-10*US%T_to_s) / BLD ) ) + MStar_N = CS%C_Ek * log( max( 1., UStar / (Abs_Coriolis + 1.e-10*US%T_to_s) / BLD ) ) else ! The limit for the balance of rotation and stabilizing is f(L_Ekman,L_Obukhov) - mstar_S = CS%MSTAR_COEF*sqrt(max(0.0, Buoyancy_Flux) / (Ustar**2 * max(Abs_Coriolis, 1.e-20*US%T_to_s))) + MStar_S = CS%MSTAR_COEF*sqrt(max(0.0, Buoyancy_Flux) / (UStar**2 * max(Abs_Coriolis, 1.e-20*US%T_to_s))) ! The limit for rotation (Ekman length) limited mixing - mstar_N = 0.0 - if (Ustar > Abs_Coriolis * BLD) mstar_N = CS%C_EK * log(Ustar / (Abs_Coriolis * BLD)) + MStar_N = 0.0 + if (UStar > Abs_Coriolis * BLD) Mstar_N = CS%C_EK * log(UStar / (Abs_Coriolis * BLD)) endif ! Here 1.25 is about .5/von Karman, which gives the Obukhov limit. @@ -1792,11 +1792,11 @@ subroutine find_mstar(CS, US, Buoyancy_Flux, UStar, UStar_Mean,& MStar_N = CS%RH18_MStar_cn1 * ( 1.0 - 1.0 / ( 1. + CS%RH18_MStar_cn2 * & exp( CS%RH18_mstar_CN3 * BLD * Abs_Coriolis / UStar) ) ) else - MSN_term = CS%RH18_MStar_cn2 * exp( CS%RH18_mstar_CN3 * BLD * Abs_Coriolis / Ustar) + MSN_term = CS%RH18_MStar_cn2 * exp( CS%RH18_mstar_CN3 * BLD * Abs_Coriolis / UStar) MStar_N = (CS%RH18_MStar_cn1 * MSN_term) / ( 1. + MSN_term) endif MStar_S = CS%RH18_MStar_CS1 * ( max(0.0, Buoyancy_Flux)**2 * BLD / & - ( Ustar**5 * max(Abs_Coriolis,1.e-20*US%T_to_s) ) )**CS%RH18_mstar_cs2 + ( UStar**5 * max(Abs_Coriolis,1.e-20*US%T_to_s) ) )**CS%RH18_mstar_cs2 MStar = MStar_N + MStar_S endif @@ -1804,11 +1804,15 @@ subroutine find_mstar(CS, US, Buoyancy_Flux, UStar, UStar_Mean,& if (CS%answers_2018) then MStar_Conv_Red = 1. - CS%MStar_Convect_coef * (-min(0.0,Buoyancy_Flux) + 1.e-10*US%T_to_s**3*US%m_to_Z**2) / & ( (-min(0.0,Buoyancy_Flux) + 1.e-10*US%T_to_s**3*US%m_to_Z**2) + & - 2.0 *MStar * Ustar**3 / BLD ) + 2.0 *MStar * UStar**3 / BLD ) else MSCR_term1 = -BLD * min(0.0, Buoyancy_Flux) - MSCR_term2 = 2.0*MStar * Ustar**3 - MStar_Conv_Red = ((1.-CS%mstar_convect_coef) * MSCR_term1 + MSCR_term2) / (MSCR_term1 + MSCR_term2) + MSCR_term2 = 2.0*MStar * UStar**3 + if ( abs(MSCR_term2) > 0.0) then + MStar_Conv_Red = ((1.-CS%mstar_convect_coef) * MSCR_term1 + MSCR_term2) / (MSCR_term1 + MSCR_term2) + else + MStar_Conv_Red = 1.-CS%mstar_convect_coef + endif endif !/3. Combine various mstar terms to get final value @@ -1816,15 +1820,15 @@ subroutine find_mstar(CS, US, Buoyancy_Flux, UStar, UStar_Mean,& if (present(Langmuir_Number)) then !### In this call, ustar was previously ustar_mean. Is this change deliberate? - call mstar_Langmuir(CS, US, abs_Coriolis, buoyancy_flux, ustar, BLD, Langmuir_number, mstar, & - mstar_LT, Convect_Langmuir_Number) + call mstar_Langmuir(CS, US, Abs_Coriolis, Buoyancy_Flux, UStar, BLD, Langmuir_Number, MStar, & + MStar_LT, Convect_Langmuir_Number) endif end subroutine Find_Mstar !> This subroutine modifies the Mstar value if the Langmuir number is present -subroutine Mstar_Langmuir(CS, US, abs_Coriolis, buoyancy_flux, ustar, BLD, Langmuir_Number, & - mstar, mstar_LT, Convect_Langmuir_Number) +subroutine Mstar_Langmuir(CS, US, Abs_Coriolis, Buoyancy_Flux, UStar, BLD, Langmuir_Number, & + Mstar, MStar_LT, Convect_Langmuir_Number) type(energetic_PBL_CS), pointer :: CS !< Energetic_PBL control structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, intent(in) :: Abs_Coriolis !< Absolute value of the Coriolis parameter [T-1 ~> s-1] From 68eb36564b79641c8d3fc636358e01121302ed8f Mon Sep 17 00:00:00 2001 From: William Cooke Date: Wed, 6 Nov 2019 15:15:06 -0500 Subject: [PATCH 7/7] add MJHarrison-GFDL salt_flux_add_fix to SPEAR codeset from https://github.com/MJHarrison-GFDL/MOM6/compare/salt_flux_add_fix --- src/core/MOM_forcing_type.F90 | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 7df4213a2f..9e99b15309 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -682,6 +682,12 @@ subroutine extractFluxes1d(G, GV, fluxes, optics, nsw, j, dt, if (do_NSR) Net_salt_rate(i) = (scale * 1. * (1000.0 * fluxes%salt_flux(i,j))) * GV%kg_m2_to_H endif + if (associated(fluxes%salt_flux_added)) then + Net_salt(i) = Net_salt(i) + (scale * dt * (1000.0 * fluxes%salt_flux_added(i,j))) * GV%kg_m2_to_H + !Repeat above code for 'rate' term + if (do_NSR) Net_salt_rate(i) = Net_salt_rate(i) + (scale * 1. * (1000.0 * fluxes%salt_flux_added(i,j))) * GV%kg_m2_to_H + endif + ! Diagnostics follow... if (calculate_diags) then