Skip to content

Fix for openmp #5

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

Merged
adcroft merged 3 commits into from
Mar 13, 2014
Merged

Fix for openmp #5

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
5d52916
Fix for openmp
Zhi-Liang Mar 10, 2014
7dc7181
Implement openmp
Zhi-Liang Mar 10, 2014
be17276
remove commented-out codes
Zhi-Liang Mar 10, 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
46 changes: 23 additions & 23 deletions src/parameterizations/vertical/MOM_kappa_shear.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -321,6 +321,9 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, &
k0dt = dt*CS%kappa_0
dz_massless = 0.1*sqrt(k0dt)

!$OMP parallel do default(private) shared(js,je,is,ie,nz,h,u_in,v_in,use_temperature,new_kappa, &
!$OMP tv,G,CS,kappa_io,dz_massless,k0dt,p_surf,gR0,g_R0,dt, &
!$OMP tol_dksrc,tol_dksrc_low,tol2,Ri_crit,dt_refinements)
do j=js,je
do k=1,nz ; do i=is,ie
h_2d(i,k) = h(i,j,k)*G%H_to_m
Expand Down Expand Up @@ -877,8 +880,7 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, &

return

contains
! end subroutine Calculate_kappa_shear
end subroutine Calculate_kappa_shear

subroutine calculate_projected_state(kappa, u0, v0, T0, S0, dt, nz, &
dz, I_dz_int, dbuoy_dT, dbuoy_dS, &
Expand Down Expand Up @@ -919,8 +921,7 @@ subroutine calculate_projected_state(kappa, u0, v0, T0, S0, dt, nz, &
! (in,opt) ke_int - The bottommost k-index with a non-zero diffusivity.

! UNCOMMENT THE FOLLOWING IF NOT CONTAINED IN THE OUTER SUBROUTINE.
! real, dimension(nz+1) :: &
! c1
real, dimension(nz+1) :: c1
real :: a_a, a_b, b1, d1, bd1, b1nz_0
integer :: k, ks, ke

Expand Down Expand Up @@ -1040,24 +1041,24 @@ subroutine find_kappa_tke(N2, S2, kappa_in, Idz, dz_Int, I_L2_bdry, f2, &
! (out,opt) local_src - The sum of all local sources for kappa, in s-1.

! UNCOMMENT THE FOLLOWING IF NOT CONTAINED IN Calculate_kappa_shear
! real, dimension(nz) :: &
! aQ, & ! aQ is the coupling between adjacent interfaces in the TKE
! ! equations, in m s-1.
! dQdz ! Half the partial derivative of TKE with depth, m s-2.
! real, dimension(nz+1) :: &
! dK, & ! The change in kappa, in m2 s-1.
! dQ, & ! The change in TKE, in m2 s-1.
! cQ, cK, & ! cQ and cK are the upward influences in the tridiagonal and
! ! hexadiagonal solvers for the TKE and kappa equations, ND.
! I_Ld2, & ! 1/Ld^2, where Ld is the effective decay length scale
! ! for kappa, in units of m-2.
! TKE_decay, & ! The local TKE decay rate in s-1.
! k_src, & ! The source term in the kappa equation, in s-1.
! dQmdK, & ! With Newton's method the change in dQ(k-1) due to dK(k), s.
! dKdQ, & ! With Newton's method the change in dK(k) due to dQ(k), s-1.
! e1 ! The fractional change in a layer TKE due to a change in the
! ! TKE of the layer above when all the kappas below are 0.
! ! e1 is nondimensional, and 0 < e1 < 1.
real, dimension(nz) :: &
aQ, & ! aQ is the coupling between adjacent interfaces in the TKE
! equations, in m s-1.
dQdz ! Half the partial derivative of TKE with depth, m s-2.
real, dimension(nz+1) :: &
dK, & ! The change in kappa, in m2 s-1.
dQ, & ! The change in TKE, in m2 s-1.
cQ, cK, & ! cQ and cK are the upward influences in the tridiagonal and
! hexadiagonal solvers for the TKE and kappa equations, ND.
I_Ld2, & ! 1/Ld^2, where Ld is the effective decay length scale
! for kappa, in units of m-2.
TKE_decay, & ! The local TKE decay rate in s-1.
k_src, & ! The source term in the kappa equation, in s-1.
dQmdK, & ! With Newton's method the change in dQ(k-1) due to dK(k), s.
dKdQ, & ! With Newton's method the change in dK(k) due to dQ(k), s-1.
e1 ! The fractional change in a layer TKE due to a change in the
! TKE of the layer above when all the kappas below are 0.
! e1 is nondimensional, and 0 < e1 < 1.
real :: tke_src ! The net source of TKE due to mixing against the shear
! and stratification, in m2 s-3. (For convenience,
! a term involving the non-dissipation of q0 is also
Expand Down Expand Up @@ -1619,7 +1620,6 @@ subroutine find_kappa_tke(N2, S2, kappa_in, Idz, dz_Int, I_L2_bdry, f2, &

end subroutine find_kappa_tke

end subroutine Calculate_kappa_shear

logical function kappa_shear_init(Time, G, param_file, diag, CS)
type(time_type), intent(in) :: Time
Expand Down
85 changes: 52 additions & 33 deletions src/parameterizations/vertical/MOM_vert_friction.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -693,7 +693,7 @@ subroutine vertvisc_coef(u, v, h, fluxes, visc, dt, G, CS)

call find_coupling_coef(a, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, h_ml, dt, j, G, CS, visc, fluxes, work_on_u=.true.)
if (allocated(hML_u)) then
do i=is,ie ; if (do_i(i)) then ; hML_u(I,j) = h_ml(I) ; endif ; enddo
do i=isq,ieq ; if (do_i(i)) then ; hML_u(I,j) = h_ml(I) ; endif ; enddo
endif
do_any_shelf = .false.
if (associated(fluxes%frac_shelf_u)) then
Expand Down Expand Up @@ -1122,8 +1122,10 @@ subroutine vertvisc_limit_vel(u, v, h, ADp, CDp, fluxes, visc, dt, G, CS)
real :: maxvel ! Velocities components greater than maxvel
real :: truncvel ! are truncated to truncvel, both in m s-1.
real :: CFL ! The local CFL number.
real :: vel_report(SZIB_(G))
logical :: trunc_any, dowrite(SZIB_(G))
real :: vel_report(SZIB_(G),SZJB_(G))
real :: u_old(SZIB_(G),SZJ_(G),SZK_(G))
real :: v_old(SZI_(G),SZJB_(G),SZK_(G))
logical :: trunc_any, dowrite(SZIB_(G),SZJB_(G))
integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
Expand All @@ -1132,13 +1134,13 @@ subroutine vertvisc_limit_vel(u, v, h, ADp, CDp, fluxes, visc, dt, G, CS)
truncvel = 0.9*maxvel

if (len_trim(CS%u_trunc_file) > 0) then
if (.not.CS%CFL_based_trunc) vel_report(:) = maxvel
!$OMP parallel do default(private) shared(js,je,Isq,Ieq,nz,CS,G,fluxes,u,h)
!$OMP parallel do default(private) shared(js,je,Isq,Ieq,nz,CS,G,fluxes,u,h,dt,maxvel,ADp,CDp,truncvel, &
!$OMP u_old,vel_report,dowrite)
do j=js,je
trunc_any = .false.
do I=Isq,Ieq ; dowrite(I) = .false. ; enddo
do I=Isq,Ieq ; dowrite(I,j) = .false. ; enddo
if (CS%CFL_based_trunc) then
do I=Isq,Ieq ; vel_report(:) = 3.0e8 ; enddo ! Speed of light default.
do I=Isq,Ieq ; vel_report(i,j) = 3.0e8 ; enddo ! Speed of light default.
do k=1,nz ; do I=Isq,Ieq
if (u(I,j,k) < 0.0) then
CFL = (-u(I,j,k) * dt) * (G%dy_Cu(I,j) * G%IareaT(i+1,j))
Expand All @@ -1147,23 +1149,21 @@ subroutine vertvisc_limit_vel(u, v, h, ADp, CDp, fluxes, visc, dt, G, CS)
endif
if (CFL > CS%CFL_trunc) trunc_any = .true.
if (CFL > CS%CFL_report) then
dowrite(I) = .true.
vel_report(I) = MIN(vel_report(I), abs(u(I,j,k)))
dowrite(I,j) = .true.
vel_report(I,j) = MIN(vel_report(I,j), abs(u(I,j,k)))
endif
enddo ; enddo
else
do I=Isq,Ieq; vel_report(I,j) = maxvel; enddo
do k=1,nz ; do I=Isq,Ieq ; if (abs(u(I,j,k)) > maxvel) then
dowrite(I) = .true. ; trunc_any = .true.
dowrite(I,j) = .true. ; trunc_any = .true.
endif ;enddo ; enddo
endif

do I=Isq,Ieq ; if (dowrite(I)) then
! Here the diagnostic reporting subroutines are called if
! unphysically large values were found.
call write_u_accel(I, j, u, h, ADp, CDp, dt, G, CS%PointAccel_CSp, &
vel_report(I), -vel_report(I), (dt*fluxes%taux(I,j)/G%Rho0), &
a=CS%a_u(:,j,:), hv=CS%h_u(:,j,:))
endif ; enddo
do I=Isq,Ieq ; if (dowrite(I,j)) then
u_old(i,j,:) = u(i,j,:)
endif; enddo

if (trunc_any) then ; if (CS%CFL_based_trunc) then
do k=1,nz ; do I=Isq,Ieq
if ((u(I,j,k) * (dt * G%dy_Cu(I,j))) * G%IareaT(i+1,j) < -CS%CFL_trunc) then
Expand All @@ -1183,7 +1183,7 @@ subroutine vertvisc_limit_vel(u, v, h, ADp, CDp, fluxes, visc, dt, G, CS)
enddo ! j-loop
else
if (CS%CFL_based_trunc) then
!$OMP parallel do default(shared)
!$OMP parallel do default(shared)
do k=1,nz ; do j=js,je ; do I=Isq,Ieq
if ((u(I,j,k) * (dt * G%dy_Cu(I,j))) * G%IareaT(i+1,j) < -CS%CFL_trunc) then
u(I,j,k) = (-0.9*CS%CFL_trunc) * (G%areaT(i+1,j) / (dt * G%dy_Cu(I,j)))
Expand All @@ -1194,22 +1194,33 @@ subroutine vertvisc_limit_vel(u, v, h, ADp, CDp, fluxes, visc, dt, G, CS)
endif
enddo ; enddo ; enddo
else
!$OMP parallel do default(shared)
!$OMP parallel do default(shared)
do k=1,nz ; do j=js,je ; do I=Isq,Ieq ; if (abs(u(I,j,k)) > maxvel) then
u(I,j,k) = SIGN(truncvel,u(I,j,k))
if (h(i,j,k) + h(i+1,j,k) > 6.0*G%Angstrom) CS%ntrunc = CS%ntrunc + 1
endif ; enddo ; enddo ; enddo
endif
endif

if (len_trim(CS%u_trunc_file) > 0) then
do j=js,je; do I=Isq,Ieq ; if (dowrite(I,j)) then
! Here the diagnostic reporting subroutines are called if
! unphysically large values were found.
call write_u_accel(I, j, u_old, h, ADp, CDp, dt, G, CS%PointAccel_CSp, &
vel_report(I,j), -vel_report(I,j), (dt*fluxes%taux(I,j)/G%Rho0), &
a=CS%a_u(:,j,:), hv=CS%h_u(:,j,:))
endif ; enddo; enddo
endif


if (len_trim(CS%v_trunc_file) > 0) then
if (.not.CS%CFL_based_trunc) vel_report(:) = maxvel
!$OMP parallel do default(private) shared(Jsq,Jeq,is,ie,nz,CS,G,fluxes,v,h)
!$OMP parallel do default(private) shared(Jsq,Jeq,is,ie,nz,CS,G,fluxes,v,h,dt,maxvel,ADp,CDp,truncvel, &
!$OMP v_old,vel_report,dowrite)
do J=Jsq,Jeq
trunc_any = .false.
do i=is,ie ; dowrite(i) = .false. ; enddo
do i=is,ie ; dowrite(i,J) = .false. ; enddo
if (CS%CFL_based_trunc) then
do i=is,ie ; vel_report(:) = 3.0e8 ; enddo ! Speed of light default.
do i=is,ie ; vel_report(i,J) = 3.0e8 ; enddo ! Speed of light default.
do k=1,nz ; do i=is,ie
if (v(i,J,k) < 0.0) then
CFL = (-v(i,J,k) * dt) * (G%dx_Cv(i,J) * G%IareaT(i,j+1))
Expand All @@ -1218,23 +1229,21 @@ subroutine vertvisc_limit_vel(u, v, h, ADp, CDp, fluxes, visc, dt, G, CS)
endif
if (CFL > CS%CFL_trunc) trunc_any = .true.
if (CFL > CS%CFL_report) then
dowrite(i) = .true.
vel_report(i) = MIN(vel_report(i), abs(v(i,J,k)))
dowrite(i,J) = .true.
vel_report(i,J) = MIN(vel_report(i,J), abs(v(i,J,k)))
endif
enddo ; enddo
else
do i=is,ie ; vel_report(i,J) = maxvel ; enddo
do k=1,nz ; do i=is,ie ; if (abs(v(i,J,k)) > maxvel) then
dowrite(i) = .true. ; trunc_any = .true.
dowrite(i,J) = .true. ; trunc_any = .true.
endif ; enddo ; enddo
endif

do i=is,ie ; if (dowrite(i)) then
! Here the diagnostic reporting subroutines are called if
! unphysically large values were found.
call write_v_accel(i, J, v, h, ADp, CDp, dt, G, CS%PointAccel_CSp, &
vel_report(I), -vel_report(I), (dt*fluxes%tauy(i,J)/G%Rho0), &
a=CS%a_v(:,J,:),hv=CS%h_v(:,J,:))
do i=is,ie ; if (dowrite(i,J)) then
v_old(i,J,:) = v(i,J,:)
endif ; enddo

if (trunc_any) then ; if (CS%CFL_based_trunc) then
do k=1,nz; do i=is,ie
if ((v(i,J,k) * (dt * G%dx_Cv(i,J))) * G%IareaT(i,j+1) < -CS%CFL_trunc) then
Expand Down Expand Up @@ -1265,14 +1274,24 @@ subroutine vertvisc_limit_vel(u, v, h, ADp, CDp, fluxes, visc, dt, G, CS)
endif
enddo ; enddo ; enddo
else
!$OMP parallel do default(shared)
!$OMP parallel do default(shared)
do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie ; if (abs(v(i,J,k)) > maxvel) then
v(i,J,k) = SIGN(truncvel,v(i,J,k))
if (h(i,j,k) + h(i,j+1,k) > 6.0*G%Angstrom) CS%ntrunc = CS%ntrunc + 1
endif ; enddo ; enddo ; enddo
endif
endif

if (len_trim(CS%v_trunc_file) > 0) then
do J=Jsq,Jeq; do i=is,ie ; if (dowrite(i,J)) then
! Here the diagnostic reporting subroutines are called if
! unphysically large values were found.
call write_v_accel(i, J, v_old, h, ADp, CDp, dt, G, CS%PointAccel_CSp, &
vel_report(i,J), -vel_report(i,J), (dt*fluxes%tauy(i,J)/G%Rho0), &
a=CS%a_v(:,J,:),hv=CS%h_v(:,J,:))
endif ; enddo; enddo
endif

end subroutine vertvisc_limit_vel

subroutine vertvisc_init(MIS, Time, G, param_file, diag, ADp, dirs, ntrunc, CS)
Expand Down