Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

+Rescale some sea-surface height variables to [Z] #39

Merged
merged 2 commits into from
Dec 14, 2021
Merged

+Rescale some sea-surface height variables to [Z] #39

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
111e441
+Rescale some sea-surface height variables to [Z]
Hallberg-NOAA Dec 11, 2021
eb81922
Merge branch 'dev/gfdl' into rescale_SSH
marshallward Dec 14, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
34 changes: 21 additions & 13 deletions src/core/MOM.F90
View file
Open in desktop
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
Open in desktop
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