diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 0e036d9d8f..db114ac3fa 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -188,10 +188,10 @@ module MOM vh, & !< vh = v * h * dx at v grid points [H L2 T-1 ~> m3 s-1 or kg s-1] vhtr !< accumulated meridional thickness fluxes to advect tracers [H L2 ~> m3 or kg] real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: ssh_rint - !< A running time integral of the sea surface height [T m ~> s m]. + !< A running time integral of the sea surface height [T Z ~> s m]. real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: ave_ssh_ibc !< time-averaged (over a forcing time step) sea surface height - !! with a correction for the inverse barometer [m] + !! with a correction for the inverse barometer [Z ~> m] real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: eta_av_bc !< free surface height or column mass time averaged over the last !! baroclinic dynamics time step [H ~> m or kg m-2] @@ -515,7 +515,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS logical :: therm_reset ! If true, reset running sums of thermodynamic quantities. real :: cycle_time ! The length of the coupled time-stepping cycle [T ~> s]. real, dimension(SZI_(CS%G),SZJ_(CS%G)) :: & - ssh ! sea surface height, which may be based on eta_av [m] + ssh ! sea surface height, which may be based on eta_av [Z ~> m] real, dimension(:,:,:), pointer :: & u => NULL(), & ! u : zonal velocity component [L T-1 ~> m s-1] @@ -868,7 +868,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS ! Determining the time-average sea surface height is part of the algorithm. ! This may be eta_av if Boussinesq, or need to be diagnosed if not. CS%time_in_cycle = CS%time_in_cycle + dt - call find_eta(h, CS%tv, G, GV, US, ssh, CS%eta_av_bc, eta_to_m=1.0, dZref=G%Z_ref) + call find_eta(h, CS%tv, G, GV, US, ssh, CS%eta_av_bc, dZref=G%Z_ref) do j=js,je ; do i=is,ie CS%ssh_rint(i,j) = CS%ssh_rint(i,j) + dt*ssh(i,j) enddo ; enddo @@ -2867,11 +2867,18 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & endif endif - if (.not.query_initialized(CS%ave_ssh_ibc,"ave_ssh",restart_CSp)) then + if (query_initialized(CS%ave_ssh_ibc,"ave_ssh",restart_CSp)) then + if ((US%m_to_Z_restart /= 0.0) .and. (US%m_to_Z /= US%m_to_Z_restart) ) then + Z_rescale = US%m_to_Z / US%m_to_Z_restart + do j=js,je ; do i=is,ie + CS%ave_ssh_ibc(i,j) = Z_rescale * CS%ave_ssh_ibc(i,j) + enddo ; enddo + endif + else if (CS%split) then - call find_eta(CS%h, CS%tv, G, GV, US, CS%ave_ssh_ibc, eta, eta_to_m=1.0, dZref=G%Z_ref) + call find_eta(CS%h, CS%tv, G, GV, US, CS%ave_ssh_ibc, eta, dZref=G%Z_ref) else - call find_eta(CS%h, CS%tv, G, GV, US, CS%ave_ssh_ibc, eta_to_m=1.0, dZref=G%Z_ref) + call find_eta(CS%h, CS%tv, G, GV, US, CS%ave_ssh_ibc, dZref=G%Z_ref) endif endif if (CS%split) deallocate(eta) @@ -2977,7 +2984,7 @@ subroutine register_diags(Time, G, GV, US, IDs, diag) 'Layer Thickness after the dynamics update', thickness_units, conversion=GV%H_to_MKS, & v_extensive=.true.) IDs%id_ssh_inst = register_diag_field('ocean_model', 'SSH_inst', diag%axesT1, & - Time, 'Instantaneous Sea Surface Height', 'm') + Time, 'Instantaneous Sea Surface Height', 'm', conversion=US%Z_to_m) end subroutine register_diags !> Set up CPU clock IDs for timing various subroutines. @@ -3097,14 +3104,14 @@ subroutine adjust_ssh_for_p_atm(tv, G, GV, US, ssh, p_atm, use_EOS) type(ocean_grid_type), intent(in) :: 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 - real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: ssh !< time mean surface height [m] + real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: ssh !< time mean surface height [Z ~> m] real, dimension(:,:), pointer :: p_atm !< Ocean surface pressure [R L2 T-2 ~> Pa] logical, intent(in) :: use_EOS !< If true, calculate the density for !! the SSH correction using the equation of state. real :: Rho_conv(SZI_(G)) ! The density used to convert surface pressure to ! a corrected effective SSH [R ~> kg m-3]. - real :: IgR0 ! The SSH conversion factor from R L2 T-2 to m [m T2 R-1 L-2 ~> m Pa-1]. + real :: IgR0 ! The SSH conversion factor from R L2 T-2 to Z [Z T2 R-1 L-2 ~> m Pa-1]. logical :: calc_rho integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state integer :: i, j, is, ie, js, je @@ -3119,12 +3126,13 @@ subroutine adjust_ssh_for_p_atm(tv, G, GV, US, ssh, p_atm, use_EOS) call calculate_density(tv%T(:,j,1), tv%S(:,j,1), 0.5*p_atm(:,j), Rho_conv, & tv%eqn_of_state, EOSdom) do i=is,ie - IgR0 = US%Z_to_m / (Rho_conv(i) * GV%g_Earth) + IgR0 = 1.0 / (Rho_conv(i) * GV%g_Earth) ssh(i,j) = ssh(i,j) + p_atm(i,j) * IgR0 enddo else + IgR0 = 1.0 / (GV%Rho0 * GV%g_Earth) do i=is,ie - ssh(i,j) = ssh(i,j) + p_atm(i,j) * (US%Z_to_m / (GV%Rho0 * GV%g_Earth)) + ssh(i,j) = ssh(i,j) + p_atm(i,j) * IgR0 enddo endif enddo @@ -3209,7 +3217,7 @@ subroutine extract_surface_state(CS, sfc_state_in) sfc_state%S_is_absS = CS%tv%S_is_absS do j=js,je ; do i=is,ie - sfc_state%sea_lev(i,j) = US%m_to_Z*CS%ave_ssh_ibc(i,j) + sfc_state%sea_lev(i,j) = CS%ave_ssh_ibc(i,j) enddo ; enddo if (allocated(sfc_state%frazil) .and. associated(CS%tv%frazil)) then ; do j=js,je ; do i=is,ie diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index e06cb235c4..8d667503d7 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -209,18 +209,23 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & ! tmp array for surface properties - real :: surface_field(SZI_(G),SZJ_(G)) + real :: surface_field(SZI_(G),SZJ_(G)) ! The surface temperature or salinity [degC] or [ppt] real :: pressure_1d(SZI_(G)) ! Temporary array for pressure when calling EOS [R L2 T-2 ~> Pa] - real :: wt, wt_p - + real :: wt, wt_p ! The fractional weights of two successive values when interpolating from + ! a list [nondim], scaled so that wt + wt_p = 1. real :: f2_h ! Squared Coriolis parameter at to h-points [T-2 ~> s-2] real :: mag_beta ! Magnitude of the gradient of f [T-1 L-1 ~> s-1 m-1] real :: absurdly_small_freq2 ! Frequency squared used to avoid division by 0 [T-2 ~> s-2] integer :: k_list - real, dimension(SZK_(GV)) :: temp_layer_ave, salt_layer_ave - real :: thetaoga, soga, masso, tosga, sosga + real, dimension(SZK_(GV)) :: temp_layer_ave ! The average temperature in a layer [degC] + real, dimension(SZK_(GV)) :: salt_layer_ave ! The average salinity in a layer [degC] + real :: thetaoga ! The volume mean potential temperature [degC] + real :: soga ! The volume mean ocean salinity [ppt] + real :: masso ! The total mass of the ocean [kg] + real :: tosga ! The area mean sea surface temperature [degC] + real :: sosga ! The area mean sea surface salinity [ppt] is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB @@ -437,7 +442,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & do j=js,je ; do i=is,ie surface_field(i,j) = tv%T(i,j,1) enddo ; enddo - tosga = global_area_mean(surface_field, G) + tosga = global_area_mean(tv%T(:,:,1), G) call post_data(CS%id_tosga, tosga, CS%diag) endif @@ -1240,7 +1245,8 @@ subroutine post_surface_dyn_diags(IDs, G, diag, sfc_state, ssh) type(diag_ctrl), intent(in) :: diag !< regulates diagnostic output type(surface), intent(in) :: sfc_state !< structure describing the ocean surface state real, dimension(SZI_(G),SZJ_(G)), & - intent(in) :: ssh !< Time mean surface height without corrections for ice displacement [m] + intent(in) :: ssh !< Time mean surface height without corrections + !! for ice displacement [Z ~> m] ! Local variables real, dimension(SZI_(G),SZJ_(G)) :: speed ! The surface speed [L T-1 ~> m s-1] @@ -1280,23 +1286,25 @@ subroutine post_surface_thermo_diags(IDs, G, GV, US, diag, dt_int, sfc_state, tv real, intent(in) :: dt_int !< total time step associated with these diagnostics [T ~> s]. type(surface), intent(in) :: sfc_state !< structure describing the ocean surface state type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables - real, dimension(SZI_(G),SZJ_(G)), & - intent(in) :: ssh !< Time mean surface height without corrections for ice displacement [m] + real, dimension(SZI_(G),SZJ_(G)), intent(in) :: ssh !< Time mean surface height without corrections + !! for ice displacement [Z ~> m] real, dimension(SZI_(G),SZJ_(G)), intent(in) :: ssh_ibc !< Time mean surface height with corrections - !! for ice displacement and the inverse barometer [m] + !! for ice displacement and the inverse barometer [Z ~> m] real, dimension(SZI_(G),SZJ_(G)) :: work_2d ! A 2-d work array real, dimension(SZI_(G),SZJ_(G)) :: & zos ! dynamic sea lev (zero area mean) from inverse-barometer adjusted ssh [m] real :: I_time_int ! The inverse of the time interval [T-1 ~> s-1]. - real :: zos_area_mean, volo, ssh_ga + real :: zos_area_mean ! Global area mean sea surface height [m] + real :: volo ! Total volume of the ocean [m3] + real :: ssh_ga ! Global ocean area weighted mean sea seaface height [m] integer :: i, j, is, ie, js, je is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ! area mean SSH if (IDs%id_ssh_ga > 0) then - ssh_ga = global_area_mean(ssh, G) + ssh_ga = global_area_mean(ssh, G, scale=US%Z_to_m) call post_data(IDs%id_ssh_ga, ssh_ga, diag) endif @@ -1306,7 +1314,7 @@ subroutine post_surface_thermo_diags(IDs, G, GV, US, diag, dt_int, sfc_state, tv if (IDs%id_zos > 0 .or. IDs%id_zossq > 0) then zos(:,:) = 0.0 do j=js,je ; do i=is,ie - zos(i,j) = ssh_ibc(i,j) + zos(i,j) = US%Z_to_m*ssh_ibc(i,j) enddo ; enddo zos_area_mean = global_area_mean(zos, G) do j=js,je ; do i=is,ie @@ -1324,9 +1332,9 @@ subroutine post_surface_thermo_diags(IDs, G, GV, US, diag, dt_int, sfc_state, tv ! post total volume of the liquid ocean if (IDs%id_volo > 0) then do j=js,je ; do i=is,ie - work_2d(i,j) = G%mask2dT(i,j)*(ssh(i,j) + US%Z_to_m*G%bathyT(i,j)) + work_2d(i,j) = G%mask2dT(i,j) * (ssh(i,j) + G%bathyT(i,j)) enddo ; enddo - volo = global_area_integral(work_2d, G) + volo = global_area_integral(work_2d, G, scale=US%Z_to_m) call post_data(IDs%id_volo, volo, diag) endif @@ -1841,7 +1849,7 @@ subroutine register_surface_diags(Time, G, US, IDs, diag, tv) standard_name='square_of_sea_surface_height_above_geoid', & long_name='Square of sea surface height above geoid', units='m2') IDs%id_ssh = register_diag_field('ocean_model', 'SSH', diag%axesT1, Time, & - 'Sea Surface Height', 'm') + 'Sea Surface Height', 'm', conversion=US%Z_to_m) IDs%id_ssh_ga = register_scalar_field('ocean_model', 'ssh_ga', Time, diag,& long_name='Area averaged sea surface height', units='m', & standard_name='area_averaged_sea_surface_height')