Skip to content

Commit

Permalink
+Replaced optional pres_scale args with US args
Browse files Browse the repository at this point in the history 
  Replaced the remaining pres_scale arguments the various calculate_density
and calc_spec_vol routines in MOM_EOS.F90 with new optional unit_scale_type
arguments.  When the scale and US arguments are present, density is scaled
by the product of the indicated scaling factors.  Calls to these routines
in 11 files were modified accordingly.  All answers are bitwise identical.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Apr 12, 2020
1 parent fb820c1 commit 6948cb7
Show file tree
Hide file tree
Showing 12 changed files with 195 additions and 206 deletions.
4 changes: 2 additions & 2 deletions src/core/MOM_isopycnal_slopes.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -176,7 +176,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, &
S_u(I) = 0.25*((S(i,j,k) + S(i+1,j,k)) + (S(i,j,k-1) + S(i+1,j,k-1)))
enddo
call calculate_density_derivs(T_u, S_u, pres_u, drho_dT_u, drho_dS_u, (is-IsdB+1)-1, ie-is+2, &
tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
tv%eqn_of_state, US=US)
endif

do I=is-1,ie
Expand Down Expand Up @@ -262,7 +262,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, &
S_v(i) = 0.25*((S(i,j,k) + S(i,j+1,k)) + (S(i,j,k-1) + S(i,j+1,k-1)))
enddo
call calculate_density_derivs(T_v, S_v, pres_v, drho_dT_v, drho_dS_v, is, ie-is+1, &
tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
tv%eqn_of_state, US=US)
endif
do i=is,ie
if (use_EOS) then
Expand Down
309 changes: 161 additions & 148 deletions src/equation_of_state/MOM_EOS.F90
View file

Large diffs are not rendered by default.

