Skip to content
New issue

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

Sign up for GitHub

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

Already on GitHub? Sign in to your account

Jump to bottom

User/z1l/openmp #24

Merged
merged 4 commits into from
Jun 16, 2014
Merged

User/z1l/openmp #24

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
8760bfc
Implement Openmp
Zhi-Liang Jun 9, 2014
2294a22
Merge branch 'user/z1l/openmp' of github.com:CommerceGov/NOAA-GFDL-MO…
Zhi-Liang Jun 9, 2014
07bbe03
Merge remote-tracking branch 'origin/dev/master' into user/z1l/openmp
Zhi-Liang Jun 9, 2014
6a04d99
fix for openmp
Zhi-Liang Jun 9, 2014
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
22 changes: 13 additions & 9 deletions src/core/MOM_CoriolisAdv.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -283,23 +283,26 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, AD, G, CS)
Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB ; nz = G%ke
h_neglect = G%H_subroundoff

!$OMP parallel default(none) shared(u,v,h,uh,vh,CAu,CAv,G,CS,AD,Area_h,Area_q,nz,RV,PV, &
!$OMP is,ie,js,je,Isq,Ieq,Jsq,Jeq,h_neglect) &
!$OMP private(relative_vorticity,absolute_vorticity,Ih,hArea_q,q, &
!$OMP abs_vort,Ih_q,fv1,fv2,fu1,fu2,max_fvq,max_fuq, &
!$OMP min_fvq,min_fuq,q2,a,b,c,d,ep_u,ep_v,Fe_m2,rat_lin, &
!$OMP min_Ihq,max_Ihq,rat_m1,AL_wt,Sad_wt,c1,c2,c3,slope, &
!$OMP uhc,uhm,uh_min,uh_max,vhc,vhm,vh_min,vh_max,KE,KEx, &
!$OMP KEy,temp1,temp2,Heff1,Heff2,Heff3,Heff4,VHeff, &
!$OMP QVHeff,max_fv,min_fv,UHeff,QUHeff,max_fu,min_fu)
!$OMP do
do j=Jsq-1,Jeq+2 ; do I=Isq-1,Ieq+2
Area_h(i,j) = G%mask2dT(i,j) * G%areaT(i,j)
enddo ; enddo
!$OMP do
do J=Jsq-1,Jeq+1 ; do I=Isq-1,Ieq+1
Area_q(i,j) = (Area_h(i,j) + Area_h(i+1,j+1)) + &
(Area_h(i+1,j) + Area_h(i,j+1))
enddo ; enddo

!$OMP parallel do default(none) shared(u,v,h,uh,vh,CAu,CAv,G,CS,AD,Area_h,Area_q,nz,RV,PV, &
!$OMP is,ie,js,je,Isq,Ieq,Jsq,Jeq,h_neglect) &
!$OMP private(relative_vorticity,absolute_vorticity,Ih,hArea_q,q, &
!$OMP abs_vort,Ih_q,fv1,fv2,fu1,fu2,max_fvq,max_fuq, &
!$OMP min_fvq,min_fuq,q2,a,b,c,d,ep_u,ep_v,Fe_m2,rat_lin, &
!$OMP min_Ihq,max_Ihq,rat_m1,AL_wt,Sad_wt,c1,c2,c3,slope, &
!$OMP uhc,uhm,uh_min,uh_max,vhc,vhm,vh_min,vh_max,KE,KEx, &
!$OMP KEy,temp1,temp2,Heff1,Heff2,Heff3,Heff4,VHeff, &
!$OMP QVHeff,max_fv,min_fv,UHeff,QUHeff,max_fu,min_fu)
!$OMP do
do k=1,nz
! Here the second order accurate layer potential vorticities, q,
! are calculated. hq is second order accurate in space. Relative
Expand Down Expand Up @@ -727,6 +730,7 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, AD, G, CS)
endif

enddo ! k-loop.
!$OMP end parallel

! Here the various Coriolis-related derived quantities are offered
! for averaging.
Expand Down
89 changes: 57 additions & 32 deletions src/core/MOM_PressureForce_Montgomery.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -469,28 +469,34 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta)
! and loading. This should really be based on bottom pressure anomalies,
! but that is not yet implemented, and the current form is correct for
! barotropic tides.
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
e(i,j,1) = -1.0*G%bathyT(i,j)
enddo ; enddo
do k=1,nz ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
e(i,j,1) = e(i,j,1) + h(i,j,k)
enddo ; enddo ; enddo
!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,e,h,G)
do j=Jsq,Jeq+1
do i=Isq,Ieq+1 ; e(i,j,1) = -1.0*G%bathyT(i,j) ; enddo
do k=1,nz ; do i=Isq,Ieq+1
e(i,j,1) = e(i,j,1) + h(i,j,k)
enddo ; enddo
enddo
call calc_tidal_forcing(CS%Time, e(:,:,1), e_tidal, G, CS%tides_CSp)
endif

