Skip to content

Commit

Permalink
+Rescale some sea-surface height variables to [Z]
Browse files Browse the repository at this point in the history 
  Rescaled the units of four internal sea surface height or related variables in
MOM.F90 and three sea surface height arguments to post_surface_dyn_diags and
post_surface_thermo_diags from [m] to [Z ~> m].  Also added a few comments
describing variables in MOM_diagnostics.F90.  All answers, diagnostics, and
output are bitwise identical.
  • Loading branch information
@Hallberg-NOAA @marshallward
Hallberg-NOAA authored and marshallward committed Dec 14, 2021
1 parent 6feadd0 commit bbb9753
Show file tree
Hide file tree
Showing 2 changed files with 45 additions and 29 deletions.
34 changes: 21 additions & 13 deletions src/core/MOM.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down Expand Up @@ -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]
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down
40 changes: 24 additions & 16 deletions src/diagnostics/MOM_diagnostics.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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]
Expand Down Expand Up @@ -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

Expand All @@ -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
Expand All @@ -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

Expand Down Expand Up @@ -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')
Expand Down

0 comments on commit bbb9753

Please sign in to comment.