Skip to content

Commit

Permalink
+Rescaled the units of fluxes%p_surf
Browse files Browse the repository at this point in the history 
  Rescaled the units of forces%p_surf, fluxes%p_surf, forces%p_surf_full and
fluxes%p_surf_full and related surface pressure variables to [R L2 T-2 ~> Pa]
for expanded dimensional consistency testing.  All answers are bitwise identical,
although there are changes to the rescaled units of elements to two transparent
data types.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Apr 8, 2020
1 parent 6d7dde4 commit 6af725f
Show file tree
Hide file tree
Showing 22 changed files with 125 additions and 125 deletions.
20 changes: 10 additions & 10 deletions config_src/coupled_driver/MOM_surface_forcing_gfdl.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,8 +71,8 @@ module MOM_surface_forcing_gfdl
real :: latent_heat_fusion !< Latent heat of fusion [J kg-1]
real :: latent_heat_vapor !< Latent heat of vaporization [J kg-1]

real :: max_p_surf !< The maximum surface pressure that can be
!! exerted by the atmosphere and floating sea-ice [Pa].
real :: max_p_surf !< The maximum surface pressure that can be exerted by
!! the atmosphere and floating sea-ice [R L2 T-2 ~> Pa].
!! This is needed because the FMS coupling structure
!! does not limit the water that can be frozen out
!! of the ocean and the ice-ocean heat fluxes are
Expand Down Expand Up @@ -548,14 +548,14 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
if (associated(IOB%p)) then
if (CS%max_p_surf >= 0.0) then
do j=js,je ; do i=is,ie
fluxes%p_surf_full(i,j) = G%mask2dT(i,j) * IOB%p(i-i0,j-j0)
fluxes%p_surf_full(i,j) = G%mask2dT(i,j) * US%kg_m3_to_R*US%m_s_to_L_T**2*IOB%p(i-i0,j-j0)
fluxes%p_surf(i,j) = MIN(fluxes%p_surf_full(i,j),CS%max_p_surf)
if (CS%check_no_land_fluxes) &
call check_mask_val_consistency(IOB%p(i-i0,j-j0), G%mask2dT(i,j), i, j, 'p', G)
enddo ; enddo
else
do j=js,je ; do i=is,ie
fluxes%p_surf_full(i,j) = G%mask2dT(i,j) * IOB%p(i-i0,j-j0)
fluxes%p_surf_full(i,j) = G%mask2dT(i,j) * US%kg_m3_to_R*US%m_s_to_L_T**2*IOB%p(i-i0,j-j0)
fluxes%p_surf(i,j) = fluxes%p_surf_full(i,j)
if (CS%check_no_land_fluxes) &
call check_mask_val_consistency(IOB%p(i-i0,j-j0), G%mask2dT(i,j), i, j, 'p', G)
Expand Down Expand Up @@ -673,7 +673,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS, dt_
net_mass_src, & ! A temporary of net mass sources [kg m-2 s-1].
ustar_tmp ! A temporary array of ustar values [Z T-1 ~> m s-1].

real :: I_GEarth ! 1.0 / G_Earth [s2 m-1]
real :: I_GEarth ! Pressure conversion factors times 1.0 / G_Earth [kg m-2 T2 R-1 L-2 ~> s2 m-1]
real :: Kv_rho_ice ! (CS%kv_sea_ice / CS%density_sea_ice) [m5 s-1 kg-1]
real :: mass_ice ! mass of sea ice at a face [kg m-2]
real :: mass_eff ! effective mass of sea ice for rigidity [kg m-2]
Expand Down Expand Up @@ -751,12 +751,12 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS, dt_
if (associated(IOB%p)) then
if (CS%max_p_surf >= 0.0) then
do j=js,je ; do i=is,ie
forces%p_surf_full(i,j) = G%mask2dT(i,j) * IOB%p(i-i0,j-j0)
forces%p_surf_full(i,j) = G%mask2dT(i,j) * US%kg_m3_to_R*US%m_s_to_L_T**2*IOB%p(i-i0,j-j0)
forces%p_surf(i,j) = MIN(forces%p_surf_full(i,j),CS%max_p_surf)
enddo ; enddo
else
do j=js,je ; do i=is,ie
forces%p_surf_full(i,j) = G%mask2dT(i,j) * IOB%p(i-i0,j-j0)
forces%p_surf_full(i,j) = G%mask2dT(i,j) * US%kg_m3_to_R*US%m_s_to_L_T**2*IOB%p(i-i0,j-j0)
forces%p_surf(i,j) = forces%p_surf_full(i,j)
enddo ; enddo
endif
Expand Down Expand Up @@ -837,7 +837,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS, dt_

