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

Fix for openmp #16

Merged
merged 3 commits into from
May 1, 2014
Merged

Fix for openmp #16

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
a175db8
Improve Openmp performance and single thread performance
Zhi-Liang Apr 29, 2014
abae894
openmp implementation
Zhi-Liang Apr 29, 2014
a94851e
Merge branch 'user/z1l/openmp' of github.com:CommerceGov/NOAA-GFDL-MO…
Zhi-Liang Apr 29, 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
3 changes: 0 additions & 3 deletions examples/solo_ocean/double_gyre/input.nml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,4 @@
parameter_filename = 'MOM_input',
'MOM_override' /

&fms_nml
domains_stack_size = 710000,
stack_size = 0 /

38 changes: 21 additions & 17 deletions src/ALE/MOM_ALE.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -320,6 +320,7 @@ subroutine pressure_gradient_plm( CS, S_t, S_b, T_t, T_b, G, tv, h )
! Local variables
integer :: i, j, k
real :: hTmp(G%ke)
real :: tmp(G%ke)
real, dimension(CS%nk,2) :: &
ppoly_linear_E !Edge value of polynomial
real, dimension(CS%nk,CS%degree_linear+1) :: &
Expand All @@ -331,19 +332,18 @@ subroutine pressure_gradient_plm( CS, S_t, S_b, T_t, T_b, G, tv, h )
! in 'ALE_memory_allocation'.

! Determine reconstruction within each column
!$OMP parallel do default(shared) private(hTmp,ppoly_linear_E,ppoly_linear_coefficients)
do i = G%isc,G%iec+1
do j = G%jsc,G%jec+1

!$OMP parallel do default(shared) private(hTmp,ppoly_linear_E,ppoly_linear_coefficients,tmp)
do j = G%jsc,G%jec+1
do i = G%isc,G%iec+1
! Build current grid
hTmp(:) = h(i,j,:)

tmp(:) = tv%S(i,j,:)
! Reconstruct salinity profile
ppoly_linear_E = 0.0
ppoly_linear_coefficients = 0.0
call PLM_reconstruction( G%ke, hTmp, tv%S(i,j,:), ppoly_linear_E, ppoly_linear_coefficients )
call PLM_reconstruction( G%ke, hTmp, tmp, ppoly_linear_E, ppoly_linear_coefficients )
if (CS%boundary_extrapolation_for_pressure) call &
PLM_boundary_extrapolation( G%ke, hTmp, tv%S(i,j,:), ppoly_linear_E, ppoly_linear_coefficients )
PLM_boundary_extrapolation( G%ke, hTmp, tmp, ppoly_linear_E, ppoly_linear_coefficients )

do k = 1,G%ke
S_t(i,j,k) = ppoly_linear_E(k,1)
Expand All @@ -353,9 +353,10 @@ subroutine pressure_gradient_plm( CS, S_t, S_b, T_t, T_b, G, tv, h )
! Reconstruct temperature profile
ppoly_linear_E = 0.0
ppoly_linear_coefficients = 0.0
call PLM_reconstruction( G%ke, hTmp, tv%T(i,j,:), ppoly_linear_E, ppoly_linear_coefficients )
tmp(:) = tv%T(i,j,:)
call PLM_reconstruction( G%ke, hTmp, tmp, ppoly_linear_E, ppoly_linear_coefficients )
if (CS%boundary_extrapolation_for_pressure) call &
PLM_boundary_extrapolation( G%ke, hTmp, tv%T(i,j,:), ppoly_linear_E, ppoly_linear_coefficients )
PLM_boundary_extrapolation( G%ke, hTmp, tmp, ppoly_linear_E, ppoly_linear_coefficients )

do k = 1,G%ke
T_t(i,j,k) = ppoly_linear_E(k,1)
Expand Down Expand Up @@ -393,6 +394,7 @@ subroutine pressure_gradient_ppm( CS, S_t, S_b, T_t, T_b, G, tv, h )
! Local variables
integer :: i, j, k
real :: hTmp(G%ke)
real :: tmp(G%ke)
real, dimension(CS%nk,2) :: &
ppoly_parab_E !Edge value of polynomial
real, dimension(CS%nk,CS%degree_parab+1) :: &
Expand All @@ -405,19 +407,20 @@ subroutine pressure_gradient_ppm( CS, S_t, S_b, T_t, T_b, G, tv, h )

