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 #12

Merged
merged 15 commits into from
Apr 15, 2014
Merged

User/z1l/openmp #12

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
feb8e15
make visc_bound_rem to be private to solve slowdown issue
Zhi-Liang Apr 8, 2014
5c741f6
implement openmp
Zhi-Liang Apr 8, 2014
b70a603
implement openmp
Zhi-Liang Apr 8, 2014
afd1b07
implement openmp
Zhi-Liang Apr 9, 2014
d93cacd
change loop structure and implement openmp
Zhi-Liang Apr 9, 2014
287ac53
rearrange loop structure and implement openmp
Zhi-Liang Apr 9, 2014
93dd435
Add more OMP directives
Zhi-Liang Apr 14, 2014
3f6a0c3
implement openmp
Zhi-Liang Apr 14, 2014
662d47e
implement openmp
Zhi-Liang Apr 14, 2014
b68270e
Add more OMP directives
Zhi-Liang Apr 14, 2014
b3a9a2f
implement openmp
Zhi-Liang Apr 14, 2014
dd7a799
implement openmp
Zhi-Liang Apr 14, 2014
5d742b8
implement openmp
Zhi-Liang Apr 14, 2014
2c162fc
openmp bug fix
Zhi-Liang Apr 14, 2014
985a905
Merge remote-tracking branch 'origin/dev/master' into user/z1l/openmp
Zhi-Liang Apr 14, 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
2 changes: 1 addition & 1 deletion src/core/MOM_legacy_barotropic.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3356,7 +3356,7 @@ subroutine legacy_bt_mass_source(h, eta, fluxes, set_cor, dt_therm, &

is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke

!$OMP parallel do default(shared) private(i,j,k,eta_h)
!$OMP parallel do default(shared) private(eta_h,h_tot,limit_dt,d_eta)
do j=js,je
do i=is,ie ; h_tot(i) = h(i,j,1) ; enddo
if (G%Boussinesq) then
Expand Down
143 changes: 83 additions & 60 deletions src/diagnostics/MOM_diagnostics.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -141,7 +141,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, dt, G, &
type(cont_diag_ptrs), intent(in) :: CDp
real, intent(in) :: dt
type(ocean_grid_type), intent(inout) :: G
type(diagnostics_CS), pointer :: CS
type(diagnostics_CS), intent(inout) :: CS
real, dimension(NIMEM_,NJMEM_), optional, intent(in) :: eta_bt
! Any diagnostic quantities that are not more naturally calculated
! in the various other subroutines are calculated here.
Expand Down Expand Up @@ -188,7 +188,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, dt, G, &
! one iteration which would break the follwoing one-line workaround!
if (nkmb==0) nkmb = nz

if (.not.associated(CS)) call MOM_error(FATAL, &
if (loc(CS)==0) call MOM_error(FATAL, &
"calculate_diagnostic_fields: Module must be initialized before it is used.")

call calculate_derivs(dt, G, CS)
Expand Down Expand Up @@ -222,6 +222,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, dt, G, &

if (associated(tv%eqn_of_state)) then
pres(:) = tv%P_Ref
!$OMP parallel do default(shared)
do k=1,nz ; do j=js,je+1
call calculate_density(tv%T(:,j,k),tv%S(:,j,k),pres, &
Rcv(:,j,k),is,ie-is+2, tv%eqn_of_state)
Expand All @@ -231,82 +232,103 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, dt, G, &
endif

if (ASSOCIATED(CS%h_Rlay)) then
do k=1,nkmb ; do j=js,je; do i=is,ie
CS%h_Rlay(i,j,k) = 0.0
enddo ; enddo ; enddo
do k=nkmb+1,nz ; do j=js,je ; do i=is,ie
CS%h_Rlay(i,j,k) = h(i,j,k)
enddo ; enddo ; enddo
k_list = nz/2
do j=js,je ; do k=1,nkmb ; do i=is,ie
call find_weights(G%Rlay, Rcv(i,j,k), k_list, nz, wt, wt_p)
CS%h_Rlay(i,j,k_list) = CS%h_Rlay(i,j,k_list) + h(i,j,k)*wt
CS%h_Rlay(i,j,k_list+1) = CS%h_Rlay(i,j,k_list+1) + h(i,j,k)*wt_p
enddo ; enddo ; enddo
!$OMP parallel do default(shared) private(wt,wt_p) firstprivate(k_list)
do j=js,je
do k=1,nkmb ; do i=is,ie
CS%h_Rlay(i,j,k) = 0.0
enddo ; enddo
do k=nkmb+1,nz ; do i=is,ie
CS%h_Rlay(i,j,k) = h(i,j,k)
enddo ; enddo
do k=1,nkmb ; do i=is,ie
call find_weights(G%Rlay, Rcv(i,j,k), k_list, nz, wt, wt_p)
CS%h_Rlay(i,j,k_list) = CS%h_Rlay(i,j,k_list) + h(i,j,k)*wt
CS%h_Rlay(i,j,k_list+1) = CS%h_Rlay(i,j,k_list+1) + h(i,j,k)*wt_p
enddo ; enddo
enddo

if (CS%id_h_Rlay > 0) call post_data(CS%id_h_Rlay, CS%h_Rlay, CS%diag)
endif

if (ASSOCIATED(CS%uh_Rlay)) then
do k=1,nkmb ; do j=js,je; do I=Isq,Ieq
CS%uh_Rlay(I,j,k) = 0.0
enddo ; enddo ; enddo
do k=nkmb+1,nz ; do j=js,je ; do I=Isq,Ieq
CS%uh_Rlay(I,j,k) = uh(I,j,k)
enddo ; enddo ; enddo
k_list = nz/2
do j=js,je ; do k=1,nkmb ; do I=Isq,Ieq
call find_weights(G%Rlay, 0.5*(Rcv(i,j,k)+Rcv(i+1,j,k)), k_list, nz, wt, wt_p)
CS%uh_Rlay(I,j,k_list) = CS%uh_Rlay(I,j,k_list) + uh(I,j,k)*wt
CS%uh_Rlay(I,j,k_list+1) = CS%uh_Rlay(I,j,k_list+1) + uh(I,j,k)*wt_p
enddo ; enddo ; enddo
!$OMP parallel do default(shared) private(wt,wt_p) firstprivate(k_list)
do j=js,je
do k=1,nkmb ; do I=Isq,Ieq
CS%uh_Rlay(I,j,k) = 0.0
enddo ; enddo
do k=nkmb+1,nz ; do I=Isq,Ieq
CS%uh_Rlay(I,j,k) = uh(I,j,k)
enddo ; enddo
k_list = nz/2
do k=1,nkmb ; do I=Isq,Ieq
call find_weights(G%Rlay, 0.5*(Rcv(i,j,k)+Rcv(i+1,j,k)), k_list, nz, wt, wt_p)
CS%uh_Rlay(I,j,k_list) = CS%uh_Rlay(I,j,k_list) + uh(I,j,k)*wt
CS%uh_Rlay(I,j,k_list+1) = CS%uh_Rlay(I,j,k_list+1) + uh(I,j,k)*wt_p
enddo ; enddo
enddo

if (CS%id_uh_Rlay > 0) call post_data(CS%id_uh_Rlay, CS%uh_Rlay, CS%diag)
endif

if (ASSOCIATED(CS%vh_Rlay)) then
do k=1,nkmb ; do J=Jsq,Jeq; do i=is,ie
CS%vh_Rlay(i,J,k) = 0.0
enddo ; enddo ; enddo
do k=nkmb+1,nz ; do J=Jsq,Jeq ; do i=is,ie
CS%vh_Rlay(i,J,k) = vh(i,J,k)
enddo ; enddo ; enddo
k_list = nz/2
do J=Jsq,Jeq ; do k=1,nkmb ; do i=is,ie
call find_weights(G%Rlay, 0.5*(Rcv(i,j,k)+Rcv(i,j+1,k)), k_list, nz, wt, wt_p)
CS%vh_Rlay(i,J,k_list) = CS%vh_Rlay(i,J,k_list) + vh(i,J,k)*wt
CS%vh_Rlay(i,J,k_list+1) = CS%vh_Rlay(i,J,k_list+1) + vh(i,J,k)*wt_p
enddo ; enddo ; enddo
!$OMP parallel do default(shared) private(wt,wt_p) firstprivate(k_list)
do J=Jsq,Jeq
do k=1,nkmb ; do i=is,ie
CS%vh_Rlay(i,J,k) = 0.0
enddo ; enddo
do k=nkmb+1,nz ; do i=is,ie
CS%vh_Rlay(i,J,k) = vh(i,J,k)
enddo ; enddo
do k=1,nkmb ; do i=is,ie
call find_weights(G%Rlay, 0.5*(Rcv(i,j,k)+Rcv(i,j+1,k)), k_list, nz, wt, wt_p)
CS%vh_Rlay(i,J,k_list) = CS%vh_Rlay(i,J,k_list) + vh(i,J,k)*wt
CS%vh_Rlay(i,J,k_list+1) = CS%vh_Rlay(i,J,k_list+1) + vh(i,J,k)*wt_p
enddo ; enddo
enddo

if (CS%id_vh_Rlay > 0) call post_data(CS%id_vh_Rlay, CS%vh_Rlay, CS%diag)
endif

if (ASSOCIATED(CS%uhGM_Rlay) .and. ASSOCIATED(CDp%uhGM)) then
do k=1,nkmb ; do j=js,je; do I=Isq,Ieq
CS%uhGM_Rlay(I,j,k) = 0.0
enddo ; enddo ; enddo
do k=nkmb+1,nz ; do j=js,je ; do I=Isq,Ieq
CS%uhGM_Rlay(I,j,k) = CDp%uhGM(I,j,k)
enddo ; enddo ; enddo
k_list = nz/2
do j=js,je ; do k=1,nkmb ; do I=Isq,Ieq
call find_weights(G%Rlay, 0.5*(Rcv(i,j,k)+Rcv(i+1,j,k)), k_list, nz, wt, wt_p)
CS%uhGM_Rlay(I,j,k_list) = CS%uhGM_Rlay(I,j,k_list) + CDp%uhGM(I,j,k)*wt
CS%uhGM_Rlay(I,j,k_list+1) = CS%uhGM_Rlay(I,j,k_list+1) + CDp%uhGM(I,j,k)*wt_p
enddo ; enddo ; enddo
!$OMP parallel do default(shared) private(wt,wt_p) firstprivate(k_list)
do j=js,je
do k=1,nkmb ; do I=Isq,Ieq
CS%uhGM_Rlay(I,j,k) = 0.0
enddo ; enddo
do k=nkmb+1,nz ; do I=Isq,Ieq
CS%uhGM_Rlay(I,j,k) = CDp%uhGM(I,j,k)
enddo ; enddo
do k=1,nkmb ; do I=Isq,Ieq
call find_weights(G%Rlay, 0.5*(Rcv(i,j,k)+Rcv(i+1,j,k)), k_list, nz, wt, wt_p)
CS%uhGM_Rlay(I,j,k_list) = CS%uhGM_Rlay(I,j,k_list) + CDp%uhGM(I,j,k)*wt
CS%uhGM_Rlay(I,j,k_list+1) = CS%uhGM_Rlay(I,j,k_list+1) + CDp%uhGM(I,j,k)*wt_p
enddo ; enddo
enddo

if (CS%id_uh_Rlay > 0) call post_data(CS%id_uhGM_Rlay, CS%uhGM_Rlay, CS%diag)
endif

if (ASSOCIATED(CS%vhGM_Rlay) .and. ASSOCIATED(CDp%vhGM)) then
do k=1,nkmb ; do J=Jsq,Jeq; do i=is,ie
CS%vhGM_Rlay(i,J,k) = 0.0
enddo ; enddo ; enddo
do k=nkmb+1,nz ; do J=Jsq,Jeq ; do i=is,ie
CS%vhGM_Rlay(i,J,k) = CDp%vhGM(i,J,k)
enddo ; enddo ; enddo
k_list = nz/2
do J=Jsq,Jeq ; do k=1,nkmb ; do i=is,ie
call find_weights(G%Rlay, 0.5*(Rcv(i,j,k)+Rcv(i,j+1,k)), k_list, nz, wt, wt_p)
CS%vhGM_Rlay(i,J,k_list) = CS%vhGM_Rlay(i,J,k_list) + CDp%vhGM(i,J,k)*wt
CS%vhGM_Rlay(i,J,k_list+1) = CS%vhGM_Rlay(i,J,k_list+1) + CDp%vhGM(i,J,k)*wt_p
enddo ; enddo ; enddo
!$OMP parallel do default(shared) private(wt,wt_p) firstprivate(k_list)
do J=Jsq,Jeq
do k=1,nkmb ; do i=is,ie
CS%vhGM_Rlay(i,J,k) = 0.0
enddo ; enddo
do k=nkmb+1,nz ; do i=is,ie
CS%vhGM_Rlay(i,J,k) = CDp%vhGM(i,J,k)
enddo ; enddo
do k=1,nkmb ; do i=is,ie
call find_weights(G%Rlay, 0.5*(Rcv(i,j,k)+Rcv(i,j+1,k)), k_list, nz, wt, wt_p)
CS%vhGM_Rlay(i,J,k_list) = CS%vhGM_Rlay(i,J,k_list) + CDp%vhGM(i,J,k)*wt
CS%vhGM_Rlay(i,J,k_list+1) = CS%vhGM_Rlay(i,J,k_list+1) + CDp%vhGM(i,J,k)*wt_p
enddo ; enddo
enddo

if (CS%id_vhGM_Rlay > 0) call post_data(CS%id_vhGM_Rlay, CS%vhGM_Rlay, CS%diag)
endif
endif
Expand All @@ -316,6 +338,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, dt, G, &
call wave_speed(h, tv, G, CS%cg1, CS%wave_speed_CSp)
if (CS%id_cg1>0) call post_data(CS%id_cg1, CS%cg1, CS%diag)
if (CS%id_Rd1>0) then
!$OMP parallel do default(shared) private(f2_h,mag_beta)
do j=js,je ; do i=is,ie
! Blend the equatorial deformation radius with the standard one.
f2_h = absurdly_small_freq2 + 0.25 * &
Expand Down Expand Up @@ -424,7 +447,7 @@ subroutine calculate_vertical_integrals(h, tv, G, CS)
real, dimension(NIMEM_,NJMEM_,NKMEM_), intent(in) :: h
type(thermo_var_ptrs), intent(in) :: tv
type(ocean_grid_type), intent(inout) :: G
type(diagnostics_CS), pointer :: CS
type(diagnostics_CS), intent(inout) :: CS
! This subroutine calculates the vertical integrals of several tracers, along
! with the mass-weight of these tracers, the total column mass, and the
! carefully calculated column height.
Expand Down Expand Up @@ -523,7 +546,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, CS)
type(accel_diag_ptrs), intent(in) :: ADp
type(cont_diag_ptrs), intent(in) :: CDp
type(ocean_grid_type), intent(inout) :: G
type(diagnostics_CS), pointer :: CS
type(diagnostics_CS), intent(inout) :: CS
! This subroutine calculates a series of terms in the energy
! balance equations.

Expand Down Expand Up @@ -739,7 +762,7 @@ end subroutine register_time_deriv
subroutine calculate_derivs(dt, G, CS)
real, intent(in) :: dt
type(ocean_grid_type), intent(inout) :: G
type(diagnostics_CS), pointer :: CS
type(diagnostics_CS), intent(inout) :: CS
! This subroutine calculates all registered time derivatives.
! Arguments: dt - the time interval in s over which differences occur
! (in) G - The ocean's grid structure.
Expand Down
2 changes: 1 addition & 1 deletion src/parameterizations/lateral/MOM_hor_visc.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -296,7 +296,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, CS, OBC)
"VarMix%Res_fn_q both need to be associated with Resoln_scaled_Kh.")
endif

!$OMP parallel do default(shared) private(i, j, k, u0, v0, sh_xx, str_xx, &
!$OMP parallel do default(shared) private(i, j, k, u0, v0, sh_xx, str_xx, visc_bound_rem,&
!$OMP sh_xy, str_xy, Ah, Kh, AhSm, KhSm, &
!$OMP Shear_mag, huq, hvq, hq, Kh_scale, hrat_min)
do k=1,nz
Expand Down
63 changes: 41 additions & 22 deletions src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -271,6 +271,8 @@ subroutine calc_slope_function_(h, tv, G, CS, e)
real :: one_meter ! One meter in thickness units of m or kg m-2.
integer :: is, ie, js, je, nz
integer :: i, j, k, kb_max
real :: SN_u_local(SZIB_(G), SZJ_(G),SZK_(G))
real :: SN_v_local(SZI_(G), SZJB_(G),SZK_(G))

if (.not. ASSOCIATED(CS)) call MOM_error(FATAL, "calc_slope_function:"// &
"Module must be initialized before it is used.")
Expand All @@ -286,6 +288,8 @@ subroutine calc_slope_function_(h, tv, G, CS, e)
h_neglect = G%H_subroundoff
H_cutoff = real(2*nz) * (G%Angstrom + h_neglect)

!$OMP parallel default(shared) private(E_x,E_y,S2,Hdn,Hup,H_geom,N2)
!$OMP do
do j=js-1,je+1 ; do i=is-1,ie+1
CS%SN_u(i,j) = 0.0
CS%SN_v(i,j) = 0.0
Expand All @@ -296,14 +300,16 @@ subroutine calc_slope_function_(h, tv, G, CS, e)
! and midlatitude deformation radii, using calc_resoln_function as a template.

! Set the length scale at u-points.
!$OMP do
do j=js,je ; do I=is-1,ie
CS%L2u(I,j) = CS%Visbeck_L_scale**2
enddo ; enddo
! Set length scale at v-points
!$OMP do
do J=js-1,je ; do i=is,ie
CS%L2v(i,J) = CS%Visbeck_L_scale**2
enddo ; enddo

!$OMP do
do k=nz,CS%VarMix_Ktop,-1

! Calculate the interface slopes E_x and E_y and u- and v- points respectively
Expand All @@ -328,7 +334,7 @@ subroutine calc_slope_function_(h, tv, G, CS, e)
N2 = G%g_prime(k) / (G%H_to_m * max(Hdn,Hup,one_meter))
if (min(h(i,j,k-1), h(i+1,j,k-1), h(i,j,k), h(i+1,j,k)) < H_cutoff) &
S2 = 0.0
CS%SN_u(I,j) = CS%SN_u(I,j) + (H_geom * G%H_to_m) * S2 * N2
SN_u_local(I,j,k) = (H_geom * G%H_to_m) * S2 * N2
enddo ; enddo
do J=js-1,je ; do i=is,ie
S2 = ( E_y(i,J)**2 + 0.25*( &
Expand All @@ -339,30 +345,43 @@ subroutine calc_slope_function_(h, tv, G, CS, e)
N2 = G%g_prime(k) / (G%H_to_m * max(Hdn,Hup,one_meter))
if (min(h(i,j,k-1), h(i,j+1,k-1), h(i,j,k), h(i,j+1,k)) < H_cutoff) &
S2 = 0.0
CS%SN_v(i,J) = CS%SN_v(i,J) + (H_geom * G%H_to_m) * S2 * N2
SN_v_local(i,J,k) = (H_geom * G%H_to_m) * S2 * N2
enddo ; enddo

enddo ! k
!$OMP do
do j = js,je;
do k=nz,CS%VarMix_Ktop,-1 ; do I=is-1,ie
CS%SN_u(I,j) = CS%SN_u(I,j) + SN_u_local(I,j,k)
enddo ; enddo
! SN above contains S^2*N^2*H, convert to vertical average of S*N
do I=is-1,ie
!SN_u(I,j) = sqrt( SN_u(I,j) / ( max(G%bathyT(I,j), G%bathyT(I+1,j)) + G%Angstrom ) )
!The code below behaves better than the line above. Not sure why? AJA
if ( min(G%bathyT(I,j), G%bathyT(I+1,j)) > H_cutoff ) then
CS%SN_u(I,j) = sqrt( CS%SN_u(I,j) / max(G%bathyT(I,j), G%bathyT(I+1,j)) )
else
CS%SN_u(I,j) = 0.0
endif
enddo
enddo
!$OMP do
do J=js-1,je
do k=nz,CS%VarMix_Ktop,-1 ; do I=is,ie
CS%SN_v(i,J) = CS%SN_v(i,J) + SN_v_local(i,J,k)
enddo ; enddo
do i=is,ie
!SN_v(i,J) = sqrt( SN_v(i,J) / ( max(G%bathyT(i,J), G%bathyT(i,J+1)) + G%Angstrom ) )
!The code below behaves better than the line above. Not sure why? AJA
if ( min(G%bathyT(I,j), G%bathyT(I+1,j)) > H_cutoff ) then
CS%SN_v(i,J) = sqrt( CS%SN_v(i,J) / max(G%bathyT(i,J), G%bathyT(i,J+1)) )
else
CS%SN_v(I,j) = 0.0
endif
enddo
enddo

! SN above contains S^2*N^2*H, convert to vertical average of S*N
do j=js,je ; do I=is-1,ie
!SN_u(I,j) = sqrt( SN_u(I,j) / ( max(G%bathyT(I,j), G%bathyT(I+1,j)) + G%Angstrom ) )
!The code below behaves better than the line above. Not sure why? AJA
if ( min(G%bathyT(I,j), G%bathyT(I+1,j)) > H_cutoff ) then
CS%SN_u(I,j) = sqrt( CS%SN_u(I,j) / max(G%bathyT(I,j), G%bathyT(I+1,j)) )
else
CS%SN_u(I,j) = 0.0
endif
enddo ; enddo
do J=js-1,je ; do i=is,ie
!SN_v(i,J) = sqrt( SN_v(i,J) / ( max(G%bathyT(i,J), G%bathyT(i,J+1)) + G%Angstrom ) )
!The code below behaves better than the line above. Not sure why? AJA
if ( min(G%bathyT(I,j), G%bathyT(I+1,j)) > H_cutoff ) then
CS%SN_v(i,J) = sqrt( CS%SN_v(i,J) / max(G%bathyT(i,J), G%bathyT(i,J+1)) )
else
CS%SN_v(I,j) = 0.0
endif
enddo ; enddo
!$OMP end parallel

! Offer diagnostic fields for averaging.
if (query_averaging_enabled(CS%diag)) then
Expand Down
Loading