if (CS%rigid_sea_ice) then
call pass_var(forces%p_surf_full, G%Domain, halo=1)
I_GEarth = 1.0 / CS%G_Earth
I_GEarth = US%RL2_T2_to_Pa / CS%G_Earth
Kv_rho_ice = (CS%kv_sea_ice / CS%density_sea_ice)
do I=is-1,ie ; do j=js,je
mass_ice = min(forces%p_surf_full(i,j), forces%p_surf_full(i+1,j)) * I_GEarth
Expand Down Expand Up @@ -1299,8 +1299,8 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS)
"needed because the FMS coupling structure does not "//&
"limit the water that can be frozen out of the ocean and "//&
"the ice-ocean heat fluxes are treated explicitly. No "//&
"limit is applied if a negative value is used.", units="Pa", &
default=-1.0)
"limit is applied if a negative value is used.", &
units="Pa", default=-1.0, scale=US%kg_m3_to_R*US%m_s_to_L_T**2)
call get_param(param_file, mdl, "RESTORE_SALINITY", CS%restore_salt, &
"If true, the coupled driver will add a globally-balanced "//&
"fresh-water flux that drives sea-surface salinity "//&
Expand Down
18 changes: 11 additions & 7 deletions config_src/mct_driver/mom_surface_forcing_mct.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -68,9 +68,9 @@ module MOM_surface_forcing_mct
real :: latent_heat_fusion !< latent heat of fusion [J kg-1]
real :: latent_heat_vapor !< latent heat of vaporization [J kg-1]

real :: max_p_surf !< maximum surface pressure that can be
!! exerted by the atmosphere and floating sea-ice,
!! [Pa]. This is needed because the FMS coupling
real :: max_p_surf !< The maximum surface pressure that can be exerted by
!! the atmosphere and floating sea-ice [R L2 T-2 ~> Pa].
!! This is needed because the FMS coupling
!! structure does not limit the water that can be
!! frozen out of the ocean and the ice-ocean heat
!! fluxes are treated explicitly.
Expand Down Expand Up @@ -528,11 +528,13 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
if (CS%max_p_surf >= 0.0) then
do j=js,je ; do i=is,ie
fluxes%p_surf_full(i,j) = G%mask2dT(i,j) * IOB%p(i-i0,j-j0)
US%kg_m3_to_R*US%m_s_to_L_T**2*IOB%p(i-i0,j-j0)
fluxes%p_surf(i,j) = MIN(fluxes%p_surf_full(i,j),CS%max_p_surf)
enddo; enddo
else
do j=js,je ; do i=is,ie
fluxes%p_surf_full(i,j) = G%mask2dT(i,j) * IOB%p(i-i0,j-j0)
US%kg_m3_to_R*US%m_s_to_L_T**2*IOB%p(i-i0,j-j0)
fluxes%p_surf(i,j) = fluxes%p_surf_full(i,j)
enddo; enddo
endif
Expand Down Expand Up @@ -621,7 +623,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
real :: tau_mag !< magnitude of the wind stress [R Z L T-2 ~> Pa]
real :: Pa_conversion ! A unit conversion factor from Pa to the internal wind stress units [R Z L T-2 Pa-1 ~> 1]
real :: stress_conversion ! A unit conversion factor from Pa times any stress multiplier [R Z L T-2 Pa-1 ~> 1]
real :: I_GEarth !< 1.0 / G%G_Earth [s2 m-1]
real :: I_GEarth !< Pressure conversion factors times 1.0 / G_Earth [kg m-2 T2 R-1 L-2 ~> s2 m-1]
real :: Kv_rho_ice !< (CS%kv_sea_ice / CS%density_sea_ice) [m5 s-1 kg-1]
real :: mass_ice !< mass of sea ice at a face [kg m-2]
real :: mass_eff !< effective mass of sea ice for rigidity [kg m-2]
Expand Down Expand Up @@ -687,11 +689,13 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
if (CS%max_p_surf >= 0.0) then
do j=js,je ; do i=is,ie
forces%p_surf_full(i,j) = G%mask2dT(i,j) * IOB%p(i-i0,j-j0)
US%kg_m3_to_R*US%m_s_to_L_T**2*IOB%p(i-i0,j-j0)
forces%p_surf(i,j) = MIN(forces%p_surf_full(i,j),CS%max_p_surf)
enddo ; enddo
else
do j=js,je ; do i=is,ie
forces%p_surf_full(i,j) = G%mask2dT(i,j) * IOB%p(i-i0,j-j0)
US%kg_m3_to_R*US%m_s_to_L_T**2*IOB%p(i-i0,j-j0)
forces%p_surf(i,j) = forces%p_surf_full(i,j)
enddo ; enddo
endif
Expand Down Expand Up @@ -845,7 +849,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)