8 changes: 3 additions & 5 deletions src/initialization/MOM_coord_initialization.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -240,8 +240,7 @@ subroutine set_coord_from_TS_ref(Rlay, g_prime, GV, US, param_file, eqn_of_state
! The uppermost layer's density is set here. Subsequent layers' !
! densities are determined from this value and the g values. !
! T0 = 28.228 ; S0 = 34.5848 ; Pref = P_Ref
call calculate_density(T_ref, S_ref, P_ref, Rlay(1), eqn_of_state, &
scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
call calculate_density(T_ref, S_ref, P_ref, Rlay(1), eqn_of_state, US=US)

! These statements set the layer densities. !
do k=2,nz ; Rlay(k) = Rlay(k-1) + g_prime(k)*(GV%Rho0/GV%g_Earth) ; enddo
Expand Down Expand Up @@ -292,7 +291,7 @@ subroutine set_coord_from_TS_profile(Rlay, g_prime, GV, US, param_file, eqn_of_s
! These statements set the interface reduced gravities. !
g_prime(1) = g_fs
do k=1,nz ; Pref(k) = P_Ref ; enddo
call calculate_density(T0, S0, Pref, Rlay, 1, nz, eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
call calculate_density(T0, S0, Pref, Rlay, 1, nz, eqn_of_state, US=US)
do k=2,nz; g_prime(k) = (GV%g_Earth/(GV%Rho0)) * (Rlay(k) - Rlay(k-1)) ; enddo

call callTree_leave(trim(mdl)//'()')
Expand Down Expand Up @@ -372,8 +371,7 @@ subroutine set_coord_from_TS_range(Rlay, g_prime, GV, US, param_file, eqn_of_sta

g_prime(1) = g_fs
do k=1,nz ; Pref(k) = P_Ref ; enddo
call calculate_density(T0, S0, Pref, Rlay, k_light, nz-k_light+1, eqn_of_state, &
scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
call calculate_density(T0, S0, Pref, Rlay, k_light, nz-k_light+1, eqn_of_state, US=US)
! Extrapolate target densities for the variable density mixed and buffer layers.
do k=k_light-1,1,-1
Rlay(k) = 2.0*Rlay(k+1) - Rlay(k+2)
Expand Down
18 changes: 6 additions & 12 deletions src/initialization/MOM_state_initialization.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1601,10 +1601,8 @@ subroutine initialize_temp_salt_fit(T, S, G, GV, US, param_file, eqn_of_state, P
T0(k) = T_Ref
enddo

call calculate_density(T0(1), S0(1), pres(1), rho_guess(1), eqn_of_state, &
scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, 1, eqn_of_state, &
scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
call calculate_density(T0(1), S0(1), pres(1), rho_guess(1), eqn_of_state, US=US)
call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, 1, eqn_of_state, US=US)

if (fit_salin) then
! A first guess of the layers' temperatures.
Expand All @@ -1613,10 +1611,8 @@ subroutine initialize_temp_salt_fit(T, S, G, GV, US, param_file, eqn_of_state, P
enddo
! Refine the guesses for each layer.
do itt=1,6
call calculate_density(T0, S0, pres, rho_guess, 1, nz, eqn_of_state, &
scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, nz, eqn_of_state, &
scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
call calculate_density(T0, S0, pres, rho_guess, 1, nz, eqn_of_state, US=US)
call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, nz, eqn_of_state, US=US)
do k=1,nz
S0(k) = max(0.0, S0(k) + (GV%Rlay(k) - rho_guess(k)) / drho_dS(k))
enddo
Expand All @@ -1627,10 +1623,8 @@ subroutine initialize_temp_salt_fit(T, S, G, GV, US, param_file, eqn_of_state, P
T0(k) = T0(1) + (GV%Rlay(k) - rho_guess(1)) / drho_dT(1)
enddo
do itt=1,6
call calculate_density(T0, S0, pres, rho_guess, 1, nz, eqn_of_state, &
scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, nz, eqn_of_state, &
scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
call calculate_density(T0, S0, pres, rho_guess, 1, nz, eqn_of_state, US=US)
call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, nz, eqn_of_state, US=US)
do k=1,nz
T0(k) = T0(k) + (GV%Rlay(k) - rho_guess(k)) / drho_dT(k)
enddo
Expand Down
3 changes: 1 addition & 2 deletions src/parameterizations/vertical/MOM_diabatic_aux.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -458,8 +458,7 @@ subroutine insert_brine(h, tv, G, GV, US, fluxes, nkmb, CS, dt, id_brine_lay)
if ((G%mask2dT(i,j) > 0.0) .and. dzbr(i) < brine_dz .and. salt(i) > 0.) then
s_new = S(i,k) + salt(i) / (GV%H_to_RZ * h_2d(i,k))
t0 = T(i,k)
call calculate_density(t0, s_new, tv%P_Ref, R_new, tv%eqn_of_state, &
scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
call calculate_density(t0, s_new, tv%P_Ref, R_new, tv%eqn_of_state, US=US)
if (R_new < 0.5*(Rcv(i,k)+Rcv(i,k+1)) .and. s_new<s_max) then
dzbr(i) = dzbr(i)+h_2d(i,k)
inject_layer(i,j) = min(inject_layer(i,j),real(k))
Expand Down
1 change: 0 additions & 1 deletion src/parameterizations/vertical/MOM_diabatic_driver.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,6 @@ module MOM_diabatic_driver
use MOM_entrain_diffusive, only : entrainment_diffusive, entrain_diffusive_init
use MOM_entrain_diffusive, only : entrain_diffusive_end, entrain_diffusive_CS
use MOM_EOS, only : calculate_density, calculate_TFreeze
use MOM_EOS, only : calculate_specific_vol_derivs
use MOM_error_handler, only : MOM_error, FATAL, WARNING, callTree_showQuery,MOM_mesg
use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint
use MOM_file_parser, only : get_param, log_version, param_file_type, read_param
Expand Down
6 changes: 3 additions & 3 deletions src/parameterizations/vertical/MOM_geothermal.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -199,7 +199,7 @@ subroutine geothermal(h, tv, dt, ea, eb, G, GV, US, CS, halo)

if (nkmb > 0) then
call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_Ref(:), Rcv_BL(:), isj, iej-isj+1, &
tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
tv%eqn_of_state, US=US)
else
Rcv_BL(:) = -1.0
endif
Expand Down Expand Up @@ -245,11 +245,11 @@ subroutine geothermal(h, tv, dt, ea, eb, G, GV, US, CS, halo)
Rcv = 0.0 ; dRcv_dT = 0.0 ! Is this OK?
else
call calculate_density(tv%T(i,j,k), tv%S(i,j,k), tv%P_Ref, &
Rcv, tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
Rcv, tv%eqn_of_state, US=US)
T2(1) = tv%T(i,j,k) ; S2(1) = tv%S(i,j,k)
T2(2) = tv%T(i,j,k_tgt) ; S2(2) = tv%S(i,j,k_tgt)
call calculate_density_derivs(T2(:), S2(:), p_Ref(:), dRcv_dT_, dRcv_dS_, 1, 2, &
tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
tv%eqn_of_state, US=US)
dRcv_dT = 0.5*(dRcv_dT_(1) + dRcv_dT_(2))
endif

Expand Down
2 changes: 1 addition & 1 deletion src/parameterizations/vertical/MOM_kappa_shear.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -911,7 +911,7 @@ subroutine kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, &
Sal_int(K) = 0.5*(Sal(k-1) + Sal(k))
enddo
call calculate_density_derivs(T_int, Sal_int, pressure, dbuoy_dT, dbuoy_dS, 2, nzc-1, &
tv%eqn_of_state, scale=-g_R0*US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
tv%eqn_of_state, US=US, scale=-g_R0)
else
do K=1,nzc+1 ; dbuoy_dT(K) = -g_R0 ; dbuoy_dS(K) = 0.0 ; enddo
endif
Expand Down
4 changes: 2 additions & 2 deletions src/parameterizations/vertical/MOM_set_viscosity.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1397,7 +1397,7 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS, symmetri

if (use_EOS) then
call calculate_density_derivs(T_EOS, S_EOS, forces%p_surf(:,j), dR_dT, dR_dS, &
Isq-G%IsdB+1, Ieq-Isq+1, tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_PA)
Isq-G%IsdB+1, Ieq-Isq+1, tv%eqn_of_state, US=US)
endif

do I=Isq,Ieq ; if (do_i(I)) then
Expand Down Expand Up @@ -1634,7 +1634,7 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS, symmetri

if (use_EOS) then
call calculate_density_derivs(T_EOS, S_EOS, forces%p_surf(:,j), dR_dT, dR_dS, &
is-G%IsdB+1, ie-is+1, tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_PA)
is-G%IsdB+1, ie-is+1, tv%eqn_of_state, US=US)
endif

do i=is,ie ; if (do_i(i)) then
Expand Down
10 changes: 4 additions & 6 deletions src/tracer/MOM_tracer_Z_init.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -801,10 +801,9 @@ subroutine determine_temperature(temp, salt, R_tgt, p_ref, niter, land_fill, h,
adjust_salt = .true.
iter_loop: do itt = 1,niter
do k=1, nz
call calculate_density(T(:,k), S(:,k), press, rho(:,k), 1, nx, eos, &
scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
call calculate_density(T(:,k), S(:,k), press, rho(:,k), 1, nx, eos, US=US)
call calculate_density_derivs(T(:,k), S(:,k), press, drho_dT(:,k), drho_dS(:,k), 1, nx, &
eos, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
eos, US=US)
enddo
do k=k_start,nz ; do i=1,nx

Expand Down Expand Up @@ -832,10 +831,9 @@ subroutine determine_temperature(temp, salt, R_tgt, p_ref, niter, land_fill, h,

if (adjust_salt .and. old_fit) then ; do itt = 1,niter
do k=1, nz
call calculate_density(T(:,k), S(:,k), press, rho(:,k), 1, nx, eos, &
scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
call calculate_density(T(:,k), S(:,k), press, rho(:,k), 1, nx, eos, US=US)
call calculate_density_derivs(T(:,k), S(:,k), press, drho_dT(:,k), drho_dS(:,k), 1, nx, &
eos, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
eos, US=US)
enddo
do k=k_start,nz ; do i=1,nx
! if (abs(rho(i,k)-R_tgt(k))>tol_rho .and. hin(i,k)>h_massless .and. abs(T(i,k)-land_fill) < epsln ) then
Expand Down
12 changes: 4 additions & 8 deletions src/user/DOME_initialization.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -359,17 +359,13 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, US, param_file, tr_Reg)
! target density and a salinity of 35 psu. This code is taken from
! USER_initialize_temp_sal.
pres(:) = tv%P_Ref ; S0(:) = 35.0 ; T0(1) = 25.0
call calculate_density(T0(1), S0(1), pres(1), rho_guess(1), tv%eqn_of_state, &
scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, 1, tv%eqn_of_state, &
scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
call calculate_density(T0(1), S0(1), pres(1), rho_guess(1), tv%eqn_of_state, US=US)
call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, 1, tv%eqn_of_state, US=US)

do k=1,nz ; T0(k) = T0(1) + (GV%Rlay(k)-rho_guess(1)) / drho_dT(1) ; enddo
do itt=1,6
call calculate_density(T0, S0, pres, rho_guess, 1, nz, tv%eqn_of_state, &
scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, nz, tv%eqn_of_state, &
scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
call calculate_density(T0, S0, pres, rho_guess, 1, nz, tv%eqn_of_state, US=US)
call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, nz, tv%eqn_of_state, US=US)
do k=1,nz ; T0(k) = T0(k) + (GV%Rlay(k)-rho_guess(k)) / drho_dT(k) ; enddo
enddo

Expand Down
24 changes: 8 additions & 16 deletions src/user/benchmark_initialization.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -152,10 +152,8 @@ subroutine benchmark_initialize_thickness(h, G, GV, US, param_file, eqn_of_state
pres(k) = P_Ref ; S0(k) = 35.0
enddo
T0(k1) = 29.0
call calculate_density(T0(k1), S0(k1), pres(k1), rho_guess(k1), eqn_of_state, &
scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, k1, 1, eqn_of_state, &
scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
call calculate_density(T0(k1), S0(k1), pres(k1), rho_guess(k1), eqn_of_state, US=US)
call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, k1, 1, eqn_of_state, US=US)

! A first guess of the layers' temperatures.
do k=1,nz
Expand All @@ -164,10 +162,8 @@ subroutine benchmark_initialize_thickness(h, G, GV, US, param_file, eqn_of_state

! Refine the guesses for each layer.
do itt=1,6
call calculate_density(T0, S0, pres, rho_guess, 1, nz, eqn_of_state, &
scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, nz, eqn_of_state, &
scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
call calculate_density(T0, S0, pres, rho_guess, 1, nz, eqn_of_state, US=US)
call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, nz, eqn_of_state, US=US)
do k=1,nz
T0(k) = T0(k) + (GV%Rlay(k) - rho_guess(k)) / drho_dT(k)
enddo
Expand Down Expand Up @@ -261,10 +257,8 @@ subroutine benchmark_init_temperature_salinity(T, S, G, GV, US, param_file, &
enddo

T0(k1) = 29.0
call calculate_density(T0(k1), S0(k1), pres(k1), rho_guess(k1), eqn_of_state, &
scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, k1, 1, eqn_of_state, &
scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
call calculate_density(T0(k1), S0(k1), pres(k1), rho_guess(k1), eqn_of_state, US=US)
call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, k1, 1, eqn_of_state, US=US)

! A first guess of the layers' temperatures. !
do k=1,nz
Expand All @@ -273,10 +267,8 @@ subroutine benchmark_init_temperature_salinity(T, S, G, GV, US, param_file, &

! Refine the guesses for each layer. !
do itt = 1,6
call calculate_density(T0,S0,pres,rho_guess,1,nz,eqn_of_state, &
scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
call calculate_density_derivs(T0,S0,pres,drho_dT,drho_dS,1,nz,eqn_of_state, &
scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
call calculate_density(T0, S0, pres, rho_guess, 1, nz, eqn_of_state, US=US)
call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, nz, eqn_of_state, US=US)
do k=1,nz
T0(k) = T0(k) + (GV%Rlay(k) - rho_guess(k)) / drho_dT(k)
enddo
Expand Down

0 comments on commit 6948cb7

Please sign in to comment.