Skip to content

User/z1l/openmp #32

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 4 commits into from
Jul 23, 2014
Merged

User/z1l/openmp #32

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
6732a07
Improve single thread performance
Zhi-Liang Jul 15, 2014
8c502e3
fix to reproduce old results
Zhi-Liang Jul 21, 2014
58230a3
fix for OBC
Zhi-Liang Jul 22, 2014
c20118f
Merge remote-tracking branch 'origin' into user/z1l/openmp
Zhi-Liang Jul 23, 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
61 changes: 40 additions & 21 deletions src/parameterizations/vertical/MOM_KPP.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -340,6 +340,8 @@ subroutine KPP_calculate(CS, G, h, Temp, Salt, u, v, EOS, uStar, buoyFlux, Kt, K
real, dimension( G%ke+1, 2) :: Kdiffusivity ! Vertical diffusivity at interfaces (m2/s)
real, dimension( G%ke+1 ) :: Kviscosity ! Vertical viscosity at interfaces (m2/s)
real, dimension( G%ke+1, 2) :: nonLocalTrans ! Non-local transport for heat/salt at interfaces (non-dimensional)
real, dimension( G%ke ) :: surfBuoyFlux2
real, dimension( 3*G%ke ) :: rho_1D, pres_1D, Temp_1D, Salt_1D
real :: kOBL, OBLdepth_0d, surfFricVel, surfBuoyFlux, Coriolis
real :: GoRho, pRef, rho1, rhoK, rhoKm1, Uk, Vk, const1, Cv, sigma
real :: zBottomMinusOffset ! Height of bottom plus a little bit (m)
Expand All @@ -358,7 +360,7 @@ subroutine KPP_calculate(CS, G, h, Temp, Salt, u, v, EOS, uStar, buoyFlux, Kt, K
real :: hTotV, surfHv, surfV ! Thickness, integral and average of V over the surface layer
real :: hTotUm1, surfHum1, surfUm1 ! Same fo i-1,j or U
real :: hTotVm1, surfHvm1, surfVm1 ! Same fo i,j-1 of V

integer :: kk
#ifdef __DO_SAFETY_CHECKS__
if (CS%debug) then
call hchksum(h, "KPP in: h",G,haloshift=0)
Expand Down Expand Up @@ -392,7 +394,8 @@ subroutine KPP_calculate(CS, G, h, Temp, Salt, u, v, EOS, uStar, buoyFlux, Kt, K
!$OMP rho1,rhoK,rhoKm1,deltaRho,N2_1d,N_1d,delH, &
!$OMP surfBuoyFlux,Ws_1d,Cv,Vt2_1d,BulkRi_1d, &
!$OMP OBLdepth_0d,zBottomMinusOffset,Kdiffusivity, &
!$OMP Kviscosity,sigma,kOBL)
!$OMP Kviscosity,sigma,kOBL,kk,pres_1D,Temp_1D, &
!$OMP Salt_1D,rho_1D,surfBuoyFlux2)
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 All @@ -417,7 +420,16 @@ subroutine KPP_calculate(CS, G, h, Temp, Salt, u, v, EOS, uStar, buoyFlux, Kt, K
pRef = 0.
do k = 1, G%ke
km1 = max(1, k-1)

kk = 3*(k-1)
pres_1D(kk+1) = pRef
pres_1D(kk+2) = pRef
pres_1D(kk+3) = pRef
Temp_1D(kk+1) = surfTemp
Temp_1D(kk+2) = Temp(i,j,k)
Temp_1D(kk+3) = Temp(i,j,km1)
Salt_1D(kk+1) = surfsalt
Salt_1D(kk+2) = Salt(i,j,k)
Salt_1D(kk+3) = Salt(i,j,km1)
! Compute heights, referenced to the surface (z=0)
cellHeight(k) = iFaceHeight(k) - 0.5 * h(i,j,k) * G%H_to_m ! cell center in metres
iFaceHeight(k+1) = iFaceHeight(k) - h(i,j,k) * G%H_to_m ! cell bottom in metres
Expand All @@ -431,16 +443,7 @@ subroutine KPP_calculate(CS, G, h, Temp, Salt, u, v, EOS, uStar, buoyFlux, Kt, K
Uk = 0.5 * ( abs( u(i,j,k) - surfU ) + abs( u(i-1,j,k) - surfUm1 ) ) ! delta_k U w/ C-grid average
Vk = 0.5 * ( abs( v(i,j,k) - surfV ) + abs( v(i,j-1,k) - surfVm1 ) ) ! delta_k V w/ C-grid average
deltaU2(k) = Uk**2 + Vk**2
call calculate_density(surfTemp, surfSalt, pRef, rho1, EOS)
call calculate_density(Temp(i,j,k), Salt(i,j,k), pRef, rhoK, EOS)
call calculate_density(Temp(i,j,km1), Salt(i,j,km1), pRef, rhoKm1, EOS)
deltaRho(k) = rhoK - rho1
! N2 is on interfaces.
N2_1d(k) = ( GoRho * (rhoK - rhoKm1) ) / (0.5*(h(i,j,km1) + h(i,j,k))+G%H_subroundoff) ! Can be negative
! N = sqrt(N^2) but because N^2 can be negative, we clip N^2 before taking the square root ??????
N_1d(k) = sqrt( max( N2_1d(k), 0.) )
! N_1d(k) = sqrt( abs( N2_1d(k)) ) ! From MOM4p1 ?????
! N_1d(k) = sign( sqrt( abs( N2_1d(k) ) ), N2_1d(k) ) ! Suggestion ????

! Pressure at bottom of level k will become pressure at top of level on next iteration
pRef = pRef + G%g_Earth * G%Rho0 * h(i,j,k) * G%H_to_m ! Boussinesq approximation!!!! ?????

Expand Down Expand Up @@ -472,23 +475,39 @@ subroutine KPP_calculate(CS, G, h, Temp, Salt, u, v, EOS, uStar, buoyFlux, Kt, K
surfHvm1 = surfHvm1 + v(i,j-1,k) * delH ; surfVm1 = surfHvm1 / hTotVm1
endif
enddo ! k
call calculate_density(Temp_1D, Salt_1D, pres_1D, rho_1D, 1, 3*G%ke, EOS)
do k = 1, G%ke
km1 = max(1, k-1)
kk = 3*(k-1)
deltaRho(k) = rho_1D(kk+2) - rho_1D(kk+1)
! N2 is on interfaces.
N2_1d(k) = ( GoRho * (rho_1D(kk+2) - rho_1D(kk+3)) ) / (0.5*(h(i,j,km1) + h(i,j,k))+G%H_subroundoff) ! Can be negative
! N = sqrt(N^2) but because N^2 can be negative, we clip N^2 before taking the square root ??????
N_1d(k) = sqrt( max( N2_1d(k), 0.) )
! N_1d(k) = sqrt( abs( N2_1d(k)) ) ! From MOM4p1 ?????
! N_1d(k) = sign( sqrt( abs( N2_1d(k) ) ), N2_1d(k) ) ! Suggestion ????
enddo
N2_1d( G%ke+1 ) = 0.
N_1d( G%ke+1 ) = 0.

do k = 1, G%ke
surfBuoyFlux2(k) = buoyFlux(i,j,1) - buoyFlux(i,j,k+1) ! This difference accounts for penetrating of SW
enddo

call CVmix_kpp_compute_turbulent_scales( &
eps, & ! (in) Normalized boundary layer depth; sigma = eps
-cellHeight, & ! (in) Guess that OBL depth (m) = -cellHeight(k)
surfBuoyFlux2, & ! (in) Buoyancy flux at surface (m2/s3)
surfFricVel, & ! (in) Turbulent friction velocity at surface (m/s)
w_s=Ws_1d, & ! (out) Turbulent velocity scale profile (m/s)
CVmix_kpp_params_user=CS%KPP_params ) ! KPP parameters

! Estimate Ws in order to estimate Vt^2 (eq. 23 in LMD94)
do k = 1, G%ke

! Note that if sigma > eps, then CVmix_kpp_compute_turbulent_scales
! computes w_s and w_m velocity scale at sigma=eps. So we only pass
! sigma=eps for this calculation.
surfBuoyFlux = buoyFlux(i,j,1) - buoyFlux(i,j,k+1) ! This difference accounts for penetrating of SW
call CVmix_kpp_compute_turbulent_scales( &
eps, & ! (in) Normalized boundary layer depth; sigma = eps
-cellHeight(k), & ! (in) Guess that OBL depth (m) = -cellHeight(k)
surfBuoyFlux, & ! (in) Buoyancy flux at surface (m2/s3)
surfFricVel, & ! (in) Turbulent friction velocity at surface (m/s)
w_s=Ws_1d(k), & ! (out) Turbulent velocity scale profile (m/s)
CVmix_kpp_params_user=CS%KPP_params ) ! KPP parameters

! Unresolved squared velocity, Vt^2, eq 23 from LMD94
Cv = max( 1.7, 2.1 - 200. * N_1d(k) ) ! Cv from eq A3 of Danabasoglu et al. 2006
Expand Down