! Determine reconstruction within each column
!$OMP parallel do default(shared) private(hTmp,ppoly_parab_E,ppoly_parab_coefficients)
do i = G%isc,G%iec+1
do j = G%jsc,G%jec+1
do j = G%jsc,G%jec+1
do i = G%isc,G%iec+1

! Build current grid
hTmp(:) = h(i,j,:)
tmp(:) = tv%S(i,j,:)

! Reconstruct salinity profile
ppoly_parab_E = 0.0
ppoly_parab_coefficients = 0.0
call edge_values_implicit_h4( G%ke, hTmp, tv%S(i,j,:), CS%edgeValueWrk, ppoly_parab_E )
call PPM_reconstruction( G%ke, hTmp, tv%S(i,j,:), ppoly_parab_E, ppoly_parab_coefficients )
call edge_values_implicit_h4( G%ke, hTmp, tmp, CS%edgeValueWrk, ppoly_parab_E )
call PPM_reconstruction( G%ke, hTmp, tmp, ppoly_parab_E, ppoly_parab_coefficients )
if (CS%boundary_extrapolation_for_pressure) call &
PPM_boundary_extrapolation( G%ke, hTmp, tv%S(i,j,:), ppoly_parab_E, ppoly_parab_coefficients )
PPM_boundary_extrapolation( G%ke, hTmp, tmp, ppoly_parab_E, ppoly_parab_coefficients )

do k = 1,G%ke
S_t(i,j,k) = ppoly_parab_E(k,1)
Expand All @@ -427,10 +430,11 @@ subroutine pressure_gradient_ppm( CS, S_t, S_b, T_t, T_b, G, tv, h )
! Reconstruct temperature profile
ppoly_parab_E = 0.0
ppoly_parab_coefficients = 0.0
call edge_values_implicit_h4( G%ke, hTmp, tv%T(i,j,:), CS%edgeValueWrk, ppoly_parab_E )
call PPM_reconstruction( G%ke, hTmp, tv%T(i,j,:), ppoly_parab_E, ppoly_parab_coefficients )
tmp(:) = tv%T(i,j,:)
call edge_values_implicit_h4( G%ke, hTmp, tmp, CS%edgeValueWrk, ppoly_parab_E )
call PPM_reconstruction( G%ke, hTmp, tmp, ppoly_parab_E, ppoly_parab_coefficients )
if (CS%boundary_extrapolation_for_pressure) call &
PPM_boundary_extrapolation( G%ke, hTmp, tv%T(i,j,:), ppoly_parab_E, ppoly_parab_coefficients )
PPM_boundary_extrapolation( G%ke, hTmp, tmp, ppoly_parab_E, ppoly_parab_coefficients )

do k = 1,G%ke
T_t(i,j,k) = ppoly_parab_E(k,1)
Expand Down
38 changes: 25 additions & 13 deletions src/core/MOM.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1257,7 +1257,7 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS)

if (showCallTree) call callTree_leave("DT cycles (step_MOM)")

call cpu_clock_end(id_clock_other)
call cpu_clock_end(id_clock_other)

enddo ! End of n loop

Expand Down Expand Up @@ -2375,6 +2375,7 @@ subroutine calculate_surface_state(state, u, v, h, ssh, G, CS, p_atm)

depth_ml = CS%Hmix
! Determine the mean properties of the uppermost depth_ml fluid.
!$OMP parallel do default(shared) private(depth,dh)
do j=js,je
do i=is,ie
depth(i) = 0.0
Expand Down Expand Up @@ -2410,7 +2411,7 @@ subroutine calculate_surface_state(state, u, v, h, ssh, G, CS, p_atm)
else
state%sfc_density(i,j) = state%sfc_density(i,j) / depth(i)
endif
state%Hml(:,:) = depth(i)
state%Hml(i,j) = depth(i)
enddo
enddo ! end of j loop
endif ! end BULKMIXEDLAYER
Expand All @@ -2421,13 +2422,6 @@ subroutine calculate_surface_state(state, u, v, h, ssh, G, CS, p_atm)
state%TempxPmE => CS%tv%TempxPmE
state%internal_heat => CS%tv%internal_heat

if (associated(state%salt_deficit) .and. associated(CS%tv%salt_deficit)) then
do j=js,je ; do i=is,ie
! Convert from gSalt to kgSalt
state%salt_deficit(i,j) = 1000.0 * CS%tv%salt_deficit(i,j)
enddo ; enddo
endif