if (CS%rigid_sea_ice) then
call pass_var(forces%p_surf_full, G%Domain, halo=1)
I_GEarth = 1.0 / CS%G_Earth
I_GEarth = US%RL2_T2_to_Pa / CS%g_Earth
Kv_rho_ice = (CS%kv_sea_ice / CS%density_sea_ice)
do I=is-1,ie ; do j=js,je
mass_ice = min(forces%p_surf_full(i,j), forces%p_surf_full(i+1,j)) * I_GEarth
Expand Down Expand Up @@ -1077,8 +1081,8 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt,
"needed because the FMS coupling structure does not "//&
"limit the water that can be frozen out of the ocean and "//&
"the ice-ocean heat fluxes are treated explicitly. No "//&
"limit is applied if a negative value is used.", units="Pa", &
default=-1.0)
"limit is applied if a negative value is used.", &
units="Pa", default=-1.0, scale=US%kg_m3_to_R*US%m_s_to_L_T**2)
call get_param(param_file, mdl, "ADJUST_NET_SRESTORE_TO_ZERO", &
CS%adjust_net_srestore_to_zero, &
"If true, adjusts the salinity restoring seen to zero "//&
Expand Down
22 changes: 11 additions & 11 deletions config_src/nuopc_driver/mom_surface_forcing_nuopc.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,9 +69,9 @@ module MOM_surface_forcing_nuopc
real :: latent_heat_fusion !< latent heat of fusion [J/kg]
real :: latent_heat_vapor !< latent heat of vaporization [J/kg]

real :: max_p_surf !< maximum surface pressure that can be
!! exerted by the atmosphere and floating sea-ice,
!! in Pa. This is needed because the FMS coupling
real :: max_p_surf !< maximum surface pressure that can be exerted by the
!! atmosphere and floating sea-ice [R L2 T-2 ~> Pa].
!! This is needed because the FMS coupling
!! structure does not limit the water that can be
!! frozen out of the ocean and the ice-ocean heat
!! fluxes are treated explicitly.
Expand Down Expand Up @@ -519,12 +519,12 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
if (associated(IOB%p)) then
if (CS%max_p_surf >= 0.0) then
do j=js,je ; do i=is,ie
fluxes%p_surf_full(i,j) = G%mask2dT(i,j) * IOB%p(i-i0,j-j0)
fluxes%p_surf_full(i,j) = G%mask2dT(i,j) * US%kg_m3_to_R*US%m_s_to_L_T**2*IOB%p(i-i0,j-j0)
fluxes%p_surf(i,j) = MIN(fluxes%p_surf_full(i,j),CS%max_p_surf)
enddo; enddo
else
do j=js,je ; do i=is,ie
fluxes%p_surf_full(i,j) = G%mask2dT(i,j) * IOB%p(i-i0,j-j0)
fluxes%p_surf_full(i,j) = G%mask2dT(i,j) * US%kg_m3_to_R*US%m_s_to_L_T**2*IOB%p(i-i0,j-j0)
fluxes%p_surf(i,j) = fluxes%p_surf_full(i,j)
enddo; enddo
endif
Expand Down Expand Up @@ -613,7 +613,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
real :: tau_mag !< magnitude of the wind stress [R Z L T-2 ~> Pa]
real :: Pa_conversion ! A unit conversion factor from Pa to the internal wind stress units [R Z L T-2 Pa-1 ~> 1]
real :: stress_conversion ! A unit conversion factor from Pa times any stress multiplier [R Z L T-2 Pa-1 ~> 1]
real :: I_GEarth !< 1.0 / G_Earth [s2 m-1]
real :: I_GEarth !< Pressure conversion factors times 1.0 / G_Earth [kg m-2 T2 R-1 L-2 ~> s2 m-1]
real :: Kv_rho_ice !< (CS%kv_sea_ice / CS%density_sea_ice) ( m^5/(s*kg) )
real :: mass_ice !< mass of sea ice at a face (kg/m^2)
real :: mass_eff !< effective mass of sea ice for rigidity (kg/m^2)
Expand Down Expand Up @@ -677,12 +677,12 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
if (associated(IOB%p)) then
if (CS%max_p_surf >= 0.0) then
do j=js,je ; do i=is,ie
forces%p_surf_full(i,j) = G%mask2dT(i,j) * IOB%p(i-i0,j-j0)
forces%p_surf_full(i,j) = G%mask2dT(i,j) * US%kg_m3_to_R*US%m_s_to_L_T**2*IOB%p(i-i0,j-j0)
forces%p_surf(i,j) = MIN(forces%p_surf_full(i,j),CS%max_p_surf)
enddo ; enddo
else
do j=js,je ; do i=is,ie
forces%p_surf_full(i,j) = G%mask2dT(i,j) * IOB%p(i-i0,j-j0)
forces%p_surf_full(i,j) = G%mask2dT(i,j) * US%kg_m3_to_R*US%m_s_to_L_T**2*IOB%p(i-i0,j-j0)
forces%p_surf(i,j) = forces%p_surf_full(i,j)
enddo ; enddo
endif
Expand Down Expand Up @@ -840,7 +840,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)

