!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
module topostress 3,13
!BOP
! !MODULE: topostress
! !DESCRIPTION:
! This module contains routines necessary for computing stress
! due to bottom topography.
!
! !REVISION HISTORY:
! SVN:$Id: topostress.F90 12674 2008-10-31 22:21:32Z njn01 $
! !USES:
use POP_KindsMod
use POP_ErrorMod
use POP_FieldMod
use POP_GridHorzMod
use POP_HaloMod
use domain
use blocks
use distribution
use constants
use io
use grid
use broadcast
use exit_mod
implicit none
private
save
! !PUBLIC MEMBER FUNCTIONS:
public :: init_topostress
! !PUBLIC DATA MEMBERS:
logical (POP_logical), public :: &
ltopostress ! true if topographic stress desired
real (POP_r8), dimension(:,:,:), allocatable, public :: &
TSU, TSV ! topographic stress velocities
!EOP
!BOC
!EOC
!***********************************************************************
contains
!***********************************************************************
!BOP
! !IROUTINE: init_topostress
! !INTERFACE:
subroutine init_topostress(errorCode) 1,6
! !DESCRIPTION:
! This routine allocates stress arrays if topographic stress is
! chosen and initializes topo stress parameters.
!
! !REVISION HISTORY:
! same as module
!
! !OUTPUT PARAMETERS:
integer (POP_i4), intent(out) :: &
errorCode
!EOP
!BOC
!-----------------------------------------------------------------------
!
! input namelist to choose topostress
!
!-----------------------------------------------------------------------
integer (POP_i4) :: &
nsmooth_topo ! number of passes of topography smoother
integer (POP_i4) :: nml_error ! namelist i/o error flag
namelist /topostress_nml/ltopostress, nsmooth_topo
!-----------------------------------------------------------------------
!
! read namelist to see if topostress desired
!
!-----------------------------------------------------------------------
errorCode = POP_Success
ltopostress = .false.
nsmooth_topo= 0
if (my_task == master_task) then
open (nml_in, file=nml_filename, status='old',iostat=nml_error)
if (nml_error /= 0) then
nml_error = -1
else
nml_error = 1
endif
do while (nml_error > 0)
read(nml_in, nml=topostress_nml,iostat=nml_error)
end do
if (nml_error == 0) close(nml_in)
endif
call broadcast_scalar(nml_error, master_task)
if (nml_error /= 0) then
call exit_POP(sigAbort,'ERROR reading topostress_nml')
endif
if (my_task == master_task) then
write(stdout,blank_fmt)
write(stdout,ndelim_fmt)
write(stdout,blank_fmt)
write(stdout,'(a27)') ' Topographic stress options'
write(stdout,blank_fmt)
write(stdout,delim_fmt)
if (ltopostress) then
write(stdout,'(a27)') ' Topographic stress enabled'
if (nsmooth_topo > 0) then
write(stdout,'(a26,i2,a8)') ' Topography smoothed with ', &
nsmooth_topo,' passes.'
else
write(stdout,'(a25)') ' Topography not smoothed.'
endif
else
write(stdout,'(a28)') ' Topographic stress disabled'
endif
endif
call broadcast_scalar(ltopostress, master_task)
call broadcast_scalar(nsmooth_topo, master_task)
!-----------------------------------------------------------------------
!
! allocate the topographic stress velocity arrays if required
!
!-----------------------------------------------------------------------
if (ltopostress) then
allocate (TSU(nx_block,ny_block,nblocks_clinic), &
TSV(nx_block,ny_block,nblocks_clinic))
call topo_stress(nsmooth_topo, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'init_topostress: error in topo_stress')
return
endif
endif
!-----------------------------------------------------------------------
!EOC
end subroutine init_topostress
!***********************************************************************
!BOP
! !IROUTINE: topo_stress
! !INTERFACE:
subroutine topo_stress(nsmooth_topo, errorCode) 1,9
! !DESCRIPTION:
! Calculate topographic stress (maximum entropy) velocities. These
! are time-independent 2d fields given by:
! \begin{eqnarray}
! \Psi^* &=& -f L^2 H \\
! H U^* &=& -\nabla_y(\Psi^*) \\
! H V^* &=& +\nabla_x(\Psi^*)
! \end{eqnarray}
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
integer (POP_i4), intent(in) :: &
nsmooth_topo ! number of passes of topography smoother
! !OUTPUT PARAMETERS:
integer (POP_i4), intent(out) :: &
errorCode ! returned error code
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (POP_i4) :: &
i,j,iter,iblock ! local iteration counters
real (POP_r8), dimension(nx_block,ny_block) :: &
TSP, &! topo stress streamfunction
SCALE, &! scale length
HTOLD ! old topography
real (POP_r8), dimension(nx_block,ny_block,nblocks_clinic) :: &
HTNEW ! smoothed topography
real (POP_r8), parameter :: &
tslse = 12.0e5_POP_r8, &!
tslsp = 3.0e5_POP_r8 !
type (block) :: &
this_block ! block information for current block
!-----------------------------------------------------------------------
!
! smooth topography if requested
!
!-----------------------------------------------------------------------
errorCode = POP_Success
do iblock = 1,nblocks_clinic
HTNEW(:,:,iblock) = HT(:,:,iblock) ! initialize
end do
do iter = 1, nsmooth_topo
do iblock = 1,nblocks_clinic
this_block = get_block(blocks_clinic(iblock),iblock)
HTOLD = HTNEW(:,:,iblock)
call smooth_topo2(HTOLD,HTNEW(:,:,iblock),this_block)
end do
call POP_HaloUpdate(HTNEW, POP_haloClinic, &
POP_gridHorzLocCenter, POP_fieldKindScalar, &
errorCode, fillValue = 0.0_POP_r8)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'topo_stress: error updating halo for htnew')
return
endif
enddo
!-----------------------------------------------------------------------
!
! calculate the topographic stress equilibrium stream function.
!
!-----------------------------------------------------------------------
do iblock = 1,nblocks_clinic
this_block = get_block(blocks_clinic(iblock),iblock)
SCALE = tslsp + (tslse - tslsp)* &
(p5 + p5*cos(c2*TLAT(:,:,iblock)))
where (KMT(:,:,iblock) /= 0)
TSP = -FCORT(:,:,iblock)*SCALE*SCALE*HTNEW(:,:,iblock)
elsewhere
TSP = c0
end where
!-----------------------------------------------------------------------
!
! calculate the topographic stress velocities from stream function.
! compute gradient with 4 point stencil
!
!-----------------------------------------------------------------------
TSU(:,:,iblock) = c0
TSV(:,:,iblock) = c0
do j=this_block%jb,this_block%je
do i=this_block%ib,this_block%ie
TSV(i,j,iblock) = DXUR(i,j,iblock)*p5*HUR(i,j,iblock)* &
(TSP(i+1,j+1) - TSP(i ,j) - &
TSP(i ,j+1) + TSP(i+1,j))
TSU(i,j,iblock) = -DYUR(i,j,iblock)*p5*HUR(i,j,iblock)* &
(TSP(i+1,j+1) - TSP(i ,j) + &
TSP(i ,j+1) - TSP(i+1,j))
end do
end do
where (KMU(:,:,iblock) == 0)
TSV(:,:,iblock) = c0 ! zero at land points
TSU(:,:,iblock) = c0
endwhere
! apply only in 'deep' water
! where (KMU(:,:,iblock) <= 3)
! TSU(:,:,iblock) = c0
! TSV(:,:,iblock) = c0
! endwhere
end do ! block loop
call POP_HaloUpdate(TSU, POP_haloClinic, &
POP_gridHorzLocNECorner, POP_fieldKindVector, &
errorCode, fillValue = 0.0_POP_r8)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'topo_stress: error updating halo for tsu')
return
endif
call POP_HaloUpdate(TSV, POP_haloClinic, &
POP_gridHorzLocNECorner, POP_fieldKindVector, &
errorCode, fillValue = 0.0_POP_r8)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'topo_stress: error updating halo for tsv')
return
endif
!-----------------------------------------------------------------------
!EOC
end subroutine topo_stress
!***********************************************************************
!BOP
! !IROUTINE: smooth_topo2
! !INTERFACE:
subroutine smooth_topo2(HTOLD,HTNEW,this_block) 1
! !DESCRIPTION:
! This routine smooths topography using a 9-point averaging stencil
! given by
! \begin{equation}
! \begin{array}{ccccc}
! 1 & -- & 2 & -- & 1 \\
! | & & | & & | \\
! 2 & -- & 4 & -- & 2 \\
! | & & | & & | \\
! 1 & -- & 2 & -- & 1
! \end{array} \nonumber
! \end{equation}
! Land points are not included in the smoothing, and the
! stencil is modified to include only ocean points in the
! averaging. This routine is nearly identical to the smooth
! topography routine in the grid module.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (POP_r8), dimension(nx_block,ny_block), intent(in) :: &
HTOLD ! old HT field to be smoothed
type (block), intent(in) :: &
this_block ! block info for this sub block
! !OUTPUT PARAMETERS:
real (POP_r8), dimension(nx_block,ny_block), intent(out) :: &
HTNEW ! smoothed HT field
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (POP_i4) :: &
i,j, &! loop counters
bid ! local block index
integer (POP_i4), dimension(nx_block,ny_block) :: &
NB, &! array to compute number of ocean neighbors
IWORK ! local work space
real (POP_r8), dimension(nx_block,ny_block) :: &
WORK ! local work space
!-----------------------------------------------------------------------
!
! smooth topography
!
!-----------------------------------------------------------------------
bid = this_block%local_id
where (KMT(:,:,bid) > 0)
IWORK = 1
HTNEW = HTOLD
elsewhere
IWORK = 0
HTNEW = c0
endwhere
do j=this_block%jb,this_block%je
do i=this_block%ib,this_block%ie
WORK(i,j) = c4*HTNEW(i,j) + &
c2*HTNEW(i+1,j ) + c2*HTNEW(i-1,j ) + &
c2*HTNEW(i ,j+1) + c2*HTNEW(i ,j-1) + &
HTNEW(i+1,j+1) + HTNEW(i+1,j-1) + &
HTNEW(i-1,j+1) + HTNEW(i-1,j-1)
NB(i,j) = c4*IWORK(i,j) + &
c2*IWORK(i+1,j ) + c2*IWORK(i-1,j ) + &
c2*IWORK(i ,j+1) + c2*IWORK(i ,j-1) + &
IWORK(i+1,j+1) + IWORK(i+1,j-1) + &
IWORK(i-1,j+1) + IWORK(i-1,j-1)
end do
end do
!-----------------------------------------------------------------------
!
! new depth field
!
!-----------------------------------------------------------------------
where ((KMT(:,:,bid) /= 0) .and. (NB /= 0))
HTNEW = WORK/real(NB)
elsewhere
HTNEW = c0
endwhere
!-----------------------------------------------------------------------
!EOC
end subroutine smooth_topo2
!***********************************************************************
end module topostress
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||