! Allocate structures for ocean_mass, ocean_heat, and ocean_salt. This could
! be wrapped in a run-time flag to disable it for economy, since the 3-d
! sums are not negligible.
Expand All @@ -2443,13 +2437,24 @@ subroutine calculate_surface_state(state, u, v, h, ssh, G, CS, p_atm)
endif
endif

!$OMP parallel default(shared) private(mass)
if (associated(state%salt_deficit) .and. associated(CS%tv%salt_deficit)) then
!$OMP do
do j=js,je ; do i=is,ie
! Convert from gSalt to kgSalt
state%salt_deficit(i,j) = 1000.0 * CS%tv%salt_deficit(i,j)
enddo ; enddo
endif

if (associated(state%ocean_mass) .and. associated(state%ocean_heat) .and. &
associated(state%ocean_salt)) then
!$OMP do
do j=js,je ; do i=is,ie
state%ocean_mass(i,j) = 0.0
state%ocean_heat(i,j) = 0.0 ; state%ocean_salt(i,j) = 0.0
enddo ; enddo
do k=1,nz ; do j=js,je ; do i=is,ie
!$OMP do
do j=js,je ; do k=1,nz; do i=is,ie
mass = G%H_to_kg_m2*h(i,j,k)
state%ocean_mass(i,j) = state%ocean_mass(i,j) + mass
state%ocean_heat(i,j) = state%ocean_heat(i,j) + mass*CS%tv%T(i,j,k)
Expand All @@ -2458,27 +2463,34 @@ subroutine calculate_surface_state(state, u, v, h, ssh, G, CS, p_atm)
enddo ; enddo ; enddo
else
if (associated(state%ocean_mass)) then
!$OMP do
do j=js,je ; do i=is,ie ; state%ocean_mass(i,j) = 0.0 ; enddo ; enddo
do k=1,nz ; do j=js,je ; do i=is,ie
!$OMP do
do j=js,je ; do k=1,nz ; do i=is,ie
state%ocean_mass(i,j) = state%ocean_mass(i,j) + G%H_to_kg_m2*h(i,j,k)
enddo ; enddo ; enddo
endif
if (associated(state%ocean_heat)) then
!$OMP do
do j=js,je ; do i=is,ie ; state%ocean_heat(i,j) = 0.0 ; enddo ; enddo
do k=1,nz ; do j=js,je ; do i=is,ie
!$OMP do
do j=js,je ; do k=1,nz ; do i=is,ie
mass = G%H_to_kg_m2*h(i,j,k)
state%ocean_heat(i,j) = state%ocean_heat(i,j) + mass*CS%tv%T(i,j,k)
enddo ; enddo ; enddo
endif
if (associated(state%ocean_salt)) then
!$OMP do
do j=js,je ; do i=is,ie ; state%ocean_salt(i,j) = 0.0 ; enddo ; enddo
do k=1,nz ; do j=js,je ; do i=is,ie
!$OMP do
do j=js,je ; do k=1,nz ; do i=is,ie
mass = G%H_to_kg_m2*h(i,j,k)
state%ocean_salt(i,j) = state%ocean_salt(i,j) + &
mass * (1.0e-3*CS%tv%S(i,j,k))
enddo ; enddo ; enddo
endif
endif
!$OMP end parallel

if (associated(CS%visc%taux_shelf)) state%taux_shelf => CS%visc%taux_shelf
if (associated(CS%visc%tauy_shelf)) state%tauy_shelf => CS%visc%tauy_shelf
Expand Down
3 changes: 3 additions & 0 deletions src/core/MOM_PressureForce_Montgomery.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -665,6 +665,7 @@ subroutine Set_pbce_Bouss(e, tv, G, g_Earth, Rho0, GFS_scale, pbce, rho_star)