if (CS%rigid_sea_ice) then
call pass_var(forces%p_surf_full, G%Domain, halo=1)
I_GEarth = 1.0 / CS%g_Earth
I_GEarth = US%RL2_T2_to_Pa / CS%g_Earth
Kv_rho_ice = (CS%kv_sea_ice / CS%density_sea_ice)
do I=is-1,ie ; do j=js,je
mass_ice = min(forces%p_surf_full(i,j), forces%p_surf_full(i+1,j)) * I_GEarth
Expand Down Expand Up @@ -1071,8 +1071,8 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt,
"needed because the FMS coupling structure does not "//&
"limit the water that can be frozen out of the ocean and "//&
"the ice-ocean heat fluxes are treated explicitly. No "//&
"limit is applied if a negative value is used.", units="Pa", &
default=-1.0)
"limit is applied if a negative value is used.", &
units="Pa", default=-1.0, scale=US%kg_m3_to_R*US%m_s_to_L_T**2)
call get_param(param_file, mdl, "ADJUST_NET_SRESTORE_TO_ZERO", &
CS%adjust_net_srestore_to_zero, &
"If true, adjusts the salinity restoring seen to zero "//&
Expand Down
25 changes: 12 additions & 13 deletions src/core/MOM.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -264,9 +264,9 @@ module MOM
!! a previous time-step or the ocean restart file.
!! This is only valid when interp_p_surf is true.
real, dimension(:,:), pointer :: &
p_surf_prev => NULL(), & !< surface pressure [Pa] at end previous call to step_MOM
p_surf_begin => NULL(), & !< surface pressure [Pa] at start of step_MOM_dyn_...
p_surf_end => NULL() !< surface pressure [Pa] at end of step_MOM_dyn_...
p_surf_prev => NULL(), & !< surface pressure [R L2 T-2 ~> Pa] at end previous call to step_MOM
p_surf_begin => NULL(), & !< surface pressure [R L2 T-2 ~> Pa] at start of step_MOM_dyn_...
p_surf_end => NULL() !< surface pressure [R L2 T-2 ~> Pa] at end of step_MOM_dyn_...

! Variables needed to reach between start and finish phases of initialization
logical :: write_IC !< If true, then the initial conditions will be written to file
Expand Down Expand Up @@ -473,7 +473,7 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_int_in, CS, &
v => NULL(), & ! v : meridional velocity component [L T-1 ~> m s-1]
h => NULL() ! h : layer thickness [H ~> m or kg m-2]
real, dimension(:,:), pointer :: &
p_surf => NULL() ! A pointer to the ocean surface pressure [Pa].
p_surf => NULL() ! A pointer to the ocean surface pressure [R L2 T-2 ~> Pa].
real :: I_wt_ssh ! The inverse of the time weights [T-1 ~> s-1]

type(time_type) :: Time_local, end_time_thermo, Time_temp
Expand Down Expand Up @@ -878,10 +878,10 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
real, dimension(:,:), pointer :: p_surf_begin !< A pointer (perhaps NULL) to the surface
!! pressure at the beginning of this dynamic
!! step, intent in [Pa].
!! step, intent in [R L2 T-2 ~> Pa].
real, dimension(:,:), pointer :: p_surf_end !< A pointer (perhaps NULL) to the surface
!! pressure at the end of this dynamic step,
!! intent in [Pa].
!! intent in [R L2 T-2 ~> Pa].
real, intent(in) :: dt !< time interval covered by this call [T ~> s].
real, intent(in) :: dt_thermo !< time interval covered by any updates that may
!! span multiple dynamics steps [T ~> s].
Expand Down Expand Up @@ -2449,8 +2449,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
endif

if (CS%interp_p_surf) then
CS%p_surf_prev_set = &
query_initialized(CS%p_surf_prev,"p_surf_prev",restart_CSp)
CS%p_surf_prev_set = query_initialized(CS%p_surf_prev,"p_surf_prev",restart_CSp)

if (CS%p_surf_prev_set) call pass_var(CS%p_surf_prev, G%domain)
endif
Expand Down Expand Up @@ -2683,13 +2682,13 @@ subroutine adjust_ssh_for_p_atm(tv, G, GV, US, ssh, p_atm, use_EOS)
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(:,:), optional, pointer :: p_atm !< atmospheric pressure [Pa]
real, dimension(:,:), optional, pointer :: p_atm !< atmospheric pressure [R L2 T-2 ~> Pa]
logical, optional, intent(in) :: use_EOS !< If true, calculate the density for
!! the SSH correction using the equation of state.

real :: Rho_conv ! The density used to convert surface pressure to
! a corrected effective SSH [R ~> kg m-3].
real :: IgR0 ! The SSH conversion factor from Pa to m [m Pa-1].
real :: IgR0 ! The SSH conversion factor from R L2 T-2 to m [m T2 R-1 L-2 ~> m Pa-1].
logical :: calc_rho
integer :: i, j, is, ie, js, je

Expand All @@ -2701,12 +2700,12 @@ subroutine adjust_ssh_for_p_atm(tv, G, GV, US, ssh, p_atm, use_EOS)
! atmospheric pressure
do j=js,je ; do i=is,ie
if (calc_rho) then
call calculate_density(tv%T(i,j,1), tv%S(i,j,1), p_atm(i,j)/2.0, &
Rho_conv, tv%eqn_of_state, scale=US%kg_m3_to_R)
call calculate_density(tv%T(i,j,1), tv%S(i,j,1), p_atm(i,j)/2.0, Rho_conv, &
tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
else
Rho_conv = GV%Rho0
endif
IgR0 = 1.0 / (Rho_conv * US%R_to_kg_m3*GV%mks_g_Earth)
IgR0 = US%Z_to_m / (Rho_conv * GV%g_Earth)
ssh(i,j) = ssh(i,j) + p_atm(i,j) * IgR0
enddo ; enddo
endif ; endif
Expand Down
4 changes: 2 additions & 2 deletions src/core/MOM_PressureForce.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,10 +59,10 @@ subroutine PressureForce(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm, pbce, e
type(ALE_CS), pointer :: ALE_CSp !< ALE control structure
real, dimension(:,:), &
optional, pointer :: p_atm !< The pressure at the ice-ocean or
!! atmosphere-ocean interface [Pa].
!! atmosphere-ocean interface [R L2 T-2 ~> Pa].
real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
optional, intent(out) :: pbce !< The baroclinic pressure anomaly in each layer
!! due to eta anomalies [m2 s-2 H-1 ~> m s-2 or m4 s-2 kg-1].
!! due to eta anomalies [L2 T-2 H-1 ~> m s-2 or m4 s-2 kg-1].
real, dimension(SZI_(G),SZJ_(G)), &
optional, intent(out) :: eta !< The bottom mass used to calculate PFu and PFv,
!! [H ~> m or kg m-2], with any tidal contributions.
Expand Down
Loading

0 comments on commit 6af725f

Please sign in to comment.