Skip to content

GitLab

G
getm
Commit e22480a8 authored by Knut's avatar Knut
Browse files
Options

phydis_2d + a lot of renamings

parent 52263789
Hide whitespace changes
Inline Side-by-side
Showing with 659 additions and 286 deletions
compilers/vs2008/getm.vfproj
View file @ e22480a8
......@@ -50,6 +50,7 @@
<File RelativePath="..\..\src\2d\depth_update.F90"/>
<File RelativePath="..\..\src\2d\m2d.F90"/>
<File RelativePath="..\..\src\2d\momentum.F90"/>
<File RelativePath="..\..\src\2d\physical_dissipation.F90"/>
<File RelativePath="..\..\src\2d\residual.F90"/>
<File RelativePath="..\..\src\2d\sealevel.F90"/>
<File RelativePath="..\..\src\2d\update_2d_bdy.F90"/>
......
compilers/vs2010/getm.vfproj
View file @ e22480a8
......@@ -41,6 +41,7 @@
<File RelativePath="..\..\src\2d\depth_update.F90"/>
<File RelativePath="..\..\src\2d\m2d.F90"/>
<File RelativePath="..\..\src\2d\momentum.F90"/>
<File RelativePath="..\..\src\2d\physical_dissipation.F90"/>
<File RelativePath="..\..\src\2d\residual.F90"/>
<File RelativePath="..\..\src\2d\sealevel.F90"/>
<File RelativePath="..\..\src\2d\update_2d_bdy.F90"/>
......
src/2d/Makefile
View file @ e22480a8
......@@ -10,7 +10,7 @@ LIB = $(LIBDIR)/lib2d${buildtype}.a
TEXSRC = m2d.F90 variables_2d.F90 \
advection.F90 adv_split_u.F90 adv_split_v.F90 adv_arakawa_j7_2dh.F90 adv_upstream_2dh.F90 adv_fct_2dh.F90 adv_tvd_limiter.F90 \
bottom_friction.F90 uv_advect.F90 uv_diffusion.F90 uv_diff_2dh.F90 momentum.F90 sealevel.F90 depth_update.F90 \
update_2d_bdy.F90 residual.F90 cfl_check.F90
update_2d_bdy.F90 physical_dissipation.F90 residual.F90 cfl_check.F90
MOD = \
${LIB}(variables_2d.o) \
......@@ -25,6 +25,8 @@ ${LIB}(momentum.o) \
${LIB}(sealevel.o) \
${LIB}(uv_advect.o) \
${LIB}(uv_diffusion.o) \
${LIB}(uv_diff_2dh.o) \
${LIB}(physical_dissipation.o) \
${LIB}(residual.o) \
${LIB}(update_2d_bdy.o) \
${LIB}(adv_split_u.o) \
......@@ -32,8 +34,7 @@ ${LIB}(adv_split_v.o) \
${LIB}(adv_upstream_2dh.o) \
${LIB}(adv_fct_2dh.o) \
${LIB}(adv_arakawa_j7_2dh.o) \
${LIB}(adv_tvd_limiter.o) \
${LIB}(uv_diff_2dh.o)
${LIB}(adv_tvd_limiter.o)
all: modules objects
......
src/2d/m2d.F90
View file @ e22480a8
......@@ -351,6 +351,15 @@
LEVEL1 'postinit_2d'
if (do_numerical_analyses_2d) then
allocate(phydis_2d(E2DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_2d: Error allocating memory (phydis_2d)'
phydis_2d = _ZERO_
allocate(numdis_2d(E2DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_2d: Error allocating memory (numdis_2d)'
numdis_2d = _ZERO_
end if
!
! It is possible that a user changes the land mask and reads an "old" hotstart file.
! In this case the "old" velocities will need to be zeroed out.
......@@ -461,6 +470,7 @@
call tic(TIM_INTEGR2D)
call uv_advect(U,V,DU,DV)
call uv_diffusion(An_method,U,V,D,DU,DV) ! Has to be called after uv_advect.
if (do_numerical_analyses_2d) call physical_dissipation(U,V,DU,DV,Am,phydis_2d)
call toc(TIM_INTEGR2D)
call momentum(loop,tausx,tausy,airp)
......
src/2d/physical_dissipation.F90 0 → 100644
View file @ e22480a8
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: physical_dissipation()
!
! !INTERFACE:
subroutine physical_dissipation(U,V,DU,DV,Am,phydiss)
!
! !DESCRIPTION:
!
! Here, the physical dissipation of mean kinetic energy is calculated as
! \begin{equation}
! D^{phy} = D^{phy}_h + D^{phy}_v = 2 \alpha A^M_h S_h^2 +(\nu_t+\nu) S_v^2
! \end{equation}
! with the horizontal eddy viscosity, $A^M_h$, the vertical eddy viscosity,
! $\nu_t$, the molecular viscosity, $\nu$, the drying parameter, $\alpha$
! (see equations (\ref{uEq}) and (\ref{vEq})), the horizontal shear tensor,
! \begin{equation}
! S_{ij} = \frac12 \left(\partial_{x_i} u_j + \partial_{x_j}u_i \right),
! \quad i,j=1,2
! \end{equation}
! (with $x_1=x$, $x_2=y$, $u_1=u$ and $u_2=u$) and the scalar vertical shear
! \begin{equation}
! S_v = \left[\left(\partial_zu\right)^2
! +\left(\partial_zv\right)^2\right]^{1/2}.
! \end{equation}
! The horizontal shear square is then denoted as
! \begin{equation}
! S_h^2 = \sum_{i,j=1}^2
! \left[\frac12 \left(\partial_{x_i} u_j + \partial_{x_j}u_i \right)\right]^2.
! \end{equation}
! Note that the dissipation term is obtained from the multiplication of the
! $u$-equation (\ref{uEq}) by $u$, the multiplicationof the $v$-equation
! (\ref{vEq}) by $v$ and subsequent addition of the two resulting equations.
!
! This routine is called only for {\tt save\_numerical\_analyses = .true.}
!
!
! !USES:
use domain, only: imin,imax,jmin,jmax,az,au,av,ax
#if defined(SPHERICAL) || defined(CURVILINEAR)
use domain, only: dxc,dyc,dxx,dyx
#else
use domain, only: dx,dy
#endif
use domain, only: dry_z
IMPLICIT NONE
! !INPUT PARAMETERS
REALTYPE,dimension(E2DFIELD),intent(in) :: U,V,DU,DV
REALTYPE,intent(in) :: Am
! !OUTPUT PARAMETERS
REALTYPE,dimension(E2DFIELD),intent(out) :: phydiss
!
! !REVISION HISTORY:
! Original author(s): Hans Burchard
!
! !LOCAL VARIABLES:
integer :: i,j,rc
REALTYPE,dimension(E2DFIELD) :: shear
REALTYPE,dimension(:,:),allocatable :: u_vel,v_vel
logical,save :: first=.true.
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
integer, save :: Ncall = 0
Ncall = Ncall+1
write(debug,*) 'physical_dissipation() # ',Ncall
#endif
#ifdef SLICE_MODEL
j = jmax/2 ! this MUST NOT be changed!!!
#endif
if (first) then
allocate(u_vel(E2DFIELD),stat=rc)
if (rc /= 0) stop 'physical_dissipation: Error allocating memory (u_vel)'
u_vel = _ZERO_
allocate(v_vel(E2DFIELD),stat=rc)
if (rc /= 0) stop 'physical_dissipation: Error allocating memory (v_vel)'
v_vel = _ZERO_
first = .false.
end if
if (Am .gt. _ZERO_) then
#ifndef SLICE_MODEL
do j=jmin-HALO,jmax+HALO
#endif
do i=imin-HALO,imax+HALO-1
if (au(i,j) .ge. 1) then
u_vel(i,j) = U(i,j)/DU(i,j)
end if
end do
#ifndef SLICE_MODEL
end do
#endif
#ifndef SLICE_MODEL
do j=jmin-HALO,jmax+HALO-1
#endif
do i=imin-HALO,imax+HALO
if (av(i,j) .ge. 1) then
v_vel(i,j) = V(i,j)/DV(i,j)
end if
end do
#ifndef SLICE_MODEL
end do
#endif
#ifndef SLICE_MODEL
do j=jmin-HALO,jmax+HALO-1
#endif
do i=imin-HALO,imax+HALO-1
shear(i,j) = _ZERO_
if (ax(i,j) .eq. 1) then
! calculate shearX
if (au(i,j).eq.1 .or. au(i,j+1).eq.1) then
! Note (KK): excludes concave and W/E open boundaries (dvdxX=0)
! includes convex open boundaries (no mirroring)
shear(i,j) = (v_vel(i+1,j) - v_vel(i,j)) / DXX
end if
#ifndef SLICE_MODEL
if (av(i,j).eq.1 .or. av(i+1,j).eq.1) then
! Note (KK): excludes concave and N/S open boundaries (dudyX=0)
! includes convex open boundaries (no mirroring)
shear(i,j) = shear(i,j) + (u_vel(i,j+1) - u_vel(i,j)) / DYX
end if
#endif
end if
end do
#ifndef SLICE_MODEL
end do
do j=jmin-HALO,jmax+HALO-1
#endif
do i=imax+HALO-1,imin-HALO+1 ! loop order MUST NOT be changed!!!
! Note (KK): slip condition dudyV(av=0)=0
! prolonged outflow condition dvdxV(av=3)=0
! shearV(av=3) would require shearX outside open boundary
! (however shearV(av=3) not needed, therefore not calculated)
if (az(i,j).eq.1 .or. az(i,j+1).eq.1) then
! calculate shearV
shear(i,j) = _HALF_ * ( shear(i-1,j) + shear(i,j) )
else
shear(i,j) = _ZERO_
end if
end do
#ifndef SLICE_MODEL
end do
do j=jmax+HALO-1,jmin-HALO+1 ! loop order MUST NOT be changed!!!
do i=imin-HALO+1,imax+HALO-1
if (az(i,j) .eq. 1) then
! calculate shearC
shear(i,j) = _HALF_ * ( shear(i,j-1) + shear(i,j) )
end if
end do
end do
#endif
! 2*Am*((du/dx)**2+0.5*(dv/dx + du/dy)**2+(dv/dy)**2) on T-points
#ifndef SLICE_MODEL
do j=jmin,jmax
#endif
do i=imin,imax
if (az(i,j) .eq. 1) then
phydiss(i,j) = _TWO_*Am* &
( &
((u_vel(i,j)-u_vel(i-1,j))/DXC)**2 &
#ifndef SLICE_MODEL
+((v_vel(i,j)-v_vel(i,j-1))/DYC)**2 &
#endif
+_HALF_*shear(i,j)**2 &
)
phydiss(i,j) = dry_z(i,j) * phydiss(i,j)
end if
end do
#ifndef SLICE_MODEL
end do
#endif
end if
#ifdef DEBUG
write(debug,*) 'Leaving physical_dissipation()'
write(debug,*)
#endif
return
end subroutine physical_dissipation
!EOC
!-----------------------------------------------------------------------
! Copyright (C) 2012 - Hans Burchard and Karsten Bolding !
!-----------------------------------------------------------------------
src/2d/uv_advect.F90
View file @ e22480a8
......@@ -24,6 +24,8 @@
#endif
use m2d, only: dtm,vel2d_adv_split,vel2d_adv_hor
use variables_2d, only: UEx,VEx
use variables_2d, only: do_numerical_analyses_2d
use variables_2d, only: numdis_2d
use advection, only: NOADV,UPSTREAM,J7,do_advection
use halo_zones, only: update_2d_halo,wait_halo,U_TAG,V_TAG
use getm_timers, only: tic,toc,TIM_UVADV,TIM_UVADVH
......@@ -38,10 +40,11 @@
! Original author(s): Hans Burchard & Karsten Bolding
!
! !LOCAL VARIABLES:
integer :: i,j
integer :: i,j
REALTYPE,dimension(E2DFIELD) :: fadv,Uadv,Vadv,DUadv,DVadv
REALTYPE,dimension(E2DFIELD),target :: Dadv
REALTYPE,dimension(:,:),pointer :: pDadv
REALTYPE,dimension(E2DFIELD) :: numdiss,vel2
!EOP
!-----------------------------------------------------------------------
!BOC
......@@ -179,12 +182,65 @@
call wait_halo(U_TAG)
call toc(TIM_UVADVH)
end if
!$OMP END SINGLE
if (do_numerical_analyses_2d) then
!$OMP DO SCHEDULE(RUNTIME)
#ifndef SLICE_MODEL
do j=jmin-HALO,jmax+HALO ! calculate square of u-velocity before advection step
#endif
do i=imin-HALO,imax+HALO
vel2(i,j) = fadv(i,j)**2
end do
#ifndef SLICE_MODEL
end do
#endif
!$OMP END DO
end if
!$OMP SINGLE
call do_advection(dtm,fadv,Uadv,Vadv,DUadv,DVadv,pDadv,pDadv, &
vel2d_adv_split,vel2d_adv_hor,_ZERO_,U_TAG, &
advres=UEx)
!$OMP END SINGLE
if (do_numerical_analyses_2d) then
!$OMP SINGLE
call do_advection(dtm,vel2,Uadv,Vadv,DUadv,DVadv,pDadv,pDadv, &
vel2d_adv_split,vel2d_adv_hor,_ZERO_,U_TAG)
!$OMP END SINGLE
!$OMP DO SCHEDULE(RUNTIME)
#ifndef SLICE_MODEL
do j=jmin,jmax ! calculate kinetic energy dissipaion rate for u-velocity
#endif
do i=imin,imax
numdiss(i,j) = ( vel2(i,j) - fadv(i,j)**2 ) / dtm
end do
#ifndef SLICE_MODEL
end do
#endif
!$OMP END DO
!$OMP SINGLE
call update_2d_halo(numdiss,numdiss,au,imin,jmin,imax,jmax,U_TAG)
call wait_halo(U_TAG)
!$OMP END SINGLE
!$OMP DO SCHEDULE(RUNTIME)
#ifndef SLICE_MODEL
do j=jmin,jmax
#endif
do i=imin,imax
numdis_2d(i,j) = _HALF_*( numdiss(i-1,j) + numdiss(i,j) )
end do
#ifndef SLICE_MODEL
end do
#endif
!$OMP END DO
end if
! Here begins dimensional split advection for v-velocity
......@@ -282,8 +338,7 @@
end if
!$OMP END PARALLEL
!$OMP SINGLE
if (vel2d_adv_hor.ne.UPSTREAM .and. vel2d_adv_hor.ne.J7) then
! we need to update fadv(imin-HALO:imax+HALO,jmax+HALO)
call tic(TIM_UVADVH)
......@@ -291,10 +346,72 @@
call wait_halo(V_TAG)
call toc(TIM_UVADVH)
end if
!$OMP END SINGLE
if (do_numerical_analyses_2d) then
!$OMP DO SCHEDULE(RUNTIME)
#ifndef SLICE_MODEL
do j=jmin-HALO,jmax+HALO ! calculate square of v-velocity before advection step
#endif
do i=imin-HALO,imax+HALO
vel2(i,j) = fadv(i,j)**2
end do
#ifndef SLICE_MODEL
end do
#endif
!$OMP END DO
end if
!$OMP SINGLE
call do_advection(dtm,fadv,Uadv,Vadv,DUadv,DVadv,pDadv,pDadv, &
vel2d_adv_split,vel2d_adv_hor,_ZERO_,V_TAG, &
advres=VEx)
!$OMP END SINGLE
if (do_numerical_analyses_2d) then
!$OMP SINGLE
call do_advection(dtm,vel2,Uadv,Vadv,DUadv,DVadv,pDadv,pDadv, &
vel2d_adv_split,vel2d_adv_hor,_ZERO_,V_TAG)
!$OMP END SINGLE
!$OMP DO SCHEDULE(RUNTIME)
#ifndef SLICE_MODEL
do j=jmin,jmax ! calculate kinetic energy dissipaion rate for v-velocity
#endif
do i=imin,imax
numdiss(i,j) = ( vel2(i,j) - fadv(i,j)**2 ) / dtm
end do
#ifndef SLICE_MODEL
end do
#endif
!$OMP END DO
!$OMP SINGLE
#ifdef SLICE_MODEL
numdiss(imin:imax,j-1) = numdiss(imin:imax,j)
#else
call update_2d_halo(numdiss,numdiss,av,imin,jmin,imax,jmax,V_TAG)
call wait_halo(V_TAG)
#endif
!$OMP END SINGLE
!$OMP DO SCHEDULE(RUNTIME)
#ifndef SLICE_MODEL
do j=jmin,jmax
#endif
do i=imin,imax
numdis_2d(i,j) = numdis_2d(i,j) &
+_HALF_*( numdiss(i,j-1) + numdiss(i,j) )
end do
#ifndef SLICE_MODEL
end do
#endif
!$OMP END DO
end if
!$OMP END PARALLEL
else ! if (vel2d_adv_hor .eq. NOADV)
......
src/2d/variables_2d.F90
View file @ e22480a8
......@@ -21,11 +21,16 @@
IMPLICIT NONE
!
! !PUBLIC DATA MEMBERS:
logical :: do_numerical_analyses_2d=.false.
#ifdef STATIC
#include "static_2d.h"
#else
#include "dynamic_declarations_2d.h"
#endif
REALTYPE,dimension(:,:),allocatable :: numdis_2d,phydis_2d
integer :: size2d_field
integer :: mem2d
!
......
src/3d/m3d.F90
View file @ e22480a8
......@@ -306,49 +306,49 @@
LEVEL1 'postinit_3d'
if (do_numerical_analyses) then
allocate(phydis3d(I3DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (phydis3d)'
phydis3d = _ZERO_
allocate(phydis2d(I2DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (phydis2d)'
phydis2d = _ZERO_
allocate(numdis3d(I3DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (numdis3d)'
numdis3d = _ZERO_
allocate(numdis2d(I2DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (numdis2d)'
numdis2d = _ZERO_
if (do_numerical_analyses_3d) then
allocate(phydis_3d(I3DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (phydis_3d)'
phydis_3d = _ZERO_
allocate(phydis_int(I2DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (phydis_int)'
phydis_int = _ZERO_
allocate(numdis_3d(I3DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (numdis_3d)'
numdis_3d = _ZERO_
allocate(numdis_int(I2DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (numdis_int)'
numdis_int = _ZERO_
if (calc_temp) then
allocate(phymix3d_T(I3DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (phymix3d_T)'
phymix3d_T = _ZERO_
allocate(phymix2d_T(I2DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (phymix2d_T)'
phymix2d_T = _ZERO_
allocate(nummix3d_T(I3DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (nummix3d_T)'
nummix3d_T = _ZERO_
allocate(nummix2d_T(I2DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (nummix2d_T)'
nummix2d_T = _ZERO_
allocate(phymix_T(I3DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (phymix_T)'
phymix_T = _ZERO_
allocate(phymix_T_int(I2DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (phymix_T_int)'
phymix_T_int = _ZERO_
allocate(nummix_T(I3DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (nummix_T)'
nummix_T = _ZERO_
allocate(nummix_T_int(I2DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (nummix_T_int)'
nummix_T_int = _ZERO_
end if
if (calc_salt) then
allocate(phymix3d_S(I3DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (phymix3d_S)'
phymix3d_S = _ZERO_
allocate(phymix2d_S(I2DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (phymix2d_S)'
phymix2d_S = _ZERO_
allocate(nummix3d_S(I3DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (nummix3d_S)'
nummix3d_S = _ZERO_
allocate(nummix2d_S(I2DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (nummix2d_S)'
nummix2d_S = _ZERO_
allocate(phymix_S(I3DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (phymix_S)'
phymix_S = _ZERO_
allocate(phymix_S_int(I2DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (phymix_S_int)'
phymix_S_int = _ZERO_
allocate(nummix_S(I3DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (nummix_S)'
nummix_S = _ZERO_
allocate(nummix_S_int(I2DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (nummix_S_int)'
nummix_S_int = _ZERO_
end if
end if
......@@ -535,7 +535,7 @@
! is present, where it can be skipped.
! We need SS: 1) #if (!defined(CONSTANT_VISCOSITY) && !defined(PARABOLIC_VISCOSITY))
! 2) adpative coordinates
! 3) if(do_numerical_analyses)
! 3) if(do_numerical_analyses_3d)
! We need NN (runtype .ge. 3):
! 1) #if (!defined(CONSTANT_VISCOSITY) && !defined(PARABOLIC_VISCOSITY))
! 2) adaptive coordinates
......@@ -562,7 +562,7 @@
end if
if (do_numerical_analyses) call physical_dissipation_3d(Am,phydis3d,phydis2d)
if (do_numerical_analyses_3d) call physical_dissipation_3d(Am,phydis_3d,phydis_int)
#ifndef NO_BAROCLINIC
if(runtype .eq. 4) then ! prognostic T and S
......
src/3d/physical_dissipation_3d.F90
View file @ e22480a8
......@@ -5,7 +5,7 @@
! !ROUTINE: physical_dissipation_3d()
!
! !INTERFACE:
subroutine physical_dissipation_3d(Am,pd3d,pd2d)
subroutine physical_dissipation_3d()
!
! !DESCRIPTION:
!
......@@ -40,25 +40,19 @@
! !USES:
use domain, only: imin,imax,jmin,jmax,kmax,H,au,av,ax,az
#if defined(SPHERICAL) || defined(CURVILINEAR)
use domain, only: dxx,dyx,dxc,dyc
use domain, only: dxc,dyc,dxx,dyx
#else
use domain, only: dx,dy,dry_z
use domain, only: dx,dy
#endif
use domain, only: dry_z
use variables_3d, only: num,uu,vv,hn,hun,hvn,SS
use m2d, only: Am
use variables_3d, only: phydis_3d,phydis_int,num,uu,vv,hn,hun,hvn,SS
use parameters, only: avmmol
IMPLICIT NONE
! !INPUT PARAMETERS
REALTYPE, intent(in) :: Am
! !INPUT/OUTPUT PARAMETERS
REALTYPE, intent(out) :: pd3d(I3DFIELD)
REALTYPE, intent(out) :: pd2d(I2DFIELD)
!
! !REVISION HISTORY:
! Original author(s): Hannes Rennau
! Original author(s): Hans Burchard
!
! !LOCAL VARIABLES:
REALTYPE :: dupper,dlower
......@@ -67,68 +61,63 @@
!EOP
!-----------------------------------------------------------------------
!BOC
pd3d=_ZERO_
pd2d=_ZERO_
#ifdef DEBUG
integer, save :: Ncall = 0
Ncall = Ncall+1
write(debug,*) 'physical_dissipation() # ',Ncall
#endif
#ifdef SLICE_MODEL
j = jmax/2 ! this MUST NOT be changed!!!
#endif
phydis_int = _ZERO_
if (Am .gt. _ZERO_) then