From f3725a74f27f59860ee68bc02fdc60fbc6a8cc5f Mon Sep 17 00:00:00 2001 From: Zhi Liang Date: Mon, 28 Jul 2014 16:45:51 -0400 Subject: [PATCH] Implement Openmp in MOM_PressureForce. Block data structure is used in core/MOM_PressureForce_analytic_FV.F90 --- src/core/MOM_PressureForce_Montgomery.F90 | 207 +++++---- src/core/MOM_PressureForce_analytic_FV.F90 | 478 +++++++++++++-------- src/equation_of_state/MOM_EOS.F90 | 28 +- src/framework/MOM_domains.F90 | 2 + 4 files changed, 448 insertions(+), 267 deletions(-) diff --git a/src/core/MOM_PressureForce_Montgomery.F90 b/src/core/MOM_PressureForce_Montgomery.F90 index 528159f86c..37ca5f866d 100644 --- a/src/core/MOM_PressureForce_Montgomery.F90 +++ b/src/core/MOM_PressureForce_Montgomery.F90 @@ -126,7 +126,8 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta) ! with any tidal contributions or compressibility compensation. real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: & M, & ! The Montgomery potential, M = (p/rho + gz) , in m2 s-2. - alpha_star ! Compression adjusted specific volume, in m3 kg-1. + alpha_star, & ! Compression adjusted specific volume, in m3 kg-1. + dz_geo ! The change in geopotential across a layer, in m2 s-2. real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: p ! Interface pressure in Pa. ! p may be adjusted (with a nonlinear equation of state) so that ! its derivative compensates for the adiabatic compressibility @@ -137,6 +138,9 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta) S_tmp ! Temporary array of salinities where layers that are lighter ! than the mixed layer have the mixed layer's properties, in psu. + real, dimension(SZI_(G)) :: Rho_cv_BL ! The coordinate potential density in the + ! deepest variable density near-surface layer, in kg m-3. + real, dimension(SZI_(G),SZJ_(G)) :: & dM, & ! A barotropic correction to the Montgomery potentials to ! enable the use of a reduced gravity form of the equations, @@ -146,9 +150,6 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta) ! astronomical sources and self-attraction and loading, in m. geopot_bot, & ! Bottom geopotential relative to time-mean sea level, ! including any tidal contributions, in units of m2 s-2. - dz_geo, & ! The change in geopotential across a layer, in m2 s-2. - Rho_cv_BL, & ! The coordinate potential density in the deepest variable - ! density near-surface layer, in kg m-3. SSH ! Sea surface height anomalies, in m. real :: p_ref(SZI_(G)) ! The pressure used to calculate the coordinate ! density, in Pa (usually 2e7 Pa = 2000 dbar). @@ -192,52 +193,76 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta) I_gEarth = 1.0 / G%g_Earth dp_neglect = G%H_to_Pa * G%H_subroundoff +!$OMP parallel default(none) shared(nz,alpha_Lay,G,dalpha_int) +!$OMP do do k=1,nz ; alpha_Lay(k) = 1.0 / G%Rlay(k) ; enddo +!$OMP do do k=2,nz ; dalpha_int(K) = alpha_Lay(k-1) - alpha_Lay(k) ; enddo +!$OMP end parallel +!$OMP parallel default(none) shared(Isq,Ieq,Jsq,Jeq,nz,p,p_atm,G,h,use_p_atm) if (use_p_atm) then +!$OMP do do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 ; p(i,j,1) = p_atm(i,j) ; enddo ; enddo else +!$OMP do do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 ; p(i,j,1) = 0.0 ; enddo ; enddo endif - do k=1,nz ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 +!$OMP do + do j=Jsq,Jeq+1 ; do k=1,nz ; do i=Isq,Ieq+1 p(i,j,K+1) = p(i,j,K) + G%H_to_Pa * h(i,j,k) enddo ; enddo ; enddo +!$OMP end parallel if (present(eta)) then Pa_to_H = 1.0 / G%H_to_Pa - if (use_p_atm) then ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - eta(i,j) = (p(i,j,nz+1) - p_atm(i,j))*Pa_to_H ! eta has the same units as h. - enddo ; enddo ; else ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - eta(i,j) = p(i,j,nz+1)*Pa_to_H ! eta has the same units as h. - enddo ; enddo ; endif + if (use_p_atm) then +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,eta,p,p_atm,Pa_to_H) + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + eta(i,j) = (p(i,j,nz+1) - p_atm(i,j))*Pa_to_H ! eta has the same units as h. + enddo ; enddo + else +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,eta,p,Pa_to_H) + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + eta(i,j) = p(i,j,nz+1)*Pa_to_H ! eta has the same units as h. + enddo ; enddo + endif endif if (CS%tides) then ! Determine the sea surface height anomalies, to enable the calculation ! of self-attraction and loading. +!$OMP parallel default(none) shared(Isq,Ieq,Jsq,Jeq,nz,SSH,G,use_EOS,tv,p,dz_geo, & +!$OMP I_gEarth,h,alpha_Lay) +!$OMP do do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 SSH(i,j) = -G%bathyT(i,j) enddo ; enddo if (use_EOS) then +!$OMP do do k=1,nz call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p(:,:,k), p(:,:,k+1), & - 0.0, G, tv%eqn_of_state, dz_geo, halo_size=1) - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - SSH(i,j) = SSH(i,j) + I_gEarth * dz_geo(i,j) - enddo ; enddo + 0.0, G, tv%eqn_of_state, dz_geo(:,:,k), halo_size=1) enddo +!$OMP do + do j=Jsq,Jeq+1 ; do k=1,nz; do i=Isq,Ieq+1 + SSH(i,j) = SSH(i,j) + I_gEarth * dz_geo(i,j,k) + enddo ; enddo ; enddo else - do k=1,nz ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 +!$OMP do + do j=Jsq,Jeq+1 ; do k=1,nz ; do i=Isq,Ieq+1 SSH(i,j) = SSH(i,j) + G%H_to_kg_m2*h(i,j,k)*alpha_Lay(k) enddo ; enddo ; enddo endif +!$OMP end parallel call calc_tidal_forcing(CS%Time, SSH, e_tidal, G, CS%tides_CSp) +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,geopot_bot,G,e_tidal) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 geopot_bot(i,j) = -G%g_Earth*(e_tidal(i,j) + G%bathyT(i,j)) enddo ; enddo else +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,geopot_bot,G) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 geopot_bot(i,j) = -G%g_Earth*G%bathyT(i,j) enddo ; enddo @@ -253,26 +278,30 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta) if (nkmb>0) then tv_tmp%T => T_tmp ; tv_tmp%S => S_tmp tv_tmp%eqn_of_state => tv%eqn_of_state - do k=1,nkmb ; do j=Jsq,Jeq+1 ; 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 ; enddo do i=Isq,Ieq+1 ; p_ref(i) = tv%P_Ref ; enddo +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,nkmb,tv_tmp,tv,p_ref,G) & +!$OMP private(Rho_cv_BL) do j=Jsq,Jeq+1 + 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(:,j), Isq, Ieq-Isq+2, tv%eqn_of_state) + Rho_cv_BL(:), Isq, Ieq-Isq+2, tv%eqn_of_state) + do k=nkmb+1,nz ; do i=Isq,Ieq+1 + if (G%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) + else + tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k) + endif + enddo ; enddo enddo - do k=nkmb+1,nz ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - if (G%Rlay(k) < Rho_cv_BL(i,j)) then - tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb) - else - tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k) - endif - enddo ; enddo ; enddo else tv_tmp%T => tv%T ; tv_tmp%S => tv%S tv_tmp%eqn_of_state => tv%eqn_of_state + do i=Isq,Ieq+1 ; p_ref(i) = 0 ; enddo endif - +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,tv_tmp,p_ref,tv,alpha_star) & +!$OMP 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,Isq,Ieq-Isq+2,tv%eqn_of_state) @@ -281,26 +310,35 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta) endif ! use_EOS if (use_EOS) then - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - M(i,j,nz) = geopot_bot(i,j) + p(i,j,nz+1) * alpha_star(i,j,nz) - enddo ; enddo - do k=nz-1,1,-1 ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - M(i,j,k) = M(i,j,k+1) + p(i,j,K+1) * (alpha_star(i,j,k) - alpha_star(i,j,k+1)) - enddo ; enddo ; enddo +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,M,geopot_bot,p,alpha_star) + do j=Jsq,Jeq+1 + do i=Isq,Ieq+1 + M(i,j,nz) = geopot_bot(i,j) + p(i,j,nz+1) * alpha_star(i,j,nz) + enddo + do k=nz-1,1,-1 ; do i=Isq,Ieq+1 + M(i,j,k) = M(i,j,k+1) + p(i,j,K+1) * (alpha_star(i,j,k) - alpha_star(i,j,k+1)) + enddo ; enddo + enddo else ! not use_EOS - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - M(i,j,nz) = geopot_bot(i,j) + p(i,j,nz+1) * alpha_Lay(nz) - enddo ; enddo - do k=nz-1,1,-1 ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - M(i,j,k) = M(i,j,k+1) + p(i,j,K+1) * dalpha_int(K+1) - enddo ; enddo ; enddo +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,M,geopot_bot,p,& +!$OMP alpha_Lay,dalpha_int) + do j=Jsq,Jeq+1 + do i=Isq,Ieq+1 + M(i,j,nz) = geopot_bot(i,j) + p(i,j,nz+1) * alpha_Lay(nz) + enddo + do k=nz-1,1,-1 ; do i=Isq,Ieq+1 + M(i,j,k) = M(i,j,k+1) + p(i,j,K+1) * dalpha_int(K+1) + enddo ; enddo + enddo endif ! use_EOS if (CS%GFS_scale < 1.0) then ! Adjust the Montgomery potential to make this a reduced gravity model. +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,dM,CS,M) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 dM(i,j) = (CS%GFS_scale - 1.0) * M(i,j,1) enddo ; enddo +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,dM,M) do k=1,nz ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 M(i,j,k) = M(i,j,k) + dM(i,j) enddo ; enddo ; enddo @@ -331,6 +369,9 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta) ! Calculate the pressure force. On a Cartesian grid, ! PFu = - dM/dx and PFv = - dM/dy. if (use_EOS) then +!$OMP parallel do default(none) shared(is,ie,js,je,Isq,Ieq,Jsq,Jeq,nz,p,dp_neglect, & +!$OMP alpha_star,G,PFu,PFv,M,CS) & +!$OMP private(dp_star,PFu_bc,PFv_bc) do k=1,nz do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 dp_star(i,j) = (p(i,j,K+1) - p(i,j,K)) + dp_neglect @@ -352,6 +393,7 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta) enddo ; enddo enddo ! k-loop else ! .not. use_EOS +!$OMP parallel do default(none) shared(is,ie,js,je,Isq,Ieq,Jsq,Jeq,nz,PFu,PFv,M,G) do k=1,nz do j=js,je ; do I=Isq,Ieq PFu(I,j,k) = -(M(i+1,j,k) - M(i,j,k)) * G%IdxCu(I,j) @@ -417,7 +459,7 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta) S_tmp ! Temporary array of salinities where layers that are lighter ! than the mixed layer have the mixed layer's properties, in psu. - real :: Rho_cv_BL(SZI_(G),SZJ_(G)) ! The coordinate potential density in + real :: Rho_cv_BL(SZI_(G)) ! The coordinate potential density in ! the deepest variable density near-surface layer, in kg m-3. real :: h_star(SZI_(G),SZJ_(G)) ! Layer thickness after compensation ! for compressibility, in m. @@ -511,17 +553,17 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta) tv_tmp%eqn_of_state => tv%eqn_of_state do i=Isq,Ieq+1 ; p_ref(i) = tv%P_Ref ; enddo -!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,nkmb,tv_tmp,tv,p_ref, & -!$OMP Rho_cv_BL,G) +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,nkmb,tv_tmp,tv,p_ref,G) & +!$OMP private(Rho_cv_BL) do j=Jsq,Jeq+1 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(:,j), Isq, Ieq-Isq+2, tv%eqn_of_state) + Rho_cv_BL(:), Isq, Ieq-Isq+2, tv%eqn_of_state) do k=nkmb+1,nz ; do i=Isq,Ieq+1 - if (G%Rlay(k) < Rho_cv_BL(i,j)) then + if (G%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) else tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k) @@ -807,15 +849,22 @@ subroutine Set_pbce_nonBouss(p, tv, G, g_Earth, GFS_scale, pbce, alpha_star) if (use_EOS) then if (present(alpha_star)) then - do j=Jsq,Jeq+1 ; 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 * alpha_star(i,j,nz) - enddo ; enddo - do k=nz-1,1,-1 ; do j=Jsq,Jeq+1 ; 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)) * & - (alpha_star(i,j,k) - alpha_star(i,j,k+1)) - enddo ; enddo ; enddo +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,C_htot,dP_dH,p,dp_neglect, & +!$OMP pbce,alpha_star) + do j=Jsq,Jeq+1 + 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 * alpha_star(i,j,nz) + enddo + do k=nz-1,1,-1 ; 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)) * & + (alpha_star(i,j,k) - alpha_star(i,j,k+1)) + enddo ; enddo + enddo else +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,tv,p,C_htot, & +!$OMP dP_dH,dp_neglect,pbce) & +!$OMP 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, Isq, Ieq-Isq+2, tv%eqn_of_state) @@ -823,42 +872,50 @@ subroutine Set_pbce_nonBouss(p, tv, G, g_Earth, GFS_scale, pbce, alpha_star) 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) enddo + do k=nz-1,1,-1 + do i=Isq,Ieq+1 + 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, & + Isq, Ieq-Isq+2, tv%eqn_of_state) + call calculate_density_derivs(T_int, S_int, p(:,j,k+1), dR_dT, dR_dS, & + Isq, Ieq-Isq+2, tv%eqn_of_state) + 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)) + & + dR_dS(i)*(tv%S(i,j,k+1)-tv%S(i,j,k))) / rho_in_situ(i)**2) + enddo + enddo enddo - do k=nz-1,1,-1 ; do j=Jsq,Jeq+1 - do i=Isq,Ieq+1 - 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, & - Isq, Ieq-Isq+2, tv%eqn_of_state) - call calculate_density_derivs(T_int, S_int, p(:,j,k+1), dR_dT, dR_dS, & - Isq, Ieq-Isq+2, tv%eqn_of_state) - 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)) + & - dR_dS(i)*(tv%S(i,j,k+1)-tv%S(i,j,k))) / rho_in_situ(i)**2) - enddo - enddo ; enddo endif else ! not use_EOS - do j=Jsq,Jeq+1 ; 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 * alpha_Lay(nz) - enddo ; enddo - do k=nz-1,1,-1 ; do j=Jsq,Jeq+1 ; 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)) * & - dalpha_int(K+1) - enddo ; enddo ; enddo +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,C_htot,dP_dH,p,dp_neglect, & +!$OMP pbce,alpha_Lay,dalpha_int) + do j=Jsq,Jeq+1 + 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 * alpha_Lay(nz) + enddo + do k=nz-1,1,-1 ; 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)) * & + dalpha_int(K+1) + enddo ; enddo + enddo endif ! use_EOS if (GFS_scale < 1.0) then ! Adjust the Montgomery potential to make this a reduced gravity model. +!$OMP parallel default(none) shared(Isq,Ieq,Jsq,Jeq,dpbce,GFS_scale,pbce,nz) +!$OMP do do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 dpbce(i,j) = (GFS_scale - 1.0) * pbce(i,j,1) enddo ; enddo +!$OMP do do k=1,nz ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 pbce(i,j,k) = pbce(i,j,k) + dpbce(i,j) enddo ; enddo ; enddo +!$OMP end parallel endif end subroutine Set_pbce_nonBouss diff --git a/src/core/MOM_PressureForce_analytic_FV.F90 b/src/core/MOM_PressureForce_analytic_FV.F90 index 31dd13e109..d640ab42d6 100644 --- a/src/core/MOM_PressureForce_analytic_FV.F90 +++ b/src/core/MOM_PressureForce_analytic_FV.F90 @@ -57,6 +57,7 @@ module MOM_PressureForce_AFV use MOM_diag_mediator, only : safe_alloc_ptr, diag_ctrl, time_type use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe use MOM_file_parser, only : get_param, log_param, log_version, param_file_type +use MOM_domains, only : compute_extent use MOM_grid, only : ocean_grid_type use MOM_PressureForce_Mont, only : set_pbce_Bouss, set_pbce_nonBouss use MOM_tidal_forcing, only : calc_tidal_forcing, tidal_forcing_CS @@ -91,6 +92,40 @@ module MOM_PressureForce_AFV type(tidal_forcing_CS), pointer :: tides_CSp => NULL() end type PressureForce_AFV_CS +type :: block_type + integer :: isc, iec, jsc, jec + integer :: isd, ied, jsd, jed + integer :: IscB, IecB, JscB, JecB + !--- The following are used in PressureForce_AFV_nonBouss + real, dimension(:,:), allocatable :: & + dp, & ! The (positive) change in pressure across a layer, in Pa. + za, & ! The geopotential anomaly (i.e. g*e + alpha_0*pressure) at the + ! interface atop a layer, in m2 s-2. + intx_za, & ! The zonal integral of the geopotential anomaly along the + ! interface below a layer, divided by the grid spacing, m2 s-2. + inty_za ! The meridional integral of the geopotential anomaly along the + ! interface below a layer, divided by the grid spacing, m2 s-2. + !--- The following are used in PressureForce_AFV_Bouss + real, dimension(:,:), allocatable :: & + pa, & ! The pressure anomaly (i.e. pressure + g*RHO_0*e) at the + ! the interface atop a layer, in Pa. + intx_pa, & ! The zonal integral of the pressure anomaly along the interface + ! atop a layer, divided by the grid spacing, in Pa. + inty_pa, & ! The meridional integral of the pressure anomaly along the + ! interface atop a layer, divided by the grid spacing, in Pa. + dpa, & ! The change in pressure anomaly between the top and bottom + ! of a layer, in Pa. + intx_dpa, & ! The change in intx_pa through a layer, in Pa. + inty_dpa, & ! The change in inty_pa through a layer, in Pa. + intz_dpa, & ! The vertical integral in depth of the pressure anomaly less + ! the pressure anomaly at the top of the layer, in m Pa. + dz ! The change in geopotential thickness through a layer, m2 s-2. +end type block_type + +type(ocean_grid_type), allocatable :: G_BLK(:) +type(block_type), allocatable :: block(:) +integer :: nblocks = 1 + contains subroutine PressureForce_AFV(h, tv, PFu, PFv, G, CS, ALE_CSp, p_atm, pbce, eta) @@ -174,18 +209,12 @@ subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta) ! astronomical sources and self-attraction and loading, in m. dM, & ! The barotropic adjustment to the Montgomery potential to ! account for a reduced gravity model, in m2 s-2. - Rho_cv_BL,& ! The coordinate potential density in the deepest variable - ! density near-surface layer, in kg m-3. za ! The geopotential anomaly (i.e. g*e + alpha_0*pressure) at the ! interface atop a layer, in m2 s-2. - real, dimension(SZIB_(G),SZJ_(G)) :: & - intx_za ! The zonal integral of the geopotential anomaly along the - ! interface below a layer, divided by the grid spacing, m2 s-2. + real, dimension(SZI_(G)) :: Rho_cv_BL ! The coordinate potential density in the deepest variable + ! density near-surface layer, in kg m-3. real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: & intx_dza ! The change in intx_za through a layer, in m2 s-2. - real, dimension(SZI_(G),SZJB_(G)) :: & - inty_za ! The meridional integral of the geopotential anomaly along the - ! interface below a layer, divided by the grid spacing, m2 s-2. real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: & inty_dza ! The change in inty_za through a layer, in m2 s-2. real :: p_ref(SZI_(G)) ! The pressure used to calculate the coordinate @@ -208,7 +237,8 @@ subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta) real :: I_gEarth real, parameter :: C1_6 = 1.0/6.0 integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb - integer :: i, j, k + integer :: i, j, k, n + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke nkmb=G%nk_rho_varies Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB @@ -223,18 +253,23 @@ subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta) dp_neglect = G%H_to_Pa * G%H_subroundoff alpha_ref = 1.0/CS%Rho0 +!$OMP parallel default(none) shared(Isq,Ieq,Jsq,Jeq,nz,use_p_atm,p,p_atm,G,h) if (use_p_atm) then +!$OMP do do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 p(i,j,1) = p_atm(i,j) enddo ; enddo else +!$OMP do do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 p(i,j,1) = 0.0 ! or oneatm enddo ; enddo endif - do k=2,nz+1 ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 +!$OMP do + do j=Jsq,Jeq+1 ; do k=2,nz+1 ; do i=Isq,Ieq+1 p(i,j,K) = p(i,j,K-1) + G%H_to_Pa * h(i,j,k-1) enddo ; enddo ; enddo +!$OMP end parallel I_gEarth = 1.0 / G%g_Earth @@ -246,27 +281,32 @@ subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta) if (nkmb>0) then tv_tmp%T => T_tmp ; tv_tmp%S => S_tmp tv_tmp%eqn_of_state => tv%eqn_of_state - do k=1,nkmb ; do j=Jsq,Jeq+1 ; 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 ; enddo do i=Isq,Ieq+1 ; p_ref(i) = tv%P_Ref ; enddo +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,nkmb,tv_tmp,tv,p_ref,G) & +!$OMP private(Rho_cv_BL) do j=Jsq,Jeq+1 + 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(:,j), Isq, Ieq-Isq+2, tv%eqn_of_state) + Rho_cv_BL(:), Isq, Ieq-Isq+2, tv%eqn_of_state) + do k=nkmb+1,nz ; do i=Isq,Ieq+1 + if (G%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) + else + tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k) + endif + enddo ; enddo enddo - do k=nkmb+1,nz ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - if (G%Rlay(k) < Rho_cv_BL(i,j)) then - tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb) - else - tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k) - endif - enddo ; enddo ; enddo else tv_tmp%T => tv%T ; tv_tmp%S => tv%S tv_tmp%eqn_of_state => tv%eqn_of_state endif endif +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,is,ie,js,je,tv_tmp,alpha_ref, & +!$OMP p,h,G,tv,dza,intp_dza,intx_dza,inty_dza,use_EOS) & +!$OMP private(alpha_anom,dp) do k=1,nz ! Calculate 4 integrals through the layer that are required in the ! subsequent calculation. @@ -299,19 +339,24 @@ subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta) ! inty_dza to be 3-D arrays. ! Sum vertically to determine the surface geopotential anomaly. - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - za(i,j) = alpha_ref*p(i,j,nz+1) - G%g_Earth*G%bathyT(i,j) - enddo ; enddo - do k=nz,1,-1 ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,za,alpha_ref,p,G,dza) + do j=Jsq,Jeq+1 + do i=Isq,Ieq+1 + za(i,j) = alpha_ref*p(i,j,nz+1) - G%g_Earth*G%bathyT(i,j) + enddo + do k=nz,1,-1 ; do i=Isq,Ieq+1 za(i,j) = za(i,j) + dza(i,j,k) - enddo ; enddo ; enddo + enddo ; enddo + enddo if (CS%tides) then ! Find and add the tidal geopotential anomaly. +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,SSH,za,alpha_ref,p,I_gEarth) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 SSH(i,j) = (za(i,j) - alpha_ref*p(i,j,1)) * I_gEarth enddo ; enddo call calc_tidal_forcing(CS%Time, SSH, e_tidal, G, CS%tides_CSp) +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,za,G,e_tidal) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 za(i,j) = za(i,j) - G%g_Earth*e_tidal(i,j) enddo ; enddo @@ -320,6 +365,8 @@ subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta) if (CS%GFS_scale < 1.0) then ! Adjust the Montgomery potential to make this a reduced gravity model. if (use_EOS) then +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,tv_tmp,p,tv,dM,CS,alpha_ref,za) & +!$OMP 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, Isq, Ieq-Isq+2, tv%eqn_of_state) @@ -330,6 +377,7 @@ subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta) enddo enddo else +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,dM,CS,p,G,alpha_ref,za) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 dM(i,j) = (CS%GFS_scale - 1.0) * & (p(i,j,1)*(1.0/G%Rlay(1) - alpha_ref) + za(i,j)) @@ -344,45 +392,57 @@ subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta) ! linearly between the values at thickness points, but the bottom ! geopotentials will not now be linear at the sub-grid-scale. Doing this ! ensures no motion with flat isopycnals, even with a nonlinear equation of state. - do j=js,je ; do I=Isq,Ieq - intx_za(I,j) = 0.5*(za(i,j) + za(i+1,j)) - enddo ; enddo - do J=Jsq,Jeq ; do i=is,ie - inty_za(i,J) = 0.5*(za(i,j) + za(i,j+1)) - enddo ; enddo - do k=1,nz - ! These expressions for the acceleration have been carefully checked in - ! a set of idealized cases, and should be bug-free. +!$OMP parallel do default(none) shared(nblocks,block,nz,za,G,dza,intx_dza,h,PFu, & +!$OMP intp_dza,p,dp_neglect,inty_dza,PFv,CS,dM) & +!$OMP private(is,ie,js,je,Isq,Ieq,Jsq,Jeq) + do n = 1, nblocks + is=block(n)%isc; ie=block(n)%iec; js=block(n)%jsc; je=block(n)%jec + Isq=block(n)%IscB; Ieq=block(n)%IecB; Jsq=block(n)%JscB; Jeq=block(n)%JecB do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - dp(i,j) = G%H_to_Pa*h(i,j,k) - za(i,j) = za(i,j) - dza(i,j,k) + block(n)%za(i,j) = za(i,j) enddo ; enddo do j=js,je ; do I=Isq,Ieq - intx_za(I,j) = intx_za(I,j) - intx_dza(I,j,k) - PFu(I,j,k) = (((za(i,j)*dp(i,j) + intp_dza(i,j,k)) - & - (za(i+1,j)*dp(i+1,j) + intp_dza(i+1,j,k))) + & - ((dp(i+1,j) - dp(i,j)) * intx_za(I,j) - & - (p(i+1,j,K) - p(i,j,K)) * intx_dza(I,j,k))) * & - (2.0*G%IdxCu(I,j) / ((dp(i,j) + dp(i+1,j)) + dp_neglect)) + block(n)%intx_za(I,j) = 0.5*(za(i,j) + za(i+1,j)) enddo ; enddo do J=Jsq,Jeq ; do i=is,ie - inty_za(i,J) = inty_za(i,J) - inty_dza(i,J,k) - PFv(i,J,k) = (((za(i,j)*dp(i,j) + intp_dza(i,j,k)) - & - (za(i,j+1)*dp(i,j+1) + intp_dza(i,j+1,k))) + & - ((dp(i,j+1) - dp(i,j)) * inty_za(i,J) - & - (p(i,j+1,K) - p(i,j,K)) * inty_dza(i,J,k))) * & - (2.0*G%IdyCv(i,J) / ((dp(i,j) + dp(i,j+1)) + dp_neglect)) + block(n)%inty_za(i,J) = 0.5*(za(i,j) + za(i,j+1)) enddo ; enddo - - if (CS%GFS_scale < 1.0) then - ! Adjust the Montgomery potential to make this a reduced gravity model. + do k=1,nz + ! These expressions for the acceleration have been carefully checked in + ! a set of idealized cases, and should be bug-free. + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + block(n)%dp(i,j) = G%H_to_Pa*h(i,j,k) + block(n)%za(i,j) = block(n)%za(i,j) - dza(i,j,k) + enddo ; enddo do j=js,je ; do I=Isq,Ieq - PFu(I,j,k) = PFu(I,j,k) - (dM(i+1,j) - dM(i,j)) * G%IdxCu(I,j) + block(n)%intx_za(I,j) = block(n)%intx_za(I,j) - intx_dza(I,j,k) + PFu(I,j,k) = (((block(n)%za(i,j)*block(n)%dp(i,j) + intp_dza(i,j,k)) - & + (block(n)%za(i+1,j)*block(n)%dp(i+1,j) + intp_dza(i+1,j,k))) + & + ((block(n)%dp(i+1,j) - block(n)%dp(i,j)) * block(n)%intx_za(I,j) - & + (p(i+1,j,K) - p(i,j,K)) * intx_dza(I,j,k))) * & + (2.0*G%IdxCu(I,j) / ((block(n)%dp(i,j) + block(n)%dp(i+1,j)) + & + dp_neglect)) enddo ; enddo do J=Jsq,Jeq ; do i=is,ie - PFv(i,J,k) = PFv(i,J,k) - (dM(i,j+1) - dM(i,j)) * G%IdyCv(i,J) + block(n)%inty_za(i,J) = block(n)%inty_za(i,J) - inty_dza(i,J,k) + PFv(i,J,k) = (((block(n)%za(i,j)*block(n)%dp(i,j) + intp_dza(i,j,k)) - & + (block(n)%za(i,j+1)*block(n)%dp(i,j+1) + intp_dza(i,j+1,k))) + & + ((block(n)%dp(i,j+1) - block(n)%dp(i,j)) * block(n)%inty_za(i,J) - & + (p(i,j+1,K) - p(i,j,K)) * inty_dza(i,J,k))) * & + (2.0*G%IdyCv(i,J) / ((block(n)%dp(i,j) + block(n)%dp(i,j+1)) + & + dp_neglect)) enddo ; enddo - endif + + if (CS%GFS_scale < 1.0) then + ! Adjust the Montgomery potential to make this a reduced gravity model. + do j=js,je ; do I=Isq,Ieq + PFu(I,j,k) = PFu(I,j,k) - (dM(i+1,j) - dM(i,j)) * G%IdxCu(I,j) + enddo ; enddo + do J=Jsq,Jeq ; do i=is,ie + PFv(i,J,k) = PFv(i,J,k) - (dM(i,j+1) - dM(i,j)) * G%IdyCv(i,J) + enddo ; enddo + endif + enddo enddo if (present(pbce)) then @@ -390,12 +450,19 @@ subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta) endif if (present(eta)) then + Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB Pa_to_H = 1.0 / G%H_to_Pa - if (use_p_atm) then ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - eta(i,j) = (p(i,j,nz+1) - p_atm(i,j))*Pa_to_H ! eta has the same units as h. - enddo ; enddo ; else ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - eta(i,j) = p(i,j,nz+1)*Pa_to_H ! eta has the same units as h. - enddo ; enddo ; endif + if (use_p_atm) then +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,eta,p,p_atm,Pa_to_H) + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + eta(i,j) = (p(i,j,nz+1) - p_atm(i,j))*Pa_to_H ! eta has the same units as h. + enddo ; enddo + else +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,eta,p,Pa_to_H) + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + eta(i,j) = p(i,j,nz+1)*Pa_to_H ! eta has the same units as h. + enddo ; enddo + endif endif if (CS%id_e_tidal>0) call post_data(CS%id_e_tidal, e_tidal, CS%diag) @@ -442,33 +509,14 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, CS, ALE_CSp, p_atm, pbce, real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: e ! Interface height in m. real, dimension(SZI_(G),SZJ_(G)) :: & - dz, & ! The change in geopotential thickness through a layer, m2 s-2. e_tidal, & ! The bottom geopotential anomaly due to tidal forces from ! astronomical sources and self-attraction and loading, in m. - dM, & ! The barotropic adjustment to the Montgomery potential to + dM ! The barotropic adjustment to the Montgomery potential to ! account for a reduced gravity model, in m2 s-2. - pa ! The pressure anomaly (i.e. pressure + g*RHO_0*e) at the - ! the interface atop a layer, in Pa. real, dimension(SZI_(G)) :: & Rho_cv_BL ! The coordinate potential density in the deepest variable ! density near-surface layer, in kg m-3. - real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: & - dpa, & ! The change in pressure anomaly between the top and bottom - ! of a layer, in Pa. - intz_dpa, & ! The vertical integral in depth of the pressure anomaly less - ! the pressure anomaly at the top of the layer, in m Pa. - pa_h_intz_dpa ! pa*h+intz_dpa - real, dimension(SZIB_(G),SZJ_(G)) :: & - intx_pa ! The zonal integral of the pressure anomaly along the interface - ! atop a layer, divided by the grid spacing, in Pa. - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: & - intx_dpa ! The change in intx_pa through a layer, in Pa. - real, dimension(SZI_(G),SZJB_(G)) :: & - inty_pa ! The meridional integral of the pressure anomaly along the - ! interface atop a layer, divided by the grid spacing, in Pa. - real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: & - inty_dpa ! The change in inty_pa through a layer, in Pa. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), target :: & T_tmp, & ! Temporary array of temperatures where layers that are lighter ! than the mixed layer have the mixed layer's properties, in C. @@ -493,7 +541,7 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, CS, ALE_CSp, p_atm, pbce, 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 - integer :: i, j, k + integer :: i, j, k, n, isd, ied, jsd, jed integer :: PRScheme is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke @@ -586,8 +634,8 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, CS, ALE_CSp, p_atm, pbce, endif !$OMP parallel default(none) shared(Jsq,Jeq,Isq,Ieq,tv_tmp,p_atm,rho_in_situ,tv, & -!$OMP p0,dM,CS,G_Rho0,e,use_p_atm,use_EOS,G,pa, & -!$OMP rho_ref,js,je,is,ie,intx_pa,inty_pa) +!$OMP p0,dM,CS,G_Rho0,e,use_p_atm,use_EOS,G, & +!$OMP rho_ref,js,je,is,ie) if (CS%GFS_scale < 1.0) then ! Adjust the Montgomery potential to make this a reduced gravity model. if (use_EOS) then @@ -611,29 +659,6 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, CS, ALE_CSp, p_atm, pbce, enddo ; enddo endif endif - - ! Set the surface boundary conditions on pressure anomaly and its horizontal - ! integrals, assuming that the surface pressure anomaly varies linearly - ! in x and y. - if (use_p_atm) then -!$OMP do - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - pa(i,j) = (rho_ref*G%g_Earth)*e(i,j,1) + p_atm(i,j) - enddo ; enddo - else -!$OMP do - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - pa(i,j) = (rho_ref*G%g_Earth)*e(i,j,1) - enddo ; enddo - endif -!$OMP do - do j=js,je ; do I=Isq,Ieq - intx_pa(I,j) = 0.5*(pa(i,j) + pa(i+1,j)) - enddo ; enddo -!$OMP do - do J=Jsq,Jeq ; do i=is,ie - inty_pa(i,J) = 0.5*(pa(i,j) + pa(i,j+1)) - enddo ; enddo !$OMP end parallel ! Have checked that rho_0 drops out and that the 1-layer case is right. RWH. @@ -653,94 +678,109 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, CS, ALE_CSp, p_atm, pbce, endif endif -!$OMP parallel default(none) shared(nz,Jsq,Jeq,Isq,Ieq,G,h,use_EOS,use_ALE,PRScheme,& -!$OMP T_t,T_b,S_t,S_b,e,rho_ref,CS,tv,js,je,is,ie, & -!$OMP dpa,intz_dpa,intx_dpa,inty_dpa,tv_tmp,I_Rho0,pa, & -!$OMP pa_h_intz_dpa,PFu,PFv,intx_pa,h_neglect,inty_pa) & -!$OMP private(dz) -!$OMP do +!$OMP parallel do default(none) shared(nblocks,block,G_BLK,use_p_atm,rho_ref,G,e, & +!$OMP p_atm,nz,use_EOS,use_ALE,PRScheme,T_t,T_b,S_t, & +!$OMP S_b,CS,tv,tv_tmp,h,PFu,I_Rho0,h_neglect,PFv,dM)& +!$OMP private(is,ie,js,je,Isq,Ieq,Jsq,Jeq,isd,ied,jsd,jed) +do n = 1, nblocks + is=block(n)%isc; ie=block(n)%iec; js=block(n)%jsc; je=block(n)%jec + Isq=block(n)%IscB; Ieq=block(n)%IecB; Jsq=block(n)%JscB; Jeq=block(n)%JecB + isd=block(n)%isd; ied=block(n)%ied; jsd=block(n)%jsd; jed=block(n)%jed + + ! Set the surface boundary conditions on pressure anomaly and its horizontal + ! integrals, assuming that the surface pressure anomaly varies linearly + ! in x and y. + if (use_p_atm) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + block(n)%pa(i,j) = (rho_ref*G%g_Earth)*e(i,j,1) + p_atm(i,j) + enddo ; enddo + else + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + block(n)%pa(i,j) = (rho_ref*G%g_Earth)*e(i,j,1) + enddo ; enddo + endif + do j=js,je ; do I=Isq,Ieq + block(n)%intx_pa(I,j) = 0.5*(block(n)%pa(i,j) + block(n)%pa(i+1,j)) + enddo ; enddo + do J=Jsq,Jeq ; do i=is,ie + block(n)%inty_pa(i,J) = 0.5*(block(n)%pa(i,j) + block(n)%pa(i,j+1)) + enddo ; enddo + do k=1,nz ! Calculate 4 integrals through the layer that are required in the ! subsequent calculation. - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - dz(i,j) = G%g_Earth*G%H_to_m*h(i,j,k) - enddo ; enddo + if (use_EOS) then - ! The following routine computes the integrals that are needed to - ! calculate the pressure gradient force. Linear profiles for T and S are - ! assumed when regridding is activated. Otherwise, the previous version - ! is used, whereby densities within each layer are constant no matter - ! where the layers are located. - if ( use_ALE ) then - if ( PRScheme == PRESSURE_RECONSTRUCTION_PLM ) then - call int_density_dz_generic_plm ( T_t(:,:,k), T_b(:,:,k), S_t(:,:,k), & - S_b(:,:,k), & - e(:,:,K), e(:,:,K+1), rho_ref, & - CS%Rho0, G%g_Earth, & - G, tv%eqn_of_state, dpa(:,:,k), intz_dpa(:,:,k), & - intx_dpa(:,:,k), inty_dpa(:,:,k), & - useMassWghtInterp = CS%useMassWghtInterp) - elseif ( PRScheme == PRESSURE_RECONSTRUCTION_PPM ) then - call int_density_dz_generic_ppm ( tv%T(:,:,k), T_t(:,:,k), T_b(:,:,k), & - tv%S(:,:,k), S_t(:,:,k), S_b(:,:,k), & - e(:,:,K), e(:,:,K+1), rho_ref, & - CS%Rho0, G%G_Earth, & - G, tv%eqn_of_state, dpa(:,:,k), intz_dpa(:,:,k), & - intx_dpa(:,:,k), inty_dpa(:,:,k)) - endif - else - call int_density_dz(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), e(:,:,K), e(:,:,K+1), & - rho_ref, CS%Rho0, G%g_Earth, G, tv%eqn_of_state, & - dpa(:,:,k), intz_dpa(:,:,k), intx_dpa(:,:,k), inty_dpa(:,:,k)) - endif + ! The following routine computes the integrals that are needed to + ! calculate the pressure gradient force. Linear profiles for T and S are + ! assumed when regridding is activated. Otherwise, the previous version + ! is used, whereby densities within each layer are constant no matter + ! where the layers are located. + if ( use_ALE ) then + if ( PRScheme == PRESSURE_RECONSTRUCTION_PLM ) then + call int_density_dz_generic_plm ( T_t(isd:ied,jsd:jed,k), & + T_b(isd:ied,jsd:jed,k), S_t(isd:ied,jsd:jed,k), & + S_b(isd:ied,jsd:jed,k), e(isd:ied,jsd:jed,K), & + e(isd:ied,jsd:jed,K+1), rho_ref, CS%Rho0, G%g_Earth, & + G_BLK(n), tv%eqn_of_state, Block(n)%dpa, Block(n)%intz_dpa, & + Block(n)%intx_dpa, Block(n)%inty_dpa, & + useMassWghtInterp = CS%useMassWghtInterp) + elseif ( PRScheme == PRESSURE_RECONSTRUCTION_PPM ) then + call int_density_dz_generic_ppm ( tv%T(isd:ied,jsd:jed,k), & + T_t(isd:ied,jsd:jed,k), T_b(isd:ied,jsd:jed,k), & + tv%S(isd:ied,jsd:jed,k), S_t(isd:ied,jsd:jed,k), & + S_b(isd:ied,jsd:jed,k), e(isd:ied,jsd:jed,K), & + e(isd:ied,jsd:jed,K+1), rho_ref, CS%Rho0, G%G_Earth, & + G, tv%eqn_of_state, Block(n)%dpa, Block(n)%intz_dpa, & + Block(n)%intx_dpa, Block(n)%inty_dpa) + endif + else + call int_density_dz(tv_tmp%T(isd:ied,jsd:jed,k), tv_tmp%S(isd:ied,jsd:jed,k), & + e(isd:ied,jsd:jed,K), e(isd:ied,jsd:jed,K+1), & + rho_ref, CS%Rho0, G%g_Earth, G, tv%eqn_of_state, & + Block(n)%dpa, Block(n)%intz_dpa, Block(n)%intx_dpa, Block(n)%inty_dpa) + endif else do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - dpa(i,j,k) = (G%Rlay(k) - rho_ref)*dz(i,j) - intz_dpa(i,j,k) = 0.5*(G%Rlay(k) - rho_ref)*dz(i,j)*h(i,j,k) + Block(n)%dz(i,j) = G%g_Earth*G%H_to_m*h(i,j,k) + Block(n)%dpa(i,j) = (G%Rlay(k) - rho_ref)*Block(n)%dz(i,j) + Block(n)%intz_dpa(i,j) = 0.5*(G%Rlay(k) - rho_ref)*Block(n)%dz(i,j)*h(i,j,k) enddo ; enddo do j=js,je ; do I=Isq,Ieq - intx_dpa(I,j,k) = 0.5*(G%Rlay(k) - rho_ref) * (dz(i,j)+dz(i+1,j)) + Block(n)%intx_dpa(I,j) = 0.5*(G%Rlay(k) - rho_ref) * (Block(n)%dz(i,j)+Block(n)%dz(i+1,j)) enddo ; enddo do J=Jsq,Jeq ; do i=is,ie - inty_dpa(i,J,k) = 0.5*(G%Rlay(k) - rho_ref) * (dz(i,j)+dz(i,j+1)) + Block(n)%inty_dpa(i,J) = 0.5*(G%Rlay(k) - rho_ref) * (Block(n)%dz(i,j)+Block(n)%dz(i,j+1)) enddo ; enddo endif - enddo -!$OMP do - do j=jsq,Jeq+1; do k=1,nz; do i=Isq,Ieq+1 - pa_h_intz_dpa(i,j,k) = pa(i,j)*h(i,j,k) + intz_dpa(i,j,k) - pa(i,j) = pa(i,j)+dpa(i,j,k) - enddo; enddo; enddo - - ! Compute pressure gradient in x direction -!$OMP do - do j =js,je - do k = 1,nz; do I=Isq,Ieq - PFu(I,j,k) = ((pa_h_intz_dpa(i,j,k) - pa_h_intz_dpa(i+1,j,k)) + & - ((h(i+1,j,k) - h(i,j,k)) * intx_pa(I,j) - & - (e(i+1,j,K+1) - e(i,j,K+1)) * intx_dpa(I,j,k))) * & + ! Compute pressure gradient in x direction + do j=js,je ; do I=Isq,Ieq + PFu(I,j,k) = (((Block(n)%pa(i,j)*h(i,j,k) + Block(n)%intz_dpa(i,j)) - & + (Block(n)%pa(i+1,j)*h(i+1,j,k) + Block(n)%intz_dpa(i+1,j))) + & + ((h(i+1,j,k) - h(i,j,k)) * Block(n)%intx_pa(I,j) - & + (e(i+1,j,K+1) - e(i,j,K+1)) * Block(n)%intx_dpa(I,j))) * & ((2.0*I_Rho0*G%IdxCu(I,j)) / & ((h(i,j,k) + h(i+1,j,k)) + h_neglect)) - intx_pa(I,j) = intx_pa(I,j) + intx_dpa(I,j,k) + Block(n)%intx_pa(I,j) = Block(n)%intx_pa(I,j) + Block(n)%intx_dpa(I,j) enddo ; enddo - enddo ! Compute pressure gradient in y direction -!$OMP do - do J=Jsq,Jeq - do k=1,nz ; do i=is,ie - PFv(i,J,k) = ((pa_h_intz_dpa(i,j,k) - pa_h_intz_dpa(i,j+1,k)) + & - ((h(i,j+1,k) - h(i,j,k)) * inty_pa(i,J) - & - (e(i,j+1,K+1) - e(i,j,K+1)) * inty_dpa(i,J,k))) * & + do J=Jsq,Jeq ; do i=is,ie + PFv(i,J,k) = (((Block(n)%pa(i,j)*h(i,j,k) + Block(n)%intz_dpa(i,j)) - & + (Block(n)%pa(i,j+1)*h(i,j+1,k) + Block(n)%intz_dpa(i,j+1))) + & + ((h(i,j+1,k) - h(i,j,k)) * Block(n)%inty_pa(i,J) - & + (e(i,j+1,K+1) - e(i,j,K+1)) * Block(n)%inty_dpa(i,J))) * & ((2.0*I_Rho0*G%IdyCv(i,J)) / & ((h(i,j,k) + h(i,j+1,k)) + h_neglect)) - inty_pa(i,J) = inty_pa(i,J) + inty_dpa(i,J,k) + Block(n)%inty_pa(i,J) = Block(n)%inty_pa(i,J) + Block(n)%inty_dpa(i,J) + enddo ; enddo + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + Block(n)%pa(i,j) = Block(n)%pa(i,j) + Block(n)%dpa(i,j) enddo ; enddo enddo -!$OMP end parallel + if (CS%GFS_scale < 1.0) then -!$OMP parallel do default(none) shared(is,ie,js,je,nz,Isq,Ieq,Jsq,Jeq,PFu,PFv,dM,G) do k=1,nz do j=js,je ; do I=Isq,Ieq PFu(I,j,k) = PFu(I,j,k) - (dM(i+1,j) - dM(i,j)) * G%IdxCu(I,j) @@ -750,12 +790,15 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, CS, ALE_CSp, p_atm, pbce, enddo ; enddo enddo endif +enddo + if (present(pbce)) then call set_pbce_Bouss(e, tv_tmp, G, G%g_Earth, CS%Rho0, CS%GFS_scale, pbce) endif if (present(eta)) then + Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB if (CS%tides) then ! eta is the sea surface height relative to a time-invariant geoid, for ! comparison with what is used for eta in btstep. See how e was calculated @@ -794,7 +837,14 @@ subroutine PressureForce_AFV_init(Time, G, param_file, diag, CS, tides_CSp) ! (in) tides_CSp - a pointer to the control structure of the tide module. ! This include declares and sets the variable "version". #include "version_variable.h" - character(len=40) :: mod ! This module's name. + character(len=40) :: modname ! This module's name. + integer, allocatable :: ibegin(:), iend(:), jbegin(:), jend(:) + integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB + integer :: isc, iec, jsc, jec, IscB, IecB, JscB, JecB + integer :: n, halo, i, j + integer :: xblock = 1 !number of blocks in x-direction on each processor(for openmp) + integer :: yblock = 1 !number of blocks in y-direction on each processor(for openmp) + if (associated(CS)) then call MOM_error(WARNING, "PressureForce_init called with an associated "// & @@ -807,17 +857,17 @@ subroutine PressureForce_AFV_init(Time, G, param_file, diag, CS, tides_CSp) if (associated(tides_CSp)) CS%tides_CSp => tides_CSp endif - mod = "MOM_PressureForce_AFV" - call log_version(param_file, mod, version) - call get_param(param_file, mod, "RHO_0", CS%Rho0, & + modname = "MOM_PressureForce_AFV" + call log_version(param_file, modname, version) + call get_param(param_file, modname, "RHO_0", CS%Rho0, & "The mean ocean density used with BOUSSINESQ true to \n"//& "calculate accelerations and the mass for conservation \n"//& "properties, or with BOUSSINSEQ false to convert some \n"//& "parameters from vertical units of m to kg m-2.", & units="kg m-3", default=1035.0) - call get_param(param_file, mod, "TIDES", CS%tides, & + call get_param(param_file, modname, "TIDES", CS%tides, & "If true, apply tidal momentum forcing.", default=.false.) - call get_param(param_file, mod, "MASS_WEIGHT_IN_PRESSURE_GRADIENT", CS%useMassWghtInterp, & + call get_param(param_file, modname, "MASS_WEIGHT_IN_PRESSURE_GRADIENT", CS%useMassWghtInterp, & "If true, use mass weighting when interpolation T/S for\n"//& "top/bottom integrals in AFV pressure gradient calculation.", default=.false.) @@ -829,7 +879,79 @@ subroutine PressureForce_AFV_init(Time, G, param_file, diag, CS, tides_CSp) CS%GFS_scale = 1.0 if (G%g_prime(1) /= G%g_Earth) CS%GFS_scale = G%g_prime(1) / G%g_Earth - call log_param(param_file, mod, "GFS / G_EARTH", CS%GFS_scale) + call log_param(param_file, modname, "GFS / G_EARTH", CS%GFS_scale) + + call get_param(param_file, modname, "NI_BLOCK", xblock, "The number of blocks "// & + "in the x-direction on each processor(for openmp).", default=1) + call get_param(param_file, modname, "NJ_BLOCK", yblock, "The number of blocks "// & + "in the y-direction on each processor(for openmp).", default=1) + + halo = G%isc - G%isd + nblocks = xblock * yblock + allocate(ibegin(xblock), iend(xblock), jbegin(yblock), jend(yblock)) + allocate(block(nblocks)) + if (G%Boussinesq) allocate(G_BLK(nblocks)) + call compute_extent(G%isc,G%iec,xblock,ibegin,iend) + call compute_extent(G%jsc,G%jec,yblock,jbegin,jend) + do n = 1,nblocks + i = mod((n-1), xblock) + 1 + j = (n-1)/xblock + 1 + isc = ibegin(i) ; iec = iend(i) ; jsc = jbegin(j) ; jec = jend(j) + IscB = isc ; IecB = iec ; JscB = jsc ; JecB = jec + !--- For symmetry domain, the first block will have the extra point. + if (G%symmetric) then + if (i==1) IscB = IscB-1 + if (j==1) JscB = JscB-1 + endif + isd = isc - halo ; ied = iec + halo + jsd = jsc - halo ; jed = jec + halo + IsdB = IscB - halo ; IedB = IecB + halo + JsdB = JscB - halo ; JedB = JecB + halo + block(n)%isc = isc ; block(n)%iec = iec + block(n)%jsc = jsc ; block(n)%jec = jec + block(n)%isd = isd ; block(n)%ied = ied + block(n)%jsd = jsd ; block(n)%jed = jed + block(n)%IscB = IscB ; block(n)%IecB = IecB + block(n)%JscB = JscB ; block(n)%JecB = JecB + if (G%Boussinesq) then + G_BLK(n)%H_subroundoff = G%H_subroundoff + G_BLK(n)%isc = isc-isd+1 ; G_BLK(n)%iec = iec-isd+1 + G_BLK(n)%jsc = jsc-jsd+1 ; G_BLK(n)%jec = jec-jsd+1 + G_BLK(n)%IscB = IscB-isd+1 ; G_BLK(n)%IecB = IecB-isd+1 + G_BLK(n)%JscB = JscB-jsd+1 ; G_BLK(n)%JecB = JecB-jsd+1 + G_BLK(n)%isd = G_BLK(n)%isc - halo + G_BLK(n)%ied = G_BLK(n)%iec + halo + G_BLK(n)%jsd = G_BLK(n)%jsc - halo + G_BLK(n)%jed = G_BLK(n)%jec + halo + G_BLK(n)%IsdB = G_BLK(n)%IscB - halo + G_BLK(n)%IedB = G_BLK(n)%Iecb + halo + G_BLK(n)%Jsdb = G_BLK(n)%JscB - halo + G_BLK(n)%Jedb = G_BLK(n)%JecB + halo + allocate(G_BLK(n)%bathyT(G_BLK(n)%isd:G_BLK(n)%ied,G_BLK(n)%jsd:G_BLK(n)%jed)) + G_BLK(n)%bathyT(:,:) = G%bathyT(isd:ied,jsd:jed) + + allocate(block(n)%dz (isd :ied ,jsd :jed )) + allocate(block(n)%pa (isd :ied ,jsd :jed )) + allocate(block(n)%intx_pa (IsdB:IedB,jsd :jed )) + allocate(block(n)%inty_pa (isd :ied ,JsdB:JedB)) + allocate(block(n)%dpa (isd :ied ,jsd :jed )) + allocate(block(n)%intx_dpa(IsdB:IedB,jsd :jed )) + allocate(block(n)%inty_dpa(isd :ied ,JsdB:JedB)) + allocate(block(n)%intz_dpa(isd :ied ,jsd :jed )) + block(n)%dz = 0 ; block(n)%pa = 0 + block(n)%intx_pa = 0 ; block(n)%inty_pa = 0 + block(n)%dpa = 0 ; block(n)%intx_dpa = 0 + block(n)%inty_dpa = 0 ; block(n)%intz_dpa = 0 + else + allocate(block(n)%dp (isd :ied ,jsd :jed )) + allocate(block(n)%za (isd :ied ,jsd :jed )) + allocate(block(n)%intx_za(IsdB:IedB,jsd :jed )) + allocate(block(n)%inty_za(isd :ied ,JsdB:JedB)) + block(n)%dp = 0 ; block(n)%za = 0 + block(n)%intx_za = 0 ; block(n)%inty_za = 0 + endif + enddo + deallocate(ibegin, iend, jbegin, jend) end subroutine PressureForce_AFV_init diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index eb6ebd8a1c..b6e6b5e43e 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -1010,13 +1010,13 @@ subroutine int_density_dz_generic_plm (T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, & ! 1. Compute vertical integrals ! ============================= do j=Jsq,Jeq+1 - do i = Isq,Ieq+1 - dz(i) = z_t(i,j) - z_b(i,j) - do n=1,5 - p5(i*5+n) = -GxRho*(z_t(i,j) - 0.25*real(n-1)*dz(i)) - ! Salinity and temperature points are linearly interpolated - S5(i*5+n) = wt_t(n) * S_t(i,j) + wt_b(n) * S_b(i,j) - T5(i*5+n) = wt_t(n) * T_t(i,j) + wt_b(n) * T_b(i,j) + do i = Isq,Ieq+1 + dz(i) = z_t(i,j) - z_b(i,j) + do n=1,5 + p5(i*5+n) = -GxRho*(z_t(i,j) - 0.25*real(n-1)*dz(i)) + ! Salinity and temperature points are linearly interpolated + S5(i*5+n) = wt_t(n) * S_t(i,j) + wt_b(n) * S_b(i,j) + T5(i*5+n) = wt_t(n) * T_t(i,j) + wt_b(n) * T_b(i,j) enddo enddo call calculate_density_array(T5, S5, p5, r5, 1, (ieq-isq+2)*5, EOS ) @@ -1024,15 +1024,15 @@ subroutine int_density_dz_generic_plm (T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, & do i = isq, ieq+1 ! Use Bode's rule to estimate the pressure anomaly change. - rho_anom = C1_90*(7.0*(r5(i*5+1)+r5(i*5+5)) + 32.0*(r5(i*5+2)+r5(i*5+4)) + 12.0*r5(i*5+3)) - & + rho_anom = C1_90*(7.0*(r5(i*5+1)+r5(i*5+5)) + 32.0*(r5(i*5+2)+r5(i*5+4)) + 12.0*r5(i*5+3)) - & rho_ref - dpa(i,j) = G_e*dz(i)*rho_anom - if(present(intz_dpa)) then - ! Use a Bode's-rule-like fifth-order accurate estimate of - ! the double integral of the pressure anomaly. - intz_dpa(i,j) = 0.5*G_e*dz(i)**2 * & + dpa(i,j) = G_e*dz(i)*rho_anom + if (present(intz_dpa)) then + ! Use a Bode's-rule-like fifth-order accurate estimate of + ! the double integral of the pressure anomaly. + intz_dpa(i,j) = 0.5*G_e*dz(i)**2 * & (rho_anom - C1_90*(16.0*(u5(i*5+4)-u5(i*5+2)) + 7.0*(u5(i*5+5)-u5(i*5+1))) ) - endif + endif enddo enddo ! end loops on j diff --git a/src/framework/MOM_domains.F90 b/src/framework/MOM_domains.F90 index 2c8b33cdde..b2f54d8d56 100644 --- a/src/framework/MOM_domains.F90 +++ b/src/framework/MOM_domains.F90 @@ -40,6 +40,7 @@ module MOM_domains use mpp_domains_mod, only : mpp_reset_group_update_field use mpp_domains_mod, only : mpp_group_update_initialized use mpp_domains_mod, only : mpp_start_group_update, mpp_complete_group_update +use mpp_domains_mod, only : compute_extent => mpp_compute_extent use mpp_parameter_mod, only : AGRID, BGRID_NE, CGRID_NE, SCALAR_PAIR, BITWISE_EXACT_SUM, CORNER use mpp_parameter_mod, only : To_East => WUPDATE, To_West => EUPDATE use mpp_parameter_mod, only : To_North => SUPDATE, To_South => NUPDATE @@ -59,6 +60,7 @@ module MOM_domains public :: To_East, To_West, To_North, To_South, To_All public :: create_group_pass, do_group_pass, group_pass_type public :: start_group_pass, complete_group_pass +public :: compute_extent interface pass_var module procedure pass_var_3d, pass_var_2d