! Here layer interface heights, e, are calculated.
!$OMP parallel default(none) shared(Isq,Ieq,Jsq,Jeq,nz,e,h,G,e_tidal,CS)
if (CS%tides) then
!$OMP do
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
e(i,j,nz+1) = -1.0*G%bathyT(i,j) - e_tidal(i,j)
enddo ; enddo
else
!$OMP do
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
e(i,j,nz+1) = -1.0*G%bathyT(i,j)
enddo ; enddo
endif
do k=nz,1,-1 ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
!$OMP do
do j=Jsq,Jeq+1 ; do k=nz,1,-1 ; do i=Isq,Ieq+1
e(i,j,K) = e(i,j,K+1) + h(i,j,k)
enddo ; enddo ; enddo
!$OMP end parallel
if (use_EOS) then

! Calculate in-situ densities (rho_star).
Expand All @@ -503,28 +509,34 @@ subroutine PressureForce_Mont_Bouss(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
do j=Jsq,Jeq+1
!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,nkmb,tv_tmp,tv,p_ref, &
!$OMP Rho_cv_BL,G)
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)

do k=nkmb+1,nz ; 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
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.0 ; enddo
endif

! This no longer includes any pressure dependency, since this routine
! will come down with a fatal error if there is any compressibility.
!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,tv_tmp,p_ref,rho_star,tv,G_Rho0)
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), &
Isq,Ieq-Isq+2,tv%eqn_of_state)
Expand All @@ -534,21 +546,28 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta)

! Here the layer Montgomery potentials, M, are calculated.
if (use_EOS) then
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
M(i,j,1) = CS%GFS_scale * (rho_star(i,j,1) * e(i,j,1))
if (use_p_atm) M(i,j,1) = M(i,j,1) + p_atm(i,j) * I_Rho0
enddo ; enddo
do k=2,nz ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
M(i,j,k) = M(i,j,k-1) + (rho_star(i,j,k) - rho_star(i,j,k-1)) * e(i,j,K)
enddo ; enddo ; enddo
!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,M,CS,rho_star,e,use_p_atm, &
!$OMP p_atm,I_Rho0)
do j=Jsq,Jeq+1
do i=Isq,Ieq+1
M(i,j,1) = CS%GFS_scale * (rho_star(i,j,1) * e(i,j,1))
if (use_p_atm) M(i,j,1) = M(i,j,1) + p_atm(i,j) * I_Rho0
enddo
do k=2,nz ; do i=Isq,Ieq+1
M(i,j,k) = M(i,j,k-1) + (rho_star(i,j,k) - rho_star(i,j,k-1)) * e(i,j,K)
enddo ; enddo
enddo
else ! not use_EOS
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
M(i,j,1) = G%g_prime(1) * e(i,j,1)
if (use_p_atm) M(i,j,1) = M(i,j,1) + p_atm(i,j) * I_Rho0
enddo ; enddo
do k=2,nz ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
M(i,j,k) = M(i,j,k-1) + G%g_prime(K) * e(i,j,K)
enddo ; enddo ; enddo
!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,M,G,e,use_p_atm,p_atm,I_Rho0)
do j=Jsq,Jeq+1
do i=Isq,Ieq+1
M(i,j,1) = G%g_prime(1) * e(i,j,1)
if (use_p_atm) M(i,j,1) = M(i,j,1) + p_atm(i,j) * I_Rho0
enddo
do k=2,nz ; do i=Isq,Ieq+1
M(i,j,k) = M(i,j,k-1) + G%g_prime(K) * e(i,j,K)
enddo ; enddo
enddo
endif ! use_EOS

if (present(pbce)) then
Expand All @@ -559,6 +578,9 @@ subroutine PressureForce_Mont_Bouss(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(Isq,Ieq,Jsq,Jeq,js,je,is,ie,nz,e,h_neglect, &
!$OMP rho_star,G,PFu,CS,PFv,M) &
!$OMP private(h_star,PFu_bc,PFv_bc)
do k=1,nz
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
h_star(i,j) = (e(i,j,K) - e(i,j,K+1)) + h_neglect
Expand All @@ -579,6 +601,7 @@ subroutine PressureForce_Mont_Bouss(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(Isq,Ieq,Jsq,Jeq,is,ie,js,je,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)
Expand All @@ -594,10 +617,12 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta)
! 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
! about 200 lines above.
!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,eta,e,e_tidal)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
eta(i,j) = e(i,j,1) + e_tidal(i,j)
enddo ; enddo
else
!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,eta,e)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
eta(i,j) = e(i,j,1)
enddo ; enddo
Expand Down
Loading