if (use_EOS) then
if (present(rho_star)) then
!$OMP parallel do default(shared) private(Ihtot)
do j=Jsq,Jeq+1
do i=Isq,Ieq+1
Ihtot(i) = 1.0 / ((e(i,j,1)-e(i,j,nz+1)) + h_neglect)
Expand All @@ -676,6 +677,7 @@ subroutine Set_pbce_Bouss(e, tv, G, g_Earth, Rho0, GFS_scale, pbce, rho_star)
enddo ; enddo
enddo ! end of j loop
else
!$OMP parallel do default(shared) private(Ihtot,press,rho_in_situ,T_int,S_int,dR_dT,dR_dS)
do j=Jsq,Jeq+1
do i=Isq,Ieq+1
Ihtot(i) = 1.0 / ((e(i,j,1)-e(i,j,nz+1)) + h_neglect)
Expand Down Expand Up @@ -704,6 +706,7 @@ subroutine Set_pbce_Bouss(e, tv, G, g_Earth, Rho0, GFS_scale, pbce, rho_star)
enddo ! end of j loop
endif
else ! not use_EOS
!$OMP parallel do default(shared) private(Ihtot)
do j=Jsq,Jeq+1
do i=Isq,Ieq+1
Ihtot(i) = 1.0 / ((e(i,j,1)-e(i,j,nz+1)) + h_neglect)
Expand Down
1 change: 1 addition & 0 deletions src/core/MOM_forcing_type.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -461,6 +461,7 @@ subroutine calculateBuoyancyFlux2d(G, fluxes, optics, h, Temp, Salt, tv, buoyanc

netT(:) = 0. ; netS(:) = 0.

!$OMP parallel do default(shared) firstprivate(netT,netS)
do j = G%jsc, G%jec
call calculateBuoyancyFlux1d(G, fluxes, optics, h, Temp, Salt, tv, j, buoyancyFlux(:,j,:), netT, netS )
if (present(netHeatMinusSW)) netHeatMinusSW(:,j) = netT(:)
Expand Down
3 changes: 2 additions & 1 deletion src/diagnostics/MOM_wave_speed.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -163,7 +163,8 @@ subroutine wave_speed(h, tv, G, cg1, CS, full_halos)
rescale = 1024.0**4 ; I_rescale = 1.0/rescale

min_h_frac = tol1 / real(nz)

!$OMP parallel do default(private) shared(is,ie,js,je,nz,h,G,min_h_frac,use_EOS,T,S, &
!$OMP H_to_pres,tv,cg1,g_Rho0,rescale,I_rescale)
do j=js,je
! First merge very thin layers with the one above (or below if they are
! at the top). This also transposes the row order so that columns can
Expand Down
6 changes: 5 additions & 1 deletion src/parameterizations/vertical/MOM_KPP.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -380,6 +380,10 @@ subroutine KPP_calculate(CS, G, h, Temp, Salt, u, v, EOS, uStar, buoyFlux, Kt, K

if (CS%id_Kd_in > 0) call post_data(CS%id_Kd_in, Kt, CS%diag)

!$OMP parallel do default(private) shared(G,CS,EOS,uStar,Temp,Salt,u,v,h,GoRho,buoyFlux, &
!$OMP const1,nonLocalTransHeat, &
!$OMP nonLocalTransScalar,Kt,Ks,Kv) &
!$OMP firstprivate(nonLocalTrans)
do j = G%jsc, G%jec
do i = G%isc, G%iec
if (G%mask2dT(i,j)==0.) cycle ! Skip calling KPP for land points
Expand Down Expand Up @@ -826,7 +830,7 @@ subroutine KPP_applyNonLocalTransport(CS, G, h, nonLocalTrans, surfFlux, dt, sca
endif

if (diagHeat .or. diagSalt) dSdt(:,:,:) = 0. ! Zero halos for diagnostics ???

!$OMP parallel do default(shared)
do k = 1, G%ke
do j = G%jsc, G%jec
do i = G%isc, G%iec
Expand Down
26 changes: 16 additions & 10 deletions src/parameterizations/vertical/MOM_diabatic_driver.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -164,6 +164,7 @@ module MOM_diabatic_driver
type(diag_ctrl), pointer :: diag ! A structure that is used to regulate the
! timing of diagnostic output.
real :: MLDdensityDifference ! Density difference used to determine MLD_user
integer :: nsw ! SW_NBANDS
integer :: id_dudt_dia = -1, id_dvdt_dia = -1, id_wd = -1
integer :: id_ea = -1 , id_eb = -1, id_Kd_z = -1, id_Kd_interface = -1
integer :: id_Tdif_z = -1, id_Tadv_z = -1, id_Sdif_z = -1, id_Sadv_z = -1
Expand Down Expand Up @@ -1413,6 +1414,7 @@ subroutine triDiagTS(G, is, ie, js, je, hold, ea, eb, T, S)
real :: h_tr, b_denom_1
integer :: i, j, k

!$OMP parallel do default(shared) private(h_tr,b1,d1,c1,b_denom_1)
do j=js,je
do i=is,ie
h_tr = hold(i,j,1) + G%H_subroundoff
Expand Down Expand Up @@ -1733,6 +1735,9 @@ subroutine diabatic_driver_init(Time, G, param_file, useALEalgorithm, diag, &
endif
endif

CS%nsw = 0
if (ASSOCIATED(CS%optics)) CS%nsw = CS%optics%nbands

end subroutine diabatic_driver_init

subroutine diabatic_driver_end(CS)
Expand Down Expand Up @@ -1915,8 +1920,8 @@ subroutine applyBoundaryFluxes(CS, G, dt, fluxes, optics, ea, h, tv)
real, dimension(SZI_(G)) :: netThickness, netHeat, netSalt, htot, Ttot
real, dimension(SZI_(G), SZK_(G)) :: h2d, T2d, eps
integer, dimension(SZI_(G), SZK_(G)) :: ksort
real, allocatable, dimension(:,:) :: Pen_SW_bnd
real, allocatable, dimension(:,:,:) :: opacityBand
real, dimension(max(CS%nsw,1),SZI_(G)) :: Pen_SW_bnd
real, dimension(max(CS%nsw,1),SZI_(G),SZK_(G)) :: opacityBand
logical :: use_riverHeatContent, useCalvingHeatContent
integer, parameter :: maxGroundings = 5
integer :: numberOfGroundings, iGround(maxGroundings), jGround(maxGroundings)
Expand All @@ -1928,10 +1933,7 @@ subroutine applyBoundaryFluxes(CS, G, dt, fluxes, optics, ea, h, tv)
! Skip applying forcing if fluxes%sw is not associated.
if (.not.ASSOCIATED(fluxes%sw)) return

nsw = 0
if (ASSOCIATED(optics)) nsw = optics%nbands
allocate(Pen_SW_bnd(max(nsw,1),SZI_(G)))
allocate(opacityBand(max(nsw,1),SZI_(G),SZK_(G)))
nsw = CS%nsw

Irho0 = 1.0 / G%Rho0
I_Cp = 1.0 / fluxes%C_p
Expand All @@ -1957,7 +1959,11 @@ subroutine applyBoundaryFluxes(CS, G, dt, fluxes, optics, ea, h, tv)
if (CS%id_createdH>0) CS%createdH(:,:) = 0.

numberOfGroundings = 0
!$OMP parallel do default(shared) private(opacityBand,h2d,T2d,eps,htot,netThickness, &
!$OMP netHeat,netSalt,Pen_SW_bnd,fractionOfForcing, &
!$OMP dThickness,dTemp,dSalt,hOld,Ithickness,Ttot)
do j=js,je ! Work in vertical slices (this is a hold over from the routines called with a j argument)
Ttot = 0
! Copy state into 2D-slice arrays
do k=1,nz
do i=is,ie
Expand Down Expand Up @@ -2026,11 +2032,13 @@ subroutine applyBoundaryFluxes(CS, G, dt, fluxes, optics, ea, h, tv)
endif
enddo ! k
if (netThickness(i)/=0.) then ! If anything remains then we grounded out
!$OMP critical
numberOfGroundings = numberOfGroundings +1
!$OMP end critical
if (numberOfGroundings<=maxGroundings) then
iGround(numberOfGroundings) = i ! Record i,j location of event for
jGround(numberOfGroundings) = j ! warning message
hGrounding = netThickness(i)
hGrounding(numberOfGroundings) = netThickness(i)
endif
if (CS%id_createdH>0) CS%createdH(i,j) = CS%createdH(i,j) - netThickness(i)/dt
endif
Expand All @@ -2057,6 +2065,7 @@ subroutine applyBoundaryFluxes(CS, G, dt, fluxes, optics, ea, h, tv)
enddo

enddo ! j

if (CS%id_createdH > 0) call post_data(CS%id_createdH, CS%createdH, CS%diag)

if (numberOfGroundings>0) then
Expand All @@ -2069,9 +2078,6 @@ subroutine applyBoundaryFluxes(CS, G, dt, fluxes, optics, ea, h, tv)
enddo
endif

deallocate(Pen_SW_bnd)
deallocate(opacityBand)

end subroutine applyBoundaryFluxes

end module MOM_diabatic_driver
Loading