From acf23a413a5b8c13f98ff387c9a3f0654414e30f Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 17 Apr 2020 13:15:36 -0400 Subject: [PATCH] Use the 'dom=' interface to calculate_density Use the new variant of calculate_density with the 'dom' argument or no array extent argument in calls in 15 files. All answers are bitwise identical. --- src/ALE/MOM_regridding.F90 | 6 +-- src/ALE/coord_adapt.F90 | 10 ++-- src/core/MOM_PressureForce_Montgomery.F90 | 54 ++++++++++--------- src/core/MOM_PressureForce_analytic_FV.F90 | 30 ++++++----- src/core/MOM_PressureForce_blocked_AFV.F90 | 30 ++++++----- src/core/MOM_isopycnal_slopes.F90 | 21 +++++--- src/diagnostics/MOM_wave_speed.F90 | 8 +-- src/diagnostics/MOM_wave_structure.F90 | 4 +- .../MOM_coord_initialization.F90 | 4 +- .../MOM_state_initialization.F90 | 8 +-- .../lateral/MOM_thickness_diffuse.F90 | 24 ++++++--- .../vertical/MOM_geothermal.F90 | 8 +-- .../vertical/MOM_set_viscosity.F90 | 7 +-- src/user/DOME_initialization.F90 | 6 +-- src/user/benchmark_initialization.F90 | 6 +-- 15 files changed, 123 insertions(+), 103 deletions(-) diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 1000ba0d32..f6791a3b73 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -1889,8 +1889,7 @@ subroutine convective_adjustment(G, GV, h, tv) do j = G%jsc-1,G%jec+1 ; do i = G%isc-1,G%iec+1 ! Compute densities within current water column - call calculate_density( tv%T(i,j,:), tv%S(i,j,:), p_col, & - densities, 1, GV%ke, tv%eqn_of_state ) + call calculate_density( tv%T(i,j,:), tv%S(i,j,:), p_col, densities, tv%eqn_of_state ) ! Repeat restratification until complete do @@ -1909,8 +1908,7 @@ subroutine convective_adjustment(G, GV, h, tv) tv%S(i,j,k) = S1 ; tv%S(i,j,k+1) = S0 h(i,j,k) = h1 ; h(i,j,k+1) = h0 ! Recompute densities at levels k and k+1 - call calculate_density( tv%T(i,j,k), tv%S(i,j,k), p_col(k), & - densities(k), tv%eqn_of_state ) + call calculate_density( tv%T(i,j,k), tv%S(i,j,k), p_col(k), densities(k), tv%eqn_of_state) call calculate_density( tv%T(i,j,k+1), tv%S(i,j,k+1), p_col(k+1), & densities(k+1), tv%eqn_of_state ) stratified = .false. diff --git a/src/ALE/coord_adapt.F90 b/src/ALE/coord_adapt.F90 index 98bbeb7b10..3a083af2db 100644 --- a/src/ALE/coord_adapt.F90 +++ b/src/ALE/coord_adapt.F90 @@ -156,7 +156,7 @@ subroutine build_adapt_column(CS, G, GV, tv, i, j, zInt, tInt, sInt, h, zNext) 0.5 * (tInt(i,j,2:nz) + tInt(i,j-1,2:nz)), & 0.5 * (sInt(i,j,2:nz) + sInt(i,j-1,2:nz)), & 0.5 * (zInt(i,j,2:nz) + zInt(i,j-1,2:nz)) * GV%H_to_Pa, & - alpha, beta, 2, nz - 1, tv%eqn_of_state) + alpha, beta, tv%eqn_of_state, dom=(/2,nz/)) del2sigma(2:nz) = del2sigma(2:nz) + & (alpha(2:nz) * (tInt(i,j-1,2:nz) - tInt(i,j,2:nz)) + & @@ -168,7 +168,7 @@ subroutine build_adapt_column(CS, G, GV, tv, i, j, zInt, tInt, sInt, h, zNext) 0.5 * (tInt(i,j,2:nz) + tInt(i,j+1,2:nz)), & 0.5 * (sInt(i,j,2:nz) + sInt(i,j+1,2:nz)), & 0.5 * (zInt(i,j,2:nz) + zInt(i,j+1,2:nz)) * GV%H_to_Pa, & - alpha, beta, 2, nz - 1, tv%eqn_of_state) + alpha, beta, tv%eqn_of_state, dom=(/2,nz/)) del2sigma(2:nz) = del2sigma(2:nz) + & (alpha(2:nz) * (tInt(i,j+1,2:nz) - tInt(i,j,2:nz)) + & @@ -180,7 +180,7 @@ subroutine build_adapt_column(CS, G, GV, tv, i, j, zInt, tInt, sInt, h, zNext) 0.5 * (tInt(i,j,2:nz) + tInt(i-1,j,2:nz)), & 0.5 * (sInt(i,j,2:nz) + sInt(i-1,j,2:nz)), & 0.5 * (zInt(i,j,2:nz) + zInt(i-1,j,2:nz)) * GV%H_to_Pa, & - alpha, beta, 2, nz - 1, tv%eqn_of_state) + alpha, beta, tv%eqn_of_state, dom=(/2,nz/)) del2sigma(2:nz) = del2sigma(2:nz) + & (alpha(2:nz) * (tInt(i-1,j,2:nz) - tInt(i,j,2:nz)) + & @@ -192,7 +192,7 @@ subroutine build_adapt_column(CS, G, GV, tv, i, j, zInt, tInt, sInt, h, zNext) 0.5 * (tInt(i,j,2:nz) + tInt(i+1,j,2:nz)), & 0.5 * (sInt(i,j,2:nz) + sInt(i+1,j,2:nz)), & 0.5 * (zInt(i,j,2:nz) + zInt(i+1,j,2:nz)) * GV%H_to_Pa, & - alpha, beta, 2, nz - 1, tv%eqn_of_state) + alpha, beta, tv%eqn_of_state, dom=(/2,nz/)) del2sigma(2:nz) = del2sigma(2:nz) + & (alpha(2:nz) * (tInt(i+1,j,2:nz) - tInt(i,j,2:nz)) + & @@ -206,7 +206,7 @@ subroutine build_adapt_column(CS, G, GV, tv, i, j, zInt, tInt, sInt, h, zNext) ! a positive curvature means we're too light relative to adjacent columns, ! so del2sigma needs to be positive too (push the interface deeper) call calculate_density_derivs(tInt(i,j,:), sInt(i,j,:), zInt(i,j,:) * GV%H_to_Pa, & - alpha, beta, 1, nz + 1, tv%eqn_of_state) + alpha, beta, tv%eqn_of_state, dom=(/1,nz/)) !### This should be (/1,nz+1/) - see 25 lines below. do K = 2, nz ! TODO make lower bound here configurable del2sigma(K) = del2sigma(K) * (0.5 * (h(i,j,k-1) + h(i,j,k))) / & diff --git a/src/core/MOM_PressureForce_Montgomery.F90 b/src/core/MOM_PressureForce_Montgomery.F90 index 0d8cf27dad..618199dde1 100644 --- a/src/core/MOM_PressureForce_Montgomery.F90 +++ b/src/core/MOM_PressureForce_Montgomery.F90 @@ -124,12 +124,13 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb real :: alpha_Lay(SZK_(G)) ! The specific volume of each layer [R-1 ~> m3 kg-1]. real :: dalpha_int(SZK_(G)+1) ! The change in specific volume across each ! interface [R-1 ~> m3 kg-1]. - integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb, start, npts + integer, dimension(2) :: EOSdom ! The computational domain for the equation of state + integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb integer :: i, j, k is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke nkmb=GV%nk_rho_varies Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB - start = Isq - (G%isd-1) ; npts = G%iec - Isq + 2 + EOSdom(1) = Isq - (G%isd-1) ; EOSdom(2) = G%iec+1 - (G%isd-1) use_p_atm = .false. if (present(p_atm)) then ; if (associated(p_atm)) use_p_atm = .true. ; endif @@ -228,8 +229,8 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb do k=1,nkmb ; do i=Isq,Ieq+1 tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k) enddo ; enddo - call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), start, npts, & - tv%eqn_of_state, US=US) + call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), & + tv%eqn_of_state, US=US, dom=EOSdom) do k=nkmb+1,nz ; do i=Isq,Ieq+1 if (GV%Rlay(k) < Rho_cv_BL(i)) then tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb) @@ -245,8 +246,8 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb endif !$OMP parallel do default(shared) private(rho_in_situ) do k=1,nz ; do j=Jsq,Jeq+1 - call calculate_density(tv_tmp%T(:,j,k), tv_tmp%S(:,j,k), p_ref, rho_in_situ, start, npts, & - tv%eqn_of_state, US=US) + call calculate_density(tv_tmp%T(:,j,k), tv_tmp%S(:,j,k), p_ref, rho_in_situ, & + tv%eqn_of_state, US=US, dom=EOSdom) do i=Isq,Ieq+1 ; alpha_star(i,j,k) = 1.0 / rho_in_situ(i) ; enddo enddo ; enddo endif ! use_EOS @@ -409,13 +410,14 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce, ! gradient terms are to be split into ! barotropic and baroclinic pieces. type(thermo_var_ptrs) :: tv_tmp! A structure of temporary T & S. - integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb, start, npts + integer, dimension(2) :: EOSdom ! The computational domain for the equation of state + integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb integer :: i, j, k is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke nkmb=GV%nk_rho_varies Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB - start = Isq - (G%isd-1) ; npts = G%iec - Isq + 2 + EOSdom(1) = Isq - (G%isd-1) ; EOSdom(2) = G%iec+1 - (G%isd-1) use_p_atm = .false. if (present(p_atm)) then ; if (associated(p_atm)) use_p_atm = .true. ; endif @@ -484,8 +486,8 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce, do k=1,nkmb ; do i=Isq,Ieq+1 tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k) enddo ; enddo - call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), start, npts, & - tv%eqn_of_state, US=US) + call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), & + tv%eqn_of_state, US=US, dom=EOSdom) do k=nkmb+1,nz ; do i=Isq,Ieq+1 if (GV%Rlay(k) < Rho_cv_BL(i)) then @@ -505,8 +507,8 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce, ! will come down with a fatal error if there is any compressibility. !$OMP parallel do default(shared) do k=1,nz+1 ; do j=Jsq,Jeq+1 - call calculate_density(tv_tmp%T(:,j,k), tv_tmp%S(:,j,k), p_ref, rho_star(:,j,k), start, npts, & - tv%eqn_of_state, US=US) + call calculate_density(tv_tmp%T(:,j,k), tv_tmp%S(:,j,k), p_ref, rho_star(:,j,k), & + tv%eqn_of_state, US=US, dom=EOSdom) do i=Isq,Ieq+1 ; rho_star(i,j,k) = G_Rho0*rho_star(i,j,k) ; enddo enddo ; enddo endif ! use_EOS @@ -632,10 +634,11 @@ subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star) ! an equation of state. real :: z_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [Z ~> m]. - integer :: Isq, Ieq, Jsq, Jeq, nz, i, j, k, start, npts + integer, dimension(2) :: EOSdom ! The computational domain for the equation of state + integer :: Isq, Ieq, Jsq, Jeq, nz, i, j, k Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB ; nz = G%ke - start = Isq - (G%isd-1) ; npts = G%iec - Isq + 2 + EOSdom(1) = Isq - (G%isd-1) ; EOSdom(2) = G%iec+1 - (G%isd-1) Rho0xG = Rho0 * GV%g_Earth G_Rho0 = GV%g_Earth / GV%Rho0 @@ -662,8 +665,8 @@ subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star) Ihtot(i) = GV%H_to_Z / ((e(i,j,1)-e(i,j,nz+1)) + z_neglect) press(i) = -Rho0xG*e(i,j,1) enddo - call calculate_density(tv%T(:,j,1), tv%S(:,j,1), press, rho_in_situ, start, npts, & - tv%eqn_of_state, US=US) + call calculate_density(tv%T(:,j,1), tv%S(:,j,1), press, rho_in_situ, & + tv%eqn_of_state, US=US, dom=EOSdom) do i=Isq,Ieq+1 pbce(i,j,1) = G_Rho0*(GFS_scale * rho_in_situ(i)) * GV%H_to_Z enddo @@ -673,8 +676,8 @@ subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star) T_int(i) = 0.5*(tv%T(i,j,k-1)+tv%T(i,j,k)) S_int(i) = 0.5*(tv%S(i,j,k-1)+tv%S(i,j,k)) enddo - call calculate_density_derivs(T_int, S_int, press, dR_dT, dR_dS, start, npts, & - tv%eqn_of_state, US=US) + call calculate_density_derivs(T_int, S_int, press, dR_dT, dR_dS, & + tv%eqn_of_state, US=US, dom=EOSdom) do i=Isq,Ieq+1 pbce(i,j,k) = pbce(i,j,k-1) + G_Rho0 * & ((e(i,j,K) - e(i,j,nz+1)) * Ihtot(i)) * & @@ -733,10 +736,11 @@ subroutine Set_pbce_nonBouss(p, tv, G, GV, US, GFS_scale, pbce, alpha_star) real :: dp_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [R L2 T-2 ~> Pa]. logical :: use_EOS ! If true, density is calculated from T & S using an equation of state. - integer :: Isq, Ieq, Jsq, Jeq, nz, i, j, k, start, npts + integer, dimension(2) :: EOSdom ! The computational domain for the equation of state + integer :: Isq, Ieq, Jsq, Jeq, nz, i, j, k Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB ; nz = G%ke - start = Isq - (G%isd-1) ; npts = G%iec - Isq + 2 + EOSdom(1) = Isq - (G%isd-1) ; EOSdom(2) = G%iec+1 - (G%isd-1) use_EOS = associated(tv%eqn_of_state) @@ -759,8 +763,8 @@ subroutine Set_pbce_nonBouss(p, tv, G, GV, US, GFS_scale, pbce, alpha_star) else !$OMP parallel do default(shared) private(T_int,S_int,dR_dT,dR_dS,rho_in_situ) do j=Jsq,Jeq+1 - call calculate_density(tv%T(:,j,nz), tv%S(:,j,nz), p(:,j,nz+1), rho_in_situ, start, npts, & - tv%eqn_of_state, US=US) + call calculate_density(tv%T(:,j,nz), tv%S(:,j,nz), p(:,j,nz+1), rho_in_situ, & + tv%eqn_of_state, US=US, dom=EOSdom) do i=Isq,Ieq+1 C_htot(i,j) = dP_dH / ((p(i,j,nz+1)-p(i,j,1)) + dp_neglect) pbce(i,j,nz) = dP_dH / (rho_in_situ(i)) @@ -770,9 +774,9 @@ subroutine Set_pbce_nonBouss(p, tv, G, GV, US, GFS_scale, pbce, alpha_star) T_int(i) = 0.5*(tv%T(i,j,k)+tv%T(i,j,k+1)) S_int(i) = 0.5*(tv%S(i,j,k)+tv%S(i,j,k+1)) enddo - call calculate_density(T_int, S_int, p(:,j,k+1), rho_in_situ, start, npts, tv%eqn_of_state, US=US) - call calculate_density_derivs(T_int, S_int, p(:,j,k+1), dR_dT, dR_dS, start, npts, & - tv%eqn_of_state, US=US) + call calculate_density(T_int, S_int, p(:,j,k+1), rho_in_situ, tv%eqn_of_state, US=US, dom=EOSdom) + call calculate_density_derivs(T_int, S_int, p(:,j,k+1), dR_dT, dR_dS, & + tv%eqn_of_state, US=US, dom=EOSdom) do i=Isq,Ieq+1 pbce(i,j,k) = pbce(i,j,k+1) + ((p(i,j,K+1)-p(i,j,1))*C_htot(i,j)) * & ((dR_dT(i)*(tv%T(i,j,k+1)-tv%T(i,j,k)) + & diff --git a/src/core/MOM_PressureForce_analytic_FV.F90 b/src/core/MOM_PressureForce_analytic_FV.F90 index f0a4485399..b32d81fdbc 100644 --- a/src/core/MOM_PressureForce_analytic_FV.F90 +++ b/src/core/MOM_PressureForce_analytic_FV.F90 @@ -176,13 +176,14 @@ subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p real :: H_to_RL2_T2 ! A factor to convert from thicknesss units (H) to pressure units [R L2 T-2 H-1 ~> Pa H-1]. ! real :: oneatm = 101325.0 ! 1 atm in [Pa] = [kg m-1 s-2] real, parameter :: C1_6 = 1.0/6.0 - integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb, start, npts + integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb + integer, dimension(2) :: EOSdom integer :: i, j, k is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke nkmb=GV%nk_rho_varies Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB - start = Isq - (G%isd-1) ; npts = G%iec - Isq + 2 + EOSdom(1) = Isq - (G%isd-1) ; EOSdom(2) = G%iec+1 - (G%isd-1) if (.not.associated(CS)) call MOM_error(FATAL, & "MOM_PressureForce_AFV_nonBouss: Module must be initialized before it is used.") @@ -228,8 +229,8 @@ subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p do k=1,nkmb ; do i=Isq,Ieq+1 tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k) enddo ; enddo - call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), start, npts, & - tv%eqn_of_state, US=US) + call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), & + tv%eqn_of_state, US=US, dom=EOSdom) do k=nkmb+1,nz ; do i=Isq,Ieq+1 if (GV%Rlay(k) < Rho_cv_BL(i)) then tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb) @@ -334,8 +335,8 @@ subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p if (use_EOS) then !$OMP parallel do default(shared) private(rho_in_situ) do j=Jsq,Jeq+1 - call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p(:,j,1), rho_in_situ, start, npts, & - tv%eqn_of_state, US=US) + call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p(:,j,1), rho_in_situ, & + tv%eqn_of_state, US=US, dom=EOSdom) do i=Isq,Ieq+1 dM(i,j) = (CS%GFS_scale - 1.0) * (p(i,j,1)*(1.0/rho_in_situ(i) - alpha_ref) + za(i,j)) @@ -504,13 +505,14 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_at type(thermo_var_ptrs) :: tv_tmp! A structure of temporary T & S. real, parameter :: C1_6 = 1.0/6.0 - integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb, start, npts + integer, dimension(2) :: EOSdom ! The computational domain for the equation of state + integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb integer :: i, j, k is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke nkmb=GV%nk_rho_varies Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB - start = Isq - (G%isd-1) ; npts = G%iec - Isq + 2 + EOSdom(1) = Isq - (G%isd-1) ; EOSdom(2) = G%iec+1 - (G%isd-1) if (.not.associated(CS)) call MOM_error(FATAL, & "MOM_PressureForce_AFV_Bouss: Module must be initialized before it is used.") @@ -578,8 +580,8 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_at do k=1,nkmb ; do i=Isq,Ieq+1 tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k) enddo ; enddo - call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), start, npts, & - tv%eqn_of_state, US=US) + call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), & + tv%eqn_of_state, US=US, dom=EOSdom) do k=nkmb+1,nz ; do i=Isq,Ieq+1 if (GV%Rlay(k) < Rho_cv_BL(i)) then @@ -601,11 +603,11 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_at !$OMP parallel do default(shared) do j=Jsq,Jeq+1 if (use_p_atm) then - call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p_atm(:,j), rho_in_situ, start, npts, & - tv%eqn_of_state, US=US) + call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p_atm(:,j), rho_in_situ, & + tv%eqn_of_state, US=US, dom=EOSdom) else - call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p0, rho_in_situ, start, npts, & - tv%eqn_of_state, US=US) + call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p0, rho_in_situ, & + tv%eqn_of_state, US=US, dom=EOSdom) endif do i=Isq,Ieq+1 dM(i,j) = (CS%GFS_scale - 1.0) * (G_Rho0 * rho_in_situ(i)) * e(i,j,1) diff --git a/src/core/MOM_PressureForce_blocked_AFV.F90 b/src/core/MOM_PressureForce_blocked_AFV.F90 index 5ac7831479..8e2f3c1405 100644 --- a/src/core/MOM_PressureForce_blocked_AFV.F90 +++ b/src/core/MOM_PressureForce_blocked_AFV.F90 @@ -173,14 +173,15 @@ subroutine PressureForce_blk_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, real :: H_to_RL2_T2 ! A factor to convert from thicknesss units (H) to pressure units [R L2 T-2 H-1 ~> Pa H-1]. ! real :: oneatm = 101325.0 ! 1 atm in [Pa] = [kg m-1 s-2] real, parameter :: C1_6 = 1.0/6.0 - integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb, start, npts + integer, dimension(2) :: EOSdom ! The computational domain for the equation of state + integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb integer :: is_bk, ie_bk, js_bk, je_bk, Isq_bk, Ieq_bk, Jsq_bk, Jeq_bk integer :: i, j, k, n, ib, jb, ioff_bk, joff_bk is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke nkmb=GV%nk_rho_varies Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB - start = Isq - (G%isd-1) ; npts = G%iec - Isq + 2 + EOSdom(1) = Isq - (G%isd-1) ; EOSdom(2) = G%iec+1 - (G%isd-1) if (.not.associated(CS)) call MOM_error(FATAL, & "MOM_PressureForce_AFV_nonBouss: Module must be initialized before it is used.") @@ -224,8 +225,8 @@ subroutine PressureForce_blk_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, do k=1,nkmb ; do i=Isq,Ieq+1 tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k) enddo ; enddo - call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), start, npts, & - tv%eqn_of_state, US=US) + call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), & + tv%eqn_of_state, US=US, dom=EOSdom) do k=nkmb+1,nz ; do i=Isq,Ieq+1 if (GV%Rlay(k) < Rho_cv_BL(i)) then tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb) @@ -302,8 +303,8 @@ subroutine PressureForce_blk_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, if (use_EOS) then !$OMP parallel do default(shared) private(rho_in_situ) do j=Jsq,Jeq+1 - call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p(:,j,1), rho_in_situ, start, npts, & - tv%eqn_of_state, US=US) + call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p(:,j,1), rho_in_situ, & + tv%eqn_of_state, US=US, dom=EOSdom) do i=Isq,Ieq+1 dM(i,j) = (CS%GFS_scale - 1.0) * (p(i,j,1)*(1.0/rho_in_situ(i) - alpha_ref) + za(i,j)) @@ -489,7 +490,8 @@ subroutine PressureForce_blk_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, type(thermo_var_ptrs) :: tv_tmp! A structure of temporary T & S. real, parameter :: C1_6 = 1.0/6.0 - integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb, start, npts + integer, dimension(2) :: EOSdom ! The computational domain for the equation of state + integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb integer :: is_bk, ie_bk, js_bk, je_bk, Isq_bk, Ieq_bk, Jsq_bk, Jeq_bk integer :: ioff_bk, joff_bk integer :: i, j, k, n, ib, jb @@ -497,7 +499,7 @@ subroutine PressureForce_blk_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke nkmb=GV%nk_rho_varies Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB - start = Isq - (G%isd-1) ; npts = G%iec - Isq + 2 + EOSdom(1) = Isq - (G%isd-1) ; EOSdom(2) = G%iec+1 - (G%isd-1) if (.not.associated(CS)) call MOM_error(FATAL, & "MOM_PressureForce_AFV_Bouss: Module must be initialized before it is used.") @@ -565,8 +567,8 @@ subroutine PressureForce_blk_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, do k=1,nkmb ; do i=Isq,Ieq+1 tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k) enddo ; enddo - call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), start, npts, & - tv%eqn_of_state, US=US) + call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), & + tv%eqn_of_state, US=US, dom=EOSdom) do k=nkmb+1,nz ; do i=Isq,Ieq+1 if (GV%Rlay(k) < Rho_cv_BL(i)) then @@ -588,11 +590,11 @@ subroutine PressureForce_blk_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, !$OMP parallel do default(shared) do j=Jsq,Jeq+1 if (use_p_atm) then - call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p_atm(:,j), rho_in_situ, start, npts, & - tv%eqn_of_state, US=US) + call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p_atm(:,j), rho_in_situ, & + tv%eqn_of_state, US=US, dom=EOSdom) else - call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p0, rho_in_situ, start, npts, & - tv%eqn_of_state, US=US) + call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p0, rho_in_situ, & + tv%eqn_of_state, US=US, dom=EOSdom) endif do i=Isq,Ieq+1 dM(i,j) = (CS%GFS_scale - 1.0) * (G_Rho0 * rho_in_situ(i)) * e(i,j,1) diff --git a/src/core/MOM_isopycnal_slopes.F90 b/src/core/MOM_isopycnal_slopes.F90 index 78fdc51077..42d8abe308 100644 --- a/src/core/MOM_isopycnal_slopes.F90 +++ b/src/core/MOM_isopycnal_slopes.F90 @@ -99,6 +99,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & real :: H_to_Z ! A conversion factor from thickness units to the units of e. logical :: present_N2_u, present_N2_v + integer, dimension(2) :: EOSdom_u, EOSdom_v ! Domains for the equation of state calculations at u and v points integer :: is, ie, js, je, nz, IsdB integer :: i, j, k @@ -155,9 +156,11 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & enddo ; enddo enddo - !$OMP parallel do default(none) shared(nz,is,ie,js,je,IsdB,use_EOS,G,GV,US,pres,T,S,tv, & - !$OMP h,h_neglect,e,dz_neglect,Z_to_L,L_to_Z,H_to_Z, & - !$OMP h_neglect2,present_N2_u,G_Rho0,N2_u,slope_x) & + EOSdom_u(1) = is-1 - (G%IsdB-1) ; EOSdom_u(2) = ie - (G%IsdB-1) + + !$OMP parallel do default(none) shared(nz,is,ie,js,je,IsdB,use_EOS,G,GV,US,pres,T,S,tv,h,e, & + !$OMP h_neglect,dz_neglect,Z_to_L,L_to_Z,H_to_Z,h_neglect2, & + !$OMP present_N2_u,G_Rho0,N2_u,slope_x,EOSdom_u) & !$OMP private(drdiA,drdiB,drdkL,drdkR,pres_u,T_u,S_u, & !$OMP drho_dT_u,drho_dS_u,hg2A,hg2B,hg2L,hg2R,haA, & !$OMP haB,haL,haR,dzaL,dzaR,wtA,wtB,wtL,wtR,drdz, & @@ -175,8 +178,8 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & T_u(I) = 0.25*((T(i,j,k) + T(i+1,j,k)) + (T(i,j,k-1) + T(i+1,j,k-1))) 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, US=US) + call calculate_density_derivs(T_u, S_u, pres_u, drho_dT_u, drho_dS_u, & + tv%eqn_of_state, US=US, dom=EOSdom_u) endif do I=is-1,ie @@ -241,10 +244,12 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & enddo ! I enddo ; enddo ! end of j-loop + EOSdom_v(1) = is - (G%isd-1) ; EOSdom_v(2) = ie - (G%isd-1) + ! Calculate the meridional isopycnal slope. !$OMP parallel do default(none) shared(nz,is,ie,js,je,IsdB,use_EOS,G,GV,US,pres,T,S,tv, & !$OMP h,h_neglect,e,dz_neglect,Z_to_L,L_to_Z,H_to_Z, & - !$OMP h_neglect2,present_N2_v,G_Rho0,N2_v,slope_y) & + !$OMP h_neglect2,present_N2_v,G_Rho0,N2_v,slope_y,EOSdom_v) & !$OMP private(drdjA,drdjB,drdkL,drdkR,pres_v,T_v,S_v, & !$OMP drho_dT_v,drho_dS_v,hg2A,hg2B,hg2L,hg2R,haA, & !$OMP haB,haL,haR,dzaL,dzaR,wtA,wtB,wtL,wtR,drdz, & @@ -261,8 +266,8 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & T_v(i) = 0.25*((T(i,j,k) + T(i,j+1,k)) + (T(i,j,k-1) + T(i,j+1,k-1))) 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, US=US) + call calculate_density_derivs(T_v, S_v, pres_v, drho_dT_v, drho_dS_v, tv%eqn_of_state, & + US=US, dom=EOSdom_v) endif do i=is,ie if (use_EOS) then diff --git a/src/diagnostics/MOM_wave_speed.F90 b/src/diagnostics/MOM_wave_speed.F90 index 85dbcdc13b..e1835261aa 100644 --- a/src/diagnostics/MOM_wave_speed.F90 +++ b/src/diagnostics/MOM_wave_speed.F90 @@ -242,8 +242,8 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, & T_int(k) = 0.5*(Tf(k,i)+Tf(k-1,i)) S_int(k) = 0.5*(Sf(k,i)+Sf(k-1,i)) enddo - call calculate_density_derivs(T_int, S_int, pres, drho_dT, drho_dS, 2, & - kf(i)-1, tv%eqn_of_state, US) + call calculate_density_derivs(T_int, S_int, pres, drho_dT, drho_dS, & + tv%eqn_of_state, US, dom=(/2,kf(i)/)) ! Sum the reduced gravities to find out how small a density difference ! is negligibly small. @@ -737,8 +737,8 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, full_halos) T_int(k) = 0.5*(Tf(k,i)+Tf(k-1,i)) S_int(k) = 0.5*(Sf(k,i)+Sf(k-1,i)) enddo - call calculate_density_derivs(T_int, S_int, pres, drho_dT, drho_dS, 2, & - kf(i)-1, tv%eqn_of_state, US) + call calculate_density_derivs(T_int, S_int, pres, drho_dT, drho_dS, & + tv%eqn_of_state, US, dom=(/2,kf(i)/)) ! Sum the reduced gravities to find out how small a density difference ! is negligibly small. diff --git a/src/diagnostics/MOM_wave_structure.F90 b/src/diagnostics/MOM_wave_structure.F90 index ceb6fd6c4f..bc5a06e0ea 100644 --- a/src/diagnostics/MOM_wave_structure.F90 +++ b/src/diagnostics/MOM_wave_structure.F90 @@ -277,8 +277,8 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo T_int(k) = 0.5*(Tf(k,i)+Tf(k-1,i)) S_int(k) = 0.5*(Sf(k,i)+Sf(k-1,i)) enddo - call calculate_density_derivs(T_int, S_int, pres, drho_dT, drho_dS, 2, & - kf(i)-1, tv%eqn_of_state, US) + call calculate_density_derivs(T_int, S_int, pres, drho_dT, drho_dS, & + tv%eqn_of_state, US, dom=(/2,kf(i)/)) ! Sum the reduced gravities to find out how small a density difference ! is negligibly small. diff --git a/src/initialization/MOM_coord_initialization.F90 b/src/initialization/MOM_coord_initialization.F90 index b0155ae603..42a02a70ff 100644 --- a/src/initialization/MOM_coord_initialization.F90 +++ b/src/initialization/MOM_coord_initialization.F90 @@ -291,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, US=US) + call calculate_density(T0, S0, Pref, Rlay, eqn_of_state, US=US, dom=(/1,nz/)) do k=2,nz; g_prime(k) = (GV%g_Earth/(GV%Rho0)) * (Rlay(k) - Rlay(k-1)) ; enddo call callTree_leave(trim(mdl)//'()') @@ -371,7 +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, US=US) + call calculate_density(T0, S0, Pref, Rlay, eqn_of_state, US=US, dom=(/k_light,nz/) ) ! 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) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 2ac8ac47bf..76e87aeed2 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -1610,8 +1610,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, US=US) - call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, nz, eqn_of_state, US=US) + call calculate_density(T0, S0, pres, rho_guess, eqn_of_state, US=US) + call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 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 @@ -1622,8 +1622,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, US=US) - call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, nz, eqn_of_state, US=US) + call calculate_density(T0, S0, pres, rho_guess, eqn_of_state, US=US) + call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, eqn_of_state, US=US) do k=1,nz T0(k) = T0(k) + (GV%Rlay(k) - rho_guess(k)) / drho_dT(k) enddo diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index 62f1dad445..49d1665b16 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -672,6 +672,10 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV real, dimension(SZI_(G), SZJB_(G), SZK_(G)+1) :: diag_sfn_y, diag_sfn_unlim_y ! Diagnostics logical :: present_int_slope_u, present_int_slope_v logical :: present_slope_x, present_slope_y, calc_derivatives + integer, dimension(2) :: EOSdom_u ! The shifted i-computational domain to use for equation of + ! state calculations at u-points. + integer, dimension(2) :: EOSdom_v ! The shifted I-computational domain to use for equation of + ! state calculations at v-points. integer :: is, ie, js, je, nz, IsdB integer :: i, j, k is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke ; IsdB = G%IsdB @@ -747,12 +751,13 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV enddo ; enddo !$OMP end parallel + EOSdom_u(1) = (is-1) - (G%IsdB-1) ; EOSdom_u(2) = ie - (G%IsdB-1) !$OMP parallel do default(none) shared(nz,is,ie,js,je,find_work,use_EOS,G,GV,US,pres,T,S, & !$OMP nk_linear,IsdB,tv,h,h_neglect,e,dz_neglect, & !$OMP I_slope_max2,h_neglect2,present_int_slope_u, & !$OMP int_slope_u,KH_u,uhtot,h_frac,h_avail_rsum, & !$OMP uhD,h_avail,G_scale,Work_u,CS,slope_x,cg1, & -!$OMP diag_sfn_x, diag_sfn_unlim_x,N2_floor, & +!$OMP diag_sfn_x, diag_sfn_unlim_x,N2_floor,EOS_dom_u, & !$OMP present_slope_x,G_rho0,Slope_x_PE,hN2_x_PE) & !$OMP private(drdiA,drdiB,drdkL,drdkR,pres_u,T_u,S_u, & !$OMP drho_dT_u,drho_dS_u,hg2A,hg2B,hg2L,hg2R,haA, & @@ -778,8 +783,8 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV T_u(I) = 0.25*((T(i,j,k) + T(i+1,j,k)) + (T(i,j,k-1) + T(i+1,j,k-1))) 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, US=US) + call calculate_density_derivs(T_u, S_u, pres_u, drho_dT_u, drho_dS_u, & + tv%eqn_of_state, US=US, dom=EOSdom_u) endif do I=is-1,ie @@ -1000,13 +1005,14 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV enddo ! end of j-loop ! Calculate the meridional fluxes and gradients. + EOSdom_v(:) = EOS_domain(G%HI) !$OMP parallel do default(none) shared(nz,is,ie,js,je,find_work,use_EOS,G,GV,US,pres,T,S, & !$OMP nk_linear,IsdB,tv,h,h_neglect,e,dz_neglect, & !$OMP I_slope_max2,h_neglect2,present_int_slope_v, & !$OMP int_slope_v,KH_v,vhtot,h_frac,h_avail_rsum, & !$OMP vhD,h_avail,G_scale,Work_v,CS,slope_y,cg1, & -!$OMP diag_sfn_y, diag_sfn_unlim_y,N2_floor, & -!$OMP present_slope_y,G_rho0,Slope_y_PE,hN2_y_PE) & +!$OMP diag_sfn_y,diag_sfn_unlim_y,N2_floor,EOSdom_v,& +!$OMP present_slope_y,G_rho0,Slope_y_PE,hN2_y_PE) & !$OMP private(drdjA,drdjB,drdkL,drdkR,pres_v,T_v,S_v, & !$OMP drho_dT_v,drho_dS_v,hg2A,hg2B,hg2L,hg2R,haA, & !$OMP haB,haL,haR,dzaL,dzaR,wtA,wtB,wtL,wtR,drdz, & @@ -1030,7 +1036,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV 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, & - tv%eqn_of_state, US, dom=EOS_domain(G%HI)) + tv%eqn_of_state, US, dom=EOSdom_v) endif do i=is,ie if (calc_derivatives) then @@ -1253,6 +1259,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV do j=js,je ; do I=is-1,ie ; uhD(I,j,1) = -uhtot(I,j) ; enddo ; enddo do J=js-1,je ; do i=is,ie ; vhD(i,J,1) = -vhtot(i,J) ; enddo ; enddo else + EOSdom_u(1) = (is-1) - (G%IsdB-1) ; EOSdom_u(2) = ie - (G%IsdB-1) !$OMP parallel do default(shared) private(pres_u,T_u,S_u,drho_dT_u,drho_dS_u,drdiB) do j=js,je if (use_EOS) then @@ -1262,7 +1269,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV S_u(I) = 0.5*(S(i,j,1) + S(i+1,j,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, US=US) + tv%eqn_of_state, US=US, dom=EOSdom_u ) endif do I=is-1,ie uhD(I,j,1) = -uhtot(I,j) @@ -1283,6 +1290,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV enddo enddo + EOSdom_v(:) = EOS_domain(G%HI) !$OMP parallel do default(shared) private(pres_v,T_v,S_v,drho_dT_v,drho_dS_v,drdjB) do J=js-1,je if (use_EOS) then @@ -1292,7 +1300,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV S_v(i) = 0.5*(S(i,j,1) + S(i,j+1,1)) enddo call calculate_density_derivs(T_v, S_v, pres_v, drho_dT_v, drho_dS_v, & - tv%eqn_of_state, US, dom=EOS_domain(G%HI)) + tv%eqn_of_state, US, dom=EOSdom_v) endif do i=is,ie vhD(i,J,1) = -vhtot(i,J) diff --git a/src/parameterizations/vertical/MOM_geothermal.F90 b/src/parameterizations/vertical/MOM_geothermal.F90 index 9b91c6453a..289afd19d2 100644 --- a/src/parameterizations/vertical/MOM_geothermal.F90 +++ b/src/parameterizations/vertical/MOM_geothermal.F90 @@ -198,8 +198,8 @@ subroutine geothermal(h, tv, dt, ea, eb, G, GV, US, CS, halo) iej = is-1 ; do i=ie,is,-1 ; if (do_i(i)) then ; iej = i ; exit ; endif ; enddo 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, US=US) + call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_Ref(:), Rcv_BL(:), & + tv%eqn_of_state, US=US, dom=(/isj-(G%isd-1),iej-(G%isd-1)/)) else Rcv_BL(:) = -1.0 endif @@ -248,8 +248,8 @@ subroutine geothermal(h, tv, dt, ea, eb, G, GV, US, CS, halo) 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, US=US) + call calculate_density_derivs(T2(:), S2(:), p_Ref(:), dRcv_dT_, dRcv_dS_, & + tv%eqn_of_state, US=US, dom=(/1,2/) ) dRcv_dT = 0.5*(dRcv_dT_(1) + dRcv_dT_(2)) endif diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index c628cfd1d3..fadc4874cd 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -273,7 +273,8 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, symmetrize) ! accuracy of a single L(:) Newton iteration logical :: use_L0, do_one_L_iter ! Control flags for L(:) Newton iteration logical :: use_BBL_EOS, do_i(SZIB_(G)) - integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, m, n, K2, nkmb, nkml, start, npts + integer, dimension(2) :: EOSdom ! The computational domain for the equation of state + integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, m, n, K2, nkmb, nkml integer :: itt, maxitt=20 type(ocean_OBC_type), pointer :: OBC => NULL() @@ -292,7 +293,6 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, symmetrize) if (present(symmetrize)) then ; if (symmetrize) then Jsq = js-1 ; Isq = is-1 endif ; endif - start = Isq - (G%isd-1) ; npts = G%iec - Isq + 2 if (CS%debug) then call uvchksum("Start set_viscous_BBL [uv]", u, v, G%HI, haloshift=1, scale=US%L_T_to_m_s) @@ -313,11 +313,12 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, symmetrize) ! if (CS%linear_drag) ustar(:) = cdrag_sqrt_Z*CS%drag_bg_vel if ((nkml>0) .and. .not.use_BBL_EOS) then + EOSdom(1) = Isq - (G%isd-1) ; EOSdom(2) = G%iec+1 - (G%isd-1) do i=Isq,Ieq+1 ; p_ref(i) = tv%P_Ref ; enddo !$OMP parallel do default(shared) do k=1,nkmb ; do j=Jsq,Jeq+1 call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_ref, Rml(:,j,k), tv%eqn_of_state, & - US=US, dom=(/start,start+npts-1/)) + US=US, dom=EOSdom) enddo ; enddo endif diff --git a/src/user/DOME_initialization.F90 b/src/user/DOME_initialization.F90 index 315e56051c..de4726dd1d 100644 --- a/src/user/DOME_initialization.F90 +++ b/src/user/DOME_initialization.F90 @@ -360,12 +360,12 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, US, param_file, tr_Reg) ! 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, US=US) - call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, 1, tv%eqn_of_state, US=US) + call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, tv%eqn_of_state, US=US, dom=(/1,1/)) 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, US=US) - call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, nz, tv%eqn_of_state, US=US) + call calculate_density(T0, S0, pres, rho_guess, tv%eqn_of_state, US=US) + call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 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 diff --git a/src/user/benchmark_initialization.F90 b/src/user/benchmark_initialization.F90 index 766474b364..ff76654b28 100644 --- a/src/user/benchmark_initialization.F90 +++ b/src/user/benchmark_initialization.F90 @@ -153,7 +153,7 @@ subroutine benchmark_initialize_thickness(h, G, GV, US, param_file, eqn_of_state enddo T0(k1) = 29.0 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) + call calculate_density_derivs(T0(k1), S0(k1), pres(k1), drho_dT(k1), drho_dS(k1), eqn_of_state, US=US) ! A first guess of the layers' temperatures. do k=1,nz @@ -267,8 +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, US=US) - call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, nz, eqn_of_state, US=US) + call calculate_density(T0, S0, pres, rho_guess, eqn_of_state, US=US) + call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, eqn_of_state, US=US) do k=1,nz T0(k) = T0(k) + (GV%Rlay(k) - rho_guess(k)) / drho_dT(k) enddo