Commit 35bf6e86 authored by bjb's avatar bjb
Browse files

OpenMP threading initial implementation

parent 1d234e88
......@@ -11,5 +11,7 @@ DEBUG_FLAGS = -g -C -check -fpe0 -traceback -mp1
#PROF_FLAGS = -qp -p -O0
PROF_FLAGS = -p -O2 -mp1
PROD_FLAGS = -O3 -i-static -mp1
OMP_FLAGS = -openmp -openmp-link static -openmp-threadprivate legacy
OMPMOD=true
DEFINES += -DREAL_4B=real\(4\)
endif
!$Id: bottom_friction.F90,v 1.9 2009-08-21 07:26:26 bjb Exp $
!$Id: bottom_friction.F90,v 1.10 2009-09-30 11:28:44 bjb Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -36,7 +36,8 @@
use parameters, only: kappa,avmmol
use domain, only: imin,imax,jmin,jmax,au,av,min_depth
use variables_2d
use getm_timers, only: tic, toc, TIM_BOTTFRICT
use getm_timers, only: tic, toc, TIM_BOTTFRICT
!$ use omp_lib
IMPLICIT NONE
!
! !INPUT PARAMETERS:
......@@ -61,7 +62,6 @@
Ncall = Ncall+1
write(debug,*) 'bottom_friction() # ',Ncall
#endif
CALL tic(TIM_BOTTFRICT)
#ifdef DEBUG
......@@ -106,6 +106,12 @@
end if
#endif
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j)
! The x-direction
!$OMP DO SCHEDULE(RUNTIME)
do j=jmin,jmax
do i=imin,imax
if (au(i,j) .gt. 0) then
......@@ -118,49 +124,61 @@
end if
end do
end do
!OMP END DO NOWAIT
! The x-direction
!$OMP DO SCHEDULE(RUNTIME)
do j=jmin-HALO,jmax+HALO
do i=imin-HALO,imax+HALO
if (au(i,j) .gt. 0) then
uloc(i,j) = U(i,j)/DU(i,j)
HH(i,j)=max(min_depth,DU(i,j))
#ifndef DEBUG
where (au .gt. 0)
uloc=U/DU
HH=max(min_depth,DU)
ruu=(kappa/log((zub+_HALF_*HH)/zub))**2
end where
ruu(i,j)=(kappa/log((zub(i,j)+_HALF_*HH(i,j))/zub(i,j)))**2
#else
uloc=U/DU
HH=max(min_depth,DU)
ruu=(zub+_HALF_*HH)/zub
do j=jmin,jmax
do i=imin,imax
if (ruu(i,j) .le. _ONE_) then
STDERR 'bottom_friction friction coefficient infinite.'
STDERR 'i = ',i,' j = ',j
STDERR 'min_depth = ',min_depth,' DU = ',DU(i,j)
stop
ruu(i,j)=(zub(i,j)+_HALF_*HH(i,j))/zub(i,j)
if (ruu(i,j) .le. _ONE_) then
!$OMP CRITICAL
STDERR 'bottom_friction friction coefficient infinite.'
STDERR 'i = ',i,' j = ',j
STDERR 'min_depth = ',min_depth,' DU = ',DU(i,j)
!$OMP END CRITICAL
stop
end if
ruu(i,j)=(kappa/log(ruu(i,j)))**2
#endif
end if
end do
end do
where (au .gt. 0)
ruu=(kappa/log(ruu))**2
end where
#endif
!$OMP END DO
if (runtype .eq. 1) then
where (au .gt. 0)
fricvel=sqrt(ruu*(uloc**2+vloc**2))
zub=min(HH,zub0+_TENTH_*avmmol/max(avmmol,fricvel))
ruu=(zub+_HALF_*HH)/zub
ruu=(kappa/log(ruu))**2
end where
!$OMP DO SCHEDULE(RUNTIME)
do j=jmin-HALO,jmax+HALO
do i=imin-HALO,imax+HALO
if (au(i,j) .gt. 0) then
fricvel(i,j)=sqrt(ruu(i,j)*(uloc(i,j)**2+vloc(i,j)**2))
zub(i,j)=min(HH(i,j),zub0(i,j)+_TENTH_*avmmol/max(avmmol,fricvel(i,j)))
ruu(i,j)=(zub(i,j)+_HALF_*HH(i,j))/zub(i,j)
ruu(i,j)=(kappa/log(ruu(i,j)))**2
ru(i,j)=ruu(i,j)*sqrt(uloc(i,j)**2+vloc(i,j)**2)
end if
end do
end do
!$OMP END DO
else
!$OMP DO SCHEDULE(RUNTIME)
do j=jmin-HALO,jmax+HALO
do i=imin-HALO,imax+HALO
if (au(i,j) .gt. 0) then
ru(i,j)=ruu(i,j)*sqrt(uloc(i,j)**2+vloc(i,j)**2)
end if
end do
end do
!$OMP END DO
end if
where (au .gt. 0)
ru=ruu*sqrt(uloc**2+vloc**2)
end where
! The y-direction
!$OMP DO SCHEDULE(RUNTIME)
do j=jmin,jmax
do i=imin,imax
if (av(i,j) .gt. 0) then
......@@ -173,49 +191,62 @@
end if
end do
end do
!$OMP END DO NOWAIT
!$OMP DO SCHEDULE(RUNTIME)
do j=jmin-HALO,jmax+HALO
do i=imin-HALO,imax+HALO
if (av(i,j) .gt. 0) then
vloc(i,j)=V(i,j)/DV(i,j)
HH(i,j)=max(min_depth,DV(i,j))
#ifndef DEBUG
where (av .gt. 0)
vloc=V/DV
HH=max(min_depth,DV)
rvv=(kappa/log((zvb+_HALF_*HH)/zvb))**2
end where
rvv(i,j)=(kappa/log((zvb(i,j)+_HALF_*HH(i,j))/zvb(i,j)))**2
#else
vloc=V/DV
HH=max(min_depth,DV)
rvv=(zvb+_HALF_*HH)/zvb
do j=jmin,jmax
do i=imin,imax
if (rvv(i,j) .le. _ONE_) then
STDERR 'bottom_friction friction coefficient infinite.'
STDERR 'i = ',i,' j = ',j
STDERR 'min_depth = ',min_depth,' DV = ',DU(i,j)
stop
rvv(i,j)=(zvb(i,j)+_HALF_*HH(i,j))/zvb(i,j)
if (rvv(i,j) .le. _ONE_) then
!$OMP CRITICAL
STDERR 'bottom_friction friction coefficient infinite.'
STDERR 'i = ',i,' j = ',j
STDERR 'min_depth = ',min_depth,' DV = ',DU(i,j)
!$OMP END CRITICAL
stop
end if
rvv(i,j)=(kappa/log(rvv(i,j)))**2
#endif
end if
end do
end do
where (av .gt. 0)
rvv=(kappa/log(rvv))**2
end where
#endif
!$OMP END DO
if (runtype .eq. 1) then
where (av .gt. 0)
fricvel=sqrt(rvv*(uloc**2+vloc**2))
zvb=min(HH,zvb0+_TENTH_*avmmol/max(avmmol,fricvel))
rvv=(zvb+_HALF_*HH)/zvb
rvv=(kappa/log(rvv))**2
end where
!$OMP DO SCHEDULE(RUNTIME)
do j=jmin-HALO,jmax+HALO
do i=imin-HALO,imax+HALO
if (av(i,j) .gt. 0) then
fricvel(i,j)=sqrt(rvv(i,j)*(uloc(i,j)**2+vloc(i,j)**2))
zvb(i,j)=min(HH(i,j),zvb0(i,j)+_TENTH_*avmmol/max(avmmol,fricvel(i,j)))
rvv(i,j)=(zvb(i,j)+_HALF_*HH(i,j))/zvb(i,j)
rvv(i,j)=(kappa/log(rvv(i,j)))**2
rv(i,j)=rvv(i,j)*sqrt(uloc(i,j)**2+vloc(i,j)**2)
end if
end do
end do
!$OMP END DO
else
!$OMP DO SCHEDULE(RUNTIME)
do j=jmin-HALO,jmax+HALO
do i=imin-HALO,imax+HALO
if (av(i,j) .gt. 0) then
rv(i,j)=rvv(i,j)*sqrt(uloc(i,j)**2+vloc(i,j)**2)
end if
end do
end do
!$OMP END DO
end if
where (av .gt. 0)
rv=rvv*sqrt(uloc**2+vloc**2)
end where
!$OMP END PARALLEL
CALL toc(TIM_BOTTFRICT)
#ifdef DEBUG
write(debug,*) 'Leaving bottom_friction()'
write(debug,*)
......
!$Id: depth_update.F90,v 1.11 2009-08-21 07:26:26 bjb Exp $
!$Id: depth_update.F90,v 1.12 2009-09-30 11:28:44 bjb Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -26,6 +26,7 @@
use domain, only: az,au,av,dry_z,dry_u,dry_v
use variables_2d, only: D,z,zo,DU,zu,DV,zv
use getm_timers, only: tic, toc, TIM_DPTHUPDATE
!$ use omp_lib
IMPLICIT NONE
!
! !INPUT PARAMETERS:
......@@ -39,7 +40,7 @@
!
! !LOCAL VARIABLES:
integer :: i,j
REALTYPE :: d1,x,d2i
REALTYPE :: d1,d2i,x
!EOP
!-----------------------------------------------------------------------
!BOC
......@@ -50,12 +51,24 @@
#endif
CALL tic(TIM_DPTHUPDATE)
! TODO/BJB: Why is this turned off?
#undef USE_MASK
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,d1,d2i,x)
! Depth in elevation points
D = z+H
!$OMP DO SCHEDULE(RUNTIME)
do j=jmin-HALO,jmax+HALO
do i=imin-HALO,imax+HALO
! TODO/BJB: Is it enough to do this on az?
D(i,j) = z(i,j)+H(i,j)
end do
end do
!$OMP END DO NOWAIT
! U-points
!$OMP DO SCHEDULE(RUNTIME)
do j=jmin-HALO,jmax+HALO
do i=imin-HALO,imax+HALO-1
#ifdef USE_MASK
......@@ -69,8 +82,10 @@
#endif
end do
end do
!$OMP END DO NOWAIT
! V-points
!$OMP DO SCHEDULE(RUNTIME)
do j=jmin-HALO,jmax+HALO-1
do i=imin-HALO,imax+HALO
#ifdef USE_MASK
......@@ -84,18 +99,28 @@
#endif
end do
end do
!$OMP END DO
d1 = 2*min_depth
d2i = _ONE_/(crit_depth-2*min_depth)
where (az .gt. 0)
dry_z = max(_ZERO_,min(_ONE_,(D-_HALF_*d1)*d2i))
end where
where (au .gt. 0)
dry_u = max(_ZERO_,min(_ONE_,(DU-d1)*d2i))
end where
where (av .gt. 0)
dry_v = max(_ZERO_,min(_ONE_,(DV-d1)*d2i))
end where
!$OMP DO SCHEDULE(RUNTIME)
do j=jmin-HALO,jmax+HALO
do i=imin-HALO,imax+HALO
if (az(i,j) .gt. 0) then
dry_z(i,j)=max(_ZERO_,min(_ONE_,(D(i,j)-_HALF_*d1)*d2i))
end if
if (au(i,j) .gt. 0) then
dry_u(i,j) = max(_ZERO_,min(_ONE_,(DU(i,j)-d1)*d2i))
end if
if (av(i,j) .gt. 0) then
dry_v(i,j) = max(_ZERO_,min(_ONE_,(DV(i,j)-d1)*d2i))
end if
end do
end do
!$OMP END DO
!$OMP END PARALLEL
#ifdef SLICE_MODEL
do i=imin,imax
......
!$Id: m2d.F90,v 1.33 2009-08-21 08:56:34 bjb Exp $
!$Id: m2d.F90,v 1.34 2009-09-30 11:28:44 bjb Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -354,6 +354,7 @@
end subroutine integrate_2d
!EOC
!-----------------------------------------------------------------------
!BOP
!
......
!$Id: momentum.F90,v 1.13 2009-08-18 10:24:43 bjb Exp $
!$Id: momentum.F90,v 1.14 2009-09-30 11:28:44 bjb Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -19,7 +19,6 @@
use domain, only: imin,imax,jmin,jmax
! For timer here: Only clock what is not taken at "next" level.
use getm_timers, only: tic, toc, TIM_MOMENTUM
IMPLICIT NONE
!
! !INPUT PARAMETERS:
......@@ -36,7 +35,7 @@
! Original author(s): Hans Burchard & Karsten Bolding
!
! !LOCAL VARIABLES:
logical :: ufirst=.false.
logical :: ufirst=.false.
!EOP
!-----------------------------------------------------------------------
!BOC
......@@ -134,6 +133,7 @@
use variables_2d, only: D,z,UEx,U,DU,fV,SlUx,SlRu,ru,fU,DV
use getm_timers, only: tic, toc, TIM_MOMENTUMH
use halo_zones, only : update_2d_halo,wait_halo,U_TAG
!$ use omp_lib
IMPLICIT NONE
!
! !INPUT PARAMETERS:
......@@ -147,9 +147,10 @@
integer :: i,j
REALTYPE :: zx(E2DFIELD)
REALTYPE :: Slr(E2DFIELD),tausu(E2DFIELD)
REALTYPE :: zp,zm,Uloc,Uold
REALTYPE :: zp,zm,Uloc
REALTYPE :: gamma=rho_0*g
REALTYPE :: cord_curv=_ZERO_
REALTYPE :: gammai
!EOP
!-----------------------------------------------------------------------
!BOC
......@@ -159,26 +160,50 @@
write(debug,*) 'umomentum() # ',Ncall
#endif
gammai = _ONE_/gamma
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,zp,zm,Uloc)
!$OMP DO SCHEDULE(RUNTIME)
do j=jmin,jmax
do i=imin,imax
if (au(i,j) .gt. 0) then
zp=max(z(i+1,j),-H(i ,j)+min(min_depth,D(i+1,j)))
zm=max(z(i ,j),-H(i+1,j)+min(min_depth,D(i ,j)))
zx(i,j)=(zp-zm+(airp(i+1,j)-airp(i,j))/gamma)/DXU
zx(i,j)=(zp-zm+(airp(i+1,j)-airp(i,j))*gammai)/DXU
! BJB-TODO: Change 0.5 -> _HALF_
tausu(i,j)=0.5*(tausx(i,j)+tausx(i+1,j))
end if
end do
end do
!$OMP END DO
where (U .gt. 0)
Slr=max(Slru, _ZERO_ )
else where
Slr=min(Slru, _ZERO_ )
end where
where ((au .eq. 1) .or. (au .eq. 2))
U=(U-dtm*(g*DU*zx+dry_u*(-tausu/rho_0-fV+UEx+SlUx+Slr)))/(1+dtm*ru/DU)
end where
! where (U .gt. 0)
! Slr=max(Slru, _ZERO_ )
! else where
! Slr=min(Slru, _ZERO_ )
! end where
! where ((au .eq. 1) .or. (au .eq. 2))
! U=(U-dtm*(g*DU*zx+dry_u*(-tausu/rho_0-fV+UEx+SlUx+Slr(i,j))))/(1+dtm*ru/DU)
! end where
!$OMP DO SCHEDULE(RUNTIME)
do j=jmin-HALO,jmax+HALO
do i=imin-HALO,imax+HALO
if (U(i,j) .gt. 0) then
Slr(i,j)=max(Slru(i,j), _ZERO_ )
else
Slr(i,j)=min(Slru(i,j), _ZERO_ )
end if
if ((au(i,j) .eq. 1) .or. (au(i,j) .eq. 2)) then
U(i,j)=(U(i,j)-dtm*(g*DU(i,j)*zx(i,j)+dry_u(i,j)*&
(-tausu(i,j)/rho_0-fV(i,j)+UEx(i,j)+SlUx(i,j)+Slr(i,j))))/&
(_ONE_+dtm*ru(i,j)/DU(i,j))
end if
end do
end do
!$OMP END DO
!$OMP END PARALLEL
! The rest of this sub is not easy to thread.
#ifdef SLICE_MODEL
do i=imin,imax
......@@ -199,12 +224,13 @@
if(av(i,j) .ge. 1) then
! Espelid et al. [2000], IJNME 49, 1521-1545
#ifdef NEW_CORI
! BJB-TODO: Change 0.25 -> _QUART_
Uloc= &
( U(i,j )/sqrt(DU(i,j ))+ U(i-1,j )/sqrt(DU(i-1,j )) &
+ U(i,j+1)/sqrt(DU(i,j+1))+ U(i-1,j+1)/sqrt(DU(i-1,j+1))) &
*0.25*sqrt(DV(i,j))
#else
Uloc=0.25*( U(i,j)+ U(i-1,j)+ U(i,j+1)+ U(i-1,j+1))
Uloc=0.25*( U(i-1,j)+U(i,j)+U(i-1,j+1)+U(i,j+1))
#endif
#if defined(SPHERICAL) || defined(CURVILINEAR)
cord_curv=(V(i,j)*(DYX-DYXIM1)-Uloc*(DXCJP1-DXC)) &
......@@ -311,6 +337,7 @@
REALTYPE :: zp,zm,Vloc
REALTYPE :: gamma=rho_0*g
REALTYPE :: cord_curv=_ZERO_
REALTYPE :: gammai
!EOP
!-----------------------------------------------------------------------
!BOC
......@@ -320,26 +347,62 @@
write(debug,*) 'vmomentum() # ',Ncall
#endif
gammai = _ONE_/gamma
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,zp,zm,Vloc)
!$OMP DO SCHEDULE(RUNTIME)
do j=jmin,jmax
do i=imin,imax
if (av(i,j) .gt. 0) then
zp=max(z(i,j+1),-H(i,j )+min(min_depth,D(i,j+1)))
zm=max(z(i,j ),-H(i,j+1)+min(min_depth,D(i,j )))
zy(i,j)=(zp-zm+(airp(i,j+1)-airp(i,j))/gamma)/DYV
zy(i,j)=(zp-zm+(airp(i,j+1)-airp(i,j))*gammai)/DYV
! BJB-TODO: Change 0.5 -> _HALF_
tausv(i,j)=0.5*(tausy(i,j)+tausy(i,j+1))
end if
end do
end do
!$OMP END DO
where (V .gt. 0)
Slr=max(Slrv, _ZERO_ )
else where
Slr=min(Slrv, _ZERO_ )
end where
! where (V .gt. 0)
! Slr=max(Slrv, _ZERO_ )
! else where
! Slr=min(Slrv, _ZERO_ )
! end where
! $OMP DO SCHEDULE(RUNTIME)
! do j=jmin-HALO,jmax+HALO
! do i=imin-HALO,imax+HALO
! if (V(i,j).gt.0) then
! Slr(i,j)=max(Slrv(i,j), _ZERO_ )
! else
! Slr(i,j)=min(Slrv(i,j), _ZERO_ )
! end if
! end do
! end do
! $OMP END DO
where ((av .eq. 1) .or. (av .eq. 2))
V=(V-dtm*(g*DV*zy+dry_v*(-tausv/rho_0+fU+VEx+SlVx+Slr)))/(1+dtm*rv/DV)
end where
! where ((av .eq. 1) .or. (av .eq. 2))
! V=(V-dtm*(g*DV*zy+dry_v*(-tausv/rho_0+fU+VEx+SlVx+Slr)))/(1+dtm*rv/DV)
! end where
!$OMP DO SCHEDULE(RUNTIME)
do j=jmin-HALO,jmax+HALO
do i=imin-HALO,imax+HALO
if (V(i,j).gt.0) then
Slr(i,j)=max(Slrv(i,j), _ZERO_ )
else
Slr(i,j)=min(Slrv(i,j), _ZERO_ )
end if
if ((av(i,j) .eq. 1) .or. (av(i,j) .eq. 2)) then
V(i,j)=(V(i,j)-dtm*(g*DV(i,j)*zy(i,j)+dry_v(i,j)*&
(-tausv(i,j)/rho_0+fU(i,j)+VEx(i,j)+SlVx(i,j)+Slr(i,j))))/&
(_ONE_+dtm*rv(i,j)/DV(i,j))
end if
end do
end do
!$OMP END DO
!$OMP END PARALLEL
! The rest of this sub is not easy to thread.
#ifdef SLICE_MODEL
do i=imin,imax
......@@ -361,12 +424,13 @@
if(au(i,j) .ge. 1) then
! Espelid et al. [2000], IJNME 49, 1521-1545
#ifdef NEW_CORI
! BJB-TODO: Change 0.25 -> _QUART_
Vloc = &
( V(i,j )/sqrt(DV(i,j ))+ V(i+1,j )/sqrt(DV(i+1,j )) + &
V(i,j-1)/sqrt(DV(i,j-1))+ V(i+1,j-1)/sqrt(DV(i+1,j-1))) &
*0.25*sqrt(DU(i,j))
#else
Vloc = 0.25*( V(i,j)+ V(i+1,j)+ V(i,j-1)+ V(i+1,j-1))
Vloc = 0.25*( V(i,j-1)+ V(i+1,j-1)+V(i,j)+V(i+1,j))
#endif
#if defined(SPHERICAL) || defined(CURVILINEAR)
cord_curv=(Vloc*(DYCIP1-DYC)-U(i,j)*(DXX-DXXJM1)) &
......
!$Id: sealevel.F90,v 1.18 2009-08-18 10:24:43 bjb Exp $
!$Id: sealevel.F90,v 1.19 2009-09-30 11:28:44 bjb Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -35,6 +35,7 @@
use halo_zones, only : nprocs,set_flag,u_TAG,v_TAG
use variables_2d, only: break_mask,break_stat
use domain, only : min_depth,au,av
!$ use omp_lib
#endif
IMPLICIT NONE
!
......@@ -65,7 +66,14 @@
#endif
call tic(TIM_SEALEVEL)
zo = z
! OMP-NOTE: This loop does not really improve from threading.
! It is bound by memory access, and everything is nicely
! lined up anyway, so we keep it in serial.
do j=jmin-HALO,jmax+HALO
do i=imin-HALO,imax+HALO
zo(i,j) = z(i,j)
end do
end do
#ifdef USE_BREAKS