Commit bf2d76ee authored by kbk's avatar kbk
Browse files

added numerical diffusion - An - to namelist - passed to uv_diffusion()

parent 6daf4ec8
!$Id: m2d.F90,v 1.7 2003-08-28 10:28:40 kbk Exp $
!$Id: m2d.F90,v 1.8 2003-09-13 10:00:51 kbk Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -26,7 +26,7 @@
!
! !PUBLIC DATA MEMBERS:
logical :: have_boundaries
REALTYPE :: dtm, z0_const=0.010, Am=-_ONE_
REALTYPE :: dtm, z0_const=0.010,Am=-_ONE_,An=-_ONE_
integer :: MM=1,residual=-1
logical :: bdy2d=.false.
integer :: bdyfmt_2d,bdytype,bdyramp_2d=-1
......@@ -41,7 +41,10 @@
! Original author(s): Karsten Bolding & Hans Burchard
!
! $Log: m2d.F90,v $
! Revision 1.7 2003-08-28 10:28:40 kbk
! Revision 1.8 2003-09-13 10:00:51 kbk
! added numerical diffusion - An - to namelist - passed to uv_diffusion()
!
! Revision 1.7 2003/08/28 10:28:40 kbk
! explict setting UEx and VEx to 0 every micro time step
!
! Revision 1.6 2003/08/15 12:47:40 kbk
......@@ -134,7 +137,7 @@
! !LOCAL VARIABLES:
integer :: rc
namelist /m2d/ &
MM,z0_const,Am,residual,bdy2d,bdyfmt_2d,bdyramp_2d,bdyfile_2d
MM,z0_const,Am,An,residual,bdy2d,bdyfmt_2d,bdyramp_2d,bdyfile_2d
!EOP
!-------------------------------------------------------------------------
!BOC
......@@ -162,7 +165,10 @@
#endif
if (Am .lt. _ZERO_) then
LEVEL2 'Am is less than zero ---> horizontal diffusion not included'
LEVEL2 'Am < 0 --> horizontal momentum diffusion not included'
end if
if (An .lt. _ZERO_) then
LEVEL2 'An < 0 --> numerical momentum diffusion not included'
end if
LEVEL2 'Open boundary=',bdy2d
......@@ -254,8 +260,8 @@
#else
#ifndef UV_ADV_DIRECT
call uv_advect()
if (Am .gt. _ZERO_) then
call uv_diffusion(Am) ! Has to be called after uv_advect.
if (Am .gt. _ZERO_ .or. An .gt. _ZERO_) then
call uv_diffusion(Am,An) ! Has to be called after uv_advect.
end if
#endif
#endif
......
!$Id: uv_diffusion.F90,v 1.3 2003-08-28 10:40:42 kbk Exp $
!$Id: uv_diffusion.F90,v 1.4 2003-09-13 10:00:51 kbk Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -6,7 +6,7 @@
! !ROUTINE: uv_diffusion() - diffusion of momentum.
!
! !INTERFACE:
subroutine uv_diffusion(Am)
subroutine uv_diffusion(Am,An)
!
! !DESCRIPTION:
!
......@@ -21,7 +21,7 @@
IMPLICIT NONE
!
! !INPUT PARAMETERS:
REALTYPE, intent(in) :: Am
REALTYPE, intent(in) :: Am,An
!
! !INPUT/OUTPUT PARAMETERS:
!
......@@ -32,11 +32,6 @@
!
! !LOCAL VARIABLES:
integer :: i,j
! REALTYPE :: AN=100.
REALTYPE :: AN=1000.
! REALTYPE :: AN=_ZERO_
! REALTYPE :: AN=10000.
! REALTYPE :: AN=100000.
!EOP
!-----------------------------------------------------------------------
!BOC
......@@ -46,15 +41,18 @@
write(debug,*) 'uv_diffusion() # ',Ncall
#endif
! Central for dx(2*Am*dx(uu^2/hun))
! Central for dx(2*Am*dx(U^2/HU))
do j=jmin,jmax
do i=imin,imax+1 ! PP defined on T-points
PP(i,j)=_ZERO_
if (az(i,j) .ge. 1) then
PP(i,j)=2.*Am*DYC*D(i,j) &
*(U(i,j)/DU(i,j)-U(i-1,j)/DU(i-1,j))/DXC
PP(i,j)=PP(i,j)+AN*DYC*(U(i,j)-U(i-1,j))/DXC
else
PP(i,j)=_ZERO_
if(Am .gt. _ZERO_) then
PP(i,j)=2.*Am*DYC*D(i,j) &
*(U(i,j)/DU(i,j)-U(i-1,j)/DU(i-1,j))/DXC
end if
if(An .gt. _ZERO_) then
PP(i,j)=PP(i,j)+An*DYC*(U(i,j)-U(i-1,j))/DXC
end if
end if
end do
end do
......@@ -69,13 +67,16 @@
! Central for dy(Am*(dy(U^2/DU)+dx(V^2/DV)))
do j=jmin-1,jmax ! PP defined on X-points
do i=imin,imax
PP(i,j)=_ZERO_
if (ax(i,j) .ge. 1) then
PP(i,j)=Am*0.5*(DU(i,j)+DU(i,j+1))*DXX &
*((U(i,j+1)/DU(i,j+1)-U(i,j)/DU(i,j))/DYX &
+(V(i+1,j)/DV(i+1,j)-V(i,j)/DV(i,j))/DXX )
PP(i,j)=PP(i,j)+AN*(V(i+1,j)-V(i,j))
else
PP(i,j)=_ZERO_
if(Am .gt. _ZERO_) then
PP(i,j)=Am*0.5*(DU(i,j)+DU(i,j+1))*DXX &
*((U(i,j+1)/DU(i,j+1)-U(i,j)/DU(i,j))/DYX &
+(V(i+1,j)/DV(i+1,j)-V(i,j)/DV(i,j))/DXX )
end if
if(An .gt. _ZERO_) then
PP(i,j)=PP(i,j)+An*(V(i+1,j)-V(i,j))
end if
end if
end do
end do
......@@ -90,13 +91,16 @@
! Central for dx(Am*(dy(U^2/DU)+dx(V^2/DV)))
do j=jmin,jmax ! PP defined on X-points
do i=imin-1,imax
PP(i,j)=_ZERO_
if (ax(i,j) .ge. 1) then
PP(i,j)=Am*0.5*(DV(i,j)+DV(i+1,j))*DYX &
*((U(i,j+1)/DU(i,j+1)-U(i,j)/DU(i,j))/DYX &
+(V(i+1,j)/DV(i+1,j)-V(i,j)/DV(i,j))/DXX )
PP(i,j)=PP(i,j)+AN*(U(i,j+1)-U(i,j))
else
PP(i,j)=_ZERO_
if(Am .gt. _ZERO_) then
PP(i,j)=Am*0.5*(DV(i,j)+DV(i+1,j))*DYX &
*((U(i,j+1)/DU(i,j+1)-U(i,j)/DU(i,j))/DYX &
+(V(i+1,j)/DV(i+1,j)-V(i,j)/DV(i,j))/DXX )
end if
if(An .gt. _ZERO_) then
PP(i,j)=PP(i,j)+An*(U(i,j+1)-U(i,j))
end if
end if
end do
end do
......@@ -111,12 +115,15 @@
! Central for dy(2*Am*dy(V^2/DV))
do j=jmin,jmax+1 ! PP defined on T-points
do i=imin,imax
PP(i,j)=_ZERO_
if (az(i,j) .ge. 1) then
PP(i,j)=2.*Am*DXC*D(i,j) &
*(V(i,j)/DV(i,j)-V(i,j-1)/DV(i,j-1))/DYC
PP(i,j)=PP(i,j)+AN*DXC*(V(i,j)-V(i,j-1))/DYC
else
PP(i,j)=_ZERO_
if(Am .gt. _ZERO_) then
PP(i,j)=2.*Am*DXC*D(i,j) &
*(V(i,j)/DV(i,j)-V(i,j-1)/DV(i,j-1))/DYC
end if
if(An .gt. _ZERO_) then
PP(i,j)=PP(i,j)+An*DXC*(V(i,j)-V(i,j-1))/DYC
end if
end if
end do
end do
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment