diff --git a/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 b/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 index cbbb2261d0..1ba3484ef9 100644 --- a/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 +++ b/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 @@ -690,13 +690,8 @@ subroutine ocean_model_restart(OS, timestamp, restartname, num_rest_files) "restart files can only be created after the buoyancy forcing is applied.") if (present(restartname)) then - if (present(num_rest_files)) then - call save_restart(OS%dirs%restart_output_dir, OS%Time, OS%grid, & - OS%restart_CSp, GV=OS%GV, filename=restartname, num_rest_files=num_rest_files) - else - call save_restart(OS%dirs%restart_output_dir, OS%Time, OS%grid, & - OS%restart_CSp, GV=OS%GV, filename=restartname) - endif + call save_restart(OS%dirs%restart_output_dir, OS%Time, OS%grid, & + OS%restart_CSp, GV=OS%GV, filename=restartname, num_rest_files=num_rest_files) call forcing_save_restart(OS%forcing_CSp, OS%grid, OS%Time, & OS%dirs%restart_output_dir) ! Is this needed? if (OS%use_ice_shelf) then diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 82be08100e..819131e51d 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -106,7 +106,8 @@ module MOM_diagnostics KE_dia => NULL() !< KE source from diapycnal diffusion [H L2 T-3 ~> m3 s-3] !>@{ Diagnostic IDs - integer :: id_u = -1, id_v = -1, id_h = -1 + integer :: id_u = -1, id_v = -1, id_h = -1 + integer :: id_usq = -1, id_vsq = -1, id_uv = -1 integer :: id_e = -1, id_e_D = -1 integer :: id_du_dt = -1, id_dv_dt = -1 integer :: id_col_ht = -1, id_dh_dt = -1 @@ -223,6 +224,10 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & !! calculating interface heights [H ~> m or kg m-2]. ! Local variables + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: usq ! squared eastward velocity [L2 T-2 ~> m2 s-2] + real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: vsq ! squared northward velocity [L2 T-2 ~> m2 s-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: uv ! u x v at h-points [L2 T-2 ~> m2 s-2] + integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state integer i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb @@ -289,6 +294,28 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & if (CS%id_h > 0) call post_data(CS%id_h, h, CS%diag) + if (CS%id_usq > 0) then + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + usq(I,j,k) = u(I,j,k) * u(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_usq, usq, CS%diag) + endif + + if (CS%id_vsq > 0) then + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + vsq(i,J,k) = v(i,J,k) * v(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_vsq, vsq, CS%diag) + endif + + if (CS%id_uv > 0) then + do k=1,nz ; do j=js,je ; do i=is,ie + uv(i,j,k) = (0.5*(u(I-1,j,k) + u(I,j,k))) * & + (0.5*(v(i,J-1,k) + v(i,J,k))) + enddo ; enddo ; enddo + call post_data(CS%id_uv, uv, CS%diag) + endif + if (associated(CS%e)) then call find_eta(h, tv, G, GV, US, CS%e, eta_bt) if (CS%id_e > 0) call post_data(CS%id_e, CS%e, CS%diag) @@ -1577,6 +1604,13 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag CS%id_v = register_diag_field('ocean_model', 'v', diag%axesCvL, Time, & 'Meridional velocity', 'm s-1', conversion=US%L_T_to_m_s, cmor_field_name='vo', & cmor_standard_name='sea_water_y_velocity', cmor_long_name='Sea Water Y Velocity') + CS%id_usq = register_diag_field('ocean_model', 'usq', diag%axesCuL, Time, & + 'Zonal velocity squared', 'm2 s-2', conversion=US%L_T_to_m_s**2) + CS%id_vsq = register_diag_field('ocean_model', 'vsq', diag%axesCvL, Time, & + 'Meridional velocity squared', 'm2 s-2', conversion=US%L_T_to_m_s**2) + CS%id_uv = register_diag_field('ocean_model', 'uv', diag%axesTL, Time, & + 'Product between zonal and meridional velocities at h-points', 'm2 s-2', & + conversion=US%L_T_to_m_s**2) CS%id_h = register_diag_field('ocean_model', 'h', diag%axesTL, Time, & 'Layer Thickness', thickness_units, v_extensive=.true., conversion=convert_H) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index a60d60bb9d..953cc6d838 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -515,7 +515,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, enddo ; enddo ! Components for the shearing strain - do J=js-2,Jeq+2 ; do I=is-2,Ieq+2 + do J=Jsq-2,Jeq+2 ; do I=Isq-2,Ieq+2 dvdx(I,J) = CS%DY_dxBu(I,J)*(v(i+1,J,k)*G%IdyCv(i+1,J) - v(i,J,k)*G%IdyCv(i,J)) dudy(I,J) = CS%DX_dyBu(I,J)*(u(I,j+1,k)*G%IdxCu(I,j+1) - u(I,j,k)*G%IdxCu(I,j)) enddo ; enddo @@ -636,7 +636,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, elseif (OBC%segment(n)%direction == OBC_DIRECTION_S) then if ((J >= js-1) .and. (J <= je+1)) then do I = max(Isq-1,OBC%segment(n)%HI%isd), min(Ieq+1,OBC%segment(n)%HI%ied) - h_u(I,j) = h_u(i,j+1) + h_u(I,j) = h_u(I,j+1) enddo endif elseif (OBC%segment(n)%direction == OBC_DIRECTION_E) then @@ -694,11 +694,11 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, ! Vorticity if (CS%no_slip) then - do J=js-2,Jeq+2 ; do I=is-2,Ieq+2 + do J=Jsq-2,Jeq+2 ; do I=Isq-2,Ieq+2 vort_xy(I,J) = (2.0-G%mask2dBu(I,J)) * ( dvdx(I,J) - dudy(I,J) ) enddo ; enddo else - do J=js-2,Jeq+2 ; do I=is-2,Ieq+2 + do J=Jsq-2,Jeq+2 ; do I=Isq-2,Ieq+2 vort_xy(I,J) = G%mask2dBu(I,J) * ( dvdx(I,J) - dudy(I,J) ) enddo ; enddo endif @@ -711,22 +711,23 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then ! Vorticity gradient - do J=js-2,Jeq+2 ; do i=is-1,Ieq+2 + do J=Jsq-1,Jeq+1 ; do i=Isq-1,Ieq+2 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+2 ; do I=is-2,Ieq+2 + do j=Jsq-1,Jeq+2 ; do I=Isq-1,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 ! Laplacian of vorticity do J=Jsq-1,Jeq+1 ; do I=Isq-1,Ieq+1 - DY_dxCv = G%dyCv(i,J) * G%IdxCv(i,J) - DX_dyCu = G%dyCu(I,j) * G%IdyCu(I,j) - Del2vort_q(I,J) = DY_dxCv * (vort_xy_dx(i+1,J) * G%IdyT(i+1,j) - vort_xy_dx(i,J) * G%IdyT(i,j)) + & - DX_dyCu * (vort_xy_dy(I,j+1) * G%IdyT(i,j+1) - vort_xy_dy(I,j) * G%IdyT(i,j)) + DY_dxBu = G%dyBu(I,J) * G%IdxBu(I,J) + DX_dyBu = G%dxBu(I,J) * G%IdyBu(I,J) + + Del2vort_q(I,J) = DY_dxBu * (vort_xy_dx(i+1,J) * G%IdyCv(i+1,J) - vort_xy_dx(i,J) * G%IdyCv(i,J)) + & + DX_dyBu * (vort_xy_dy(I,j+1) * G%IdyCu(I,j+1) - vort_xy_dy(I,j) * G%IdyCu(I,j)) enddo ; enddo do J=Jsq,Jeq+1 ; do I=Isq,Ieq+1 Del2vort_h(i,j) = 0.25*(Del2vort_q(I,J) + Del2vort_q(I-1,J) + Del2vort_q(I,J-1) + Del2vort_q(I-1,J-1)) @@ -1086,12 +1087,12 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, 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)) ) + 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%Re_Ah > 0.0) then - KE = 0.125*((u(I,j,k)+u(I-1,j,k))**2 + (v(i,J,k)+v(i,J-1,k))**2) + KE = 0.125*((u(I,j,k)+u(I,j+1,k))**2 + (v(i,J,k)+v(i+1,J,k))**2) Ah = sqrt(KE) * CS%Re_Ah_const_xy(i,j) endif @@ -1193,7 +1194,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, CS%dy2h(i+1,j)*str_xx(i+1,j)) + & G%IdxCu(I,j)*(CS%dx2q(I,J-1)*str_xy(I,J-1) - & CS%dx2q(I,J) *str_xy(I,J))) * & - G%IareaCu(I,j)) / (h_u(i,j) + h_neglect) + G%IareaCu(I,j)) / (h_u(I,j) + h_neglect) enddo ; enddo if (apply_OBC) then @@ -1922,7 +1923,7 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) endif endif if (CS%Leith_Ah) then - CS%biharm6_const_xy(i,j) = Leith_bi_const * (grid_sp_q3 * grid_sp_q3) + CS%biharm6_const_xy(I,J) = Leith_bi_const * (grid_sp_q3 * grid_sp_q3) endif CS%Ah_bg_xy(I,J) = MAX(Ah, Ah_vel_scale * grid_sp_q2 * sqrt(grid_sp_q2)) if (CS%Re_Ah > 0.0) CS%Re_Ah_const_xy(i,j) = grid_sp_q3 / CS%Re_Ah diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90 index 951170c039..708d6a7f46 100644 --- a/src/parameterizations/vertical/MOM_tidal_mixing.F90 +++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90 @@ -783,8 +783,8 @@ subroutine calculate_CVMix_tidal(h, j, G, GV, US, CS, N2_int, Kd_lay, Kd_int, Kv ! XXX: Temporary de-scaling of N2_int(i,:) into a temporary variable - do k=1,G%ke+1 - N2_int_i(k) = US%s_to_T**2 * N2_int(i,k) + do K=1,G%ke+1 + N2_int_i(K) = US%s_to_T**2 * N2_int(i,K) enddo call CVMix_coeffs_tidal( Mdiff_out = Kv_tidal, & @@ -803,14 +803,14 @@ subroutine calculate_CVMix_tidal(h, j, G, GV, US, CS, N2_int, Kd_lay, Kd_int, Kv Kd_lay(i,j,k) = Kd_lay(i,j,k) + 0.5 * US%m2_s_to_Z2_T * (Kd_tidal(k) + Kd_tidal(k+1)) enddo if (present(Kd_int)) then - do k=1,G%ke+1 - Kd_int(i,j,k) = Kd_int(i,j,k) + (US%m2_s_to_Z2_T * Kd_tidal(k)) + do K=1,G%ke+1 + Kd_int(i,j,K) = Kd_int(i,j,K) + (US%m2_s_to_Z2_T * Kd_tidal(K)) enddo endif ! Update viscosity with the proper unit conversion. if (associated(Kv)) then - do k=1,G%ke+1 - Kv(i,j,k) = Kv(i,j,k) + US%m2_s_to_Z2_T * Kv_tidal(k) ! Rescale from m2 s-1 to Z2 T-1. + do K=1,G%ke+1 + Kv(i,j,K) = Kv(i,j,K) + US%m2_s_to_Z2_T * Kv_tidal(K) ! Rescale from m2 s-1 to Z2 T-1. enddo endif @@ -903,15 +903,15 @@ subroutine calculate_CVMix_tidal(h, j, G, GV, US, CS, N2_int, Kd_lay, Kd_int, Kv Kd_lay(i,j,k) = Kd_lay(i,j,k) + 0.5 * US%m2_s_to_Z2_T * (Kd_tidal(k) + Kd_tidal(k+1)) enddo if (present(Kd_int)) then - do k=1,G%ke+1 - Kd_int(i,j,k) = Kd_int(i,j,k) + (US%m2_s_to_Z2_T * Kd_tidal(k)) + do K=1,G%ke+1 + Kd_int(i,j,K) = Kd_int(i,j,K) + (US%m2_s_to_Z2_T * Kd_tidal(K)) enddo endif ! Update viscosity if (associated(Kv)) then - do k=1,G%ke+1 - Kv(i,j,k) = Kv(i,j,k) + US%m2_s_to_Z2_T * Kv_tidal(k) ! Rescale from m2 s-1 to Z2 T-1. + do K=1,G%ke+1 + Kv(i,j,K) = Kv(i,j,K) + US%m2_s_to_Z2_T * Kv_tidal(K) ! Rescale from m2 s-1 to Z2 T-1. enddo endif