Commit 3fe49cfb authored by kbk's avatar kbk
Browse files

removed some temporary devel files

parent e0a37edf
!$Id: coordinates.F90.devel,v 1.1 2003-04-07 17:07:06 kbk Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: coordinates() - defines the vertical coordinate system.
!
! !INTERFACE:
subroutine coordinates(cord_type,cord_relax)
!
! !DESCRIPTION:
!
! !USES:
use domain, only: iimin,iimax,jjmin,jjmax,kmax
use domain, only: H,HU,HV,az,au,av,min_depth
use domain, only: ga,ddu,ddl,d_gamma
use variables_3d, only: dt,kmin,kumin,kvmin,ho,hn,huo,hun,hvo,hvn
use variables_3d, only: sseo,ssen,ssuo,ssun,ssvo,ssvn
use halo_zones, only: update_3d_halo,wait_halo,H_TAG,HU_TAG,HV_TAG
IMPLICIT NONE
!
! !INPUT PARAMETERS:
integer, intent(in) :: cord_type
REALTYPE, intent(in) :: cord_relax
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !REVISION HISTORY:
! Original author(s): Hans Burchard & Karsten Bolding
!
! $Log: coordinates.F90.devel,v $
! Revision 1.1 2003-04-07 17:07:06 kbk
! should be fixed and used later
!
!
! !LOCAL VARIABLES:
integer :: i,j,k,rc,kk
REALTYPE :: alpha,maxdepth
LOGICAL :: gamma_surf=.true.
REALTYPE :: HH,zz,r
logical, save :: first=.true.,equiv_sigma=.false.
REALTYPE, save :: kmaxm1
REALTYPE, save, dimension(:), allocatable :: dga,be,sig
REALTYPE, save, dimension(:,:,:), allocatable :: gga
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
integer, save :: Ncall = 0
Ncall = Ncall+1
write(debug,*) 'coordinates() # ',Ncall
#endif
if (first) then
first = .false.
allocate(ga(0:kmax),stat=rc)
if (rc /= 0) stop 'coordinates: Error allocating (ga)'
allocate(dga(0:kmax),stat=rc)
if (rc /= 0) STOP 'coordinates: Error allocating (dga)'
ga = _ZERO_; dga = _ZERO_
select case (cord_type)
case (1) ! sigma coordinates
if (ddu .le. _ZERO_ .and. ddl .le. _ZERO_) then
equiv_sigma=.true.
kmaxm1= _ONE_/float(kmax)
ga(0) = -_ONE_
do k=1,kmax
ga(k) = ga(k-1) + _ONE_/kmax
dga(k)=ga(k)-ga(k-1)
end do
ga(kmax) = _ZERO_
else
! Non-equidistant sigma coordinates
! This zooming routine is from Antoine Garapon, ICCH, DK
ga(0)= -_ONE_
dga(0)= _ZERO_
do k=1,kmax
ga(k)=tanh((ddl+ddu)*k/float(kmax)-ddl)+tanh(ddl)
ga(k)=ga(k)/(tanh(ddl)+tanh(ddu)) - _ONE_
dga(k)=ga(k)-ga(k-1)
end do
end if
kmin=1
kumin=1
kvmin=1
do j=jjmin-1,jjmax+1
do i=iimin-1,iimax+1
hn(i,j,1:kmax)=max(sseo(i,j)+H(i,j),min_depth)*dga(1:kmax)
hun(i,j,1:kmax)=max(ssuo(i,j)+HU(i,j),min_depth)*dga(1:kmax)
hvn(i,j,1:kmax)=max(ssvo(i,j)+HV(i,j),min_depth)*dga(1:kmax)
end do
end do
case (2) ! z-level
STDERR 'z-level not coded yet'
stop 'coordinates'
case (3) ! general vertical coordinates
do k=0,kmax
ga(k) = k
end do
allocate(sig(0:kmax),stat=rc) ! dimensionless sigma-coordinate
if (rc /= 0) STOP 'coordinates: Error allocating (sig)'
allocate(be(0:kmax),stat=rc) ! dimensionless beta-coordinate
if (rc /= 0) STOP 'coordinates: Error allocating (be)'
allocate(gga(I3DFIELD),stat=rc) ! dimensionless gamma-coordinate
if (rc /= 0) stop 'coordinates: Error allocating memory (gga)'
sig = _ZERO_ ; be = _ZERO_ ; gga = _ZERO_
do k=1,kmax
be(k)=tanh((ddl+ddu)*k/float(kmax)-ddl)+tanh(ddl)
be(k)=be(k)/(tanh(ddl)+tanh(ddu))-_ONE_
sig(k)=k/float(kmax)-_ONE_
end do
if (gamma_surf) then
kk=kmax
else
kk=1
end if
do j=jjmin-1,jjmax+1
do i=iimin-1,iimax+1
HH=max(H(i,j),min_depth)
alpha=min( &
((be(kk)-be(kk-1))-D_gamma/HH &
*(sig(kk)-sig(kk-1))) &
/((be(kk)-be(kk-1))-(sig(kk)-sig(kk-1))),_ONE_)
gga(i,j,0)=-_ONE_
do k=1,kmax
gga(i,j,k)=alpha*sig(k)+(1.-alpha)*be(k)
if (gga(i,j,k) .lt. gga(i,j,k-1)) then
STDERR kk,(be(kk)-be(kk-1)),(sig(kk)-sig(kk-1))
STDERR D_gamma,HH
STDERR alpha
STDERR k-1,gga(i,j,k-1),be(k-1),sig(k-1)
STDERR k,gga(i,j,k),be(k),sig(k)
stop 'coordinates'
end if
end do
end do
end do
! call update_3d_halo(gga,gga,az,iimin,jjmin,iimax,jjmax,kmax,H_TAG)
kmin=1
kumin=1
kvmin=1
! call wait_halo(H_TAG)
! Initial layer distribution - uses new as temporary var
do k=1,kmax
do j=jjmin,jjmax
do i=iimin-1,iimax
HH=max(ssuo(i,j)+HU(i,j),min_depth)
hun(i,j,k)=HH*0.5* &
(gga(i,j,k)-gga(i,j,k-1)+gga(i+1,j,k)-gga(i+1,j,k-1))
end do
end do
end do
call update_3d_halo(hun,hun,au,iimin,jjmin,iimax,jjmax,kmax,HU_TAG)
do k=1,kmax
do j=jjmin-1,jjmax
do i=iimin,iimax
HH=max(ssvo(i,j)+HV(i,j),min_depth)
hvn(i,j,k)=HH*0.5* &
(gga(i,j,k)-gga(i,j,k-1)+gga(i,j+1,k)-gga(i,j+1,k-1))
end do
end do
end do
call wait_halo(HU_TAG)
call update_3d_halo(hvn,hvn,av,iimin,jjmin,iimax,jjmax,kmax,HV_TAG)
do k=1,kmax
do j=jjmin-1,jjmax+1
do i=iimin-1,iimax+1
HH=max(sseo(i,j)+H(i,j),min_depth)
hn(i,j,k)=HH*(gga(i,j,k)-gga(i,j,k-1))
end do
end do
end do
call wait_halo(HV_TAG)
case default
end select
end if ! first
ho = hn
huo = hun
hvo = hvn
#if 0
STDERR H(iimin,jjmin),H(iimax,jjmin),H(iimax,jjmax),H(iimin,jjmax)
STDERR H(iimin-1,jjmin-1),H(iimax+1,jjmin-1),H(iimax+1,jjmax+1),H(iimin-1,jjmax+1)
STDERR hn(iimin-1,jjmin-1,1:kmax:2)
STDERR hn(iimax+1,jjmin-1,1:kmax:2)
STDERR hn(iimax+1,jjmax+1,1:kmax:2)
STDERR hn(iimin-1,jjmax+1,1:kmax:2)
#endif
select case (cord_type)
case (1) ! sigma coordinates
if (equiv_sigma) then
do j=jjmin-1,jjmax+1
do i=iimin-1,iimax+1
hn(i,j,:)=max(ssen(i,j)+H(i,j),min_depth)*kmaxm1
end do
end do
do j=jjmin-1,jjmax+1
do i=iimin-1,iimax+1
hun(i,j,:)=max(ssun(i,j)+HU(i,j),min_depth)*kmaxm1
end do
end do
do j=jjmin-1,jjmax+1
do i=iimin-1,iimax+1
hvn(i,j,:)=max(ssvn(i,j)+HV(i,j),min_depth)*kmaxm1
end do
end do
else ! non-equivdistant
do j=jjmin-1,jjmax+1
do i=iimin-1,iimax+1
hn(i,j,1:kmax)=max(ssen(i,j)+H(i,j),min_depth)*dga(1:kmax)
end do
end do
do j=jjmin-1,jjmax+1
do i=iimin-1,iimax+1
hun(i,j,1:kmax)=max(ssun(i,j)+HU(i,j),min_depth)*dga(1:kmax)
end do
end do
do j=jjmin-1,jjmax+1
do i=iimin-1,iimax+1
hvn(i,j,1:kmax)=max(ssvn(i,j)+HV(i,j),min_depth)*dga(1:kmax)
end do
end do
end if
case (2) ! z-level
case (3) ! general vertical coordinates
! The general vertical coordinates can be relaxed towards the new layer
! thicknesses by the following relaxation time scale r. This should
! later be generalised also for sigma coordinates.
! hun and hvn needs to communicate with neighbors therefore they are
! done befoe hn
maxdepth=600. ! needs to be calculated later ... HB
do j=jjmin,jjmax
do i=iimin-1,iimax
if (au(i,j).gt.0) then
r = cord_relax/dt*HU(i,j)/maxdepth
zz = -HU(i,j)
do k=1,kmax-1
HH=max(ssun(i,j)+HU(i,j),min_depth)
hun(i,j,k)=(huo(i,j,k)*r+HH*0.5*(gga(i,j,k)-gga(i,j,k-1) &
+gga(i+1,j,k)-gga(i+1,j,k-1)))/(r+1.)
zz = zz + hun(i,j,k)
end do
hun(i,j,kmax)=ssun(i,j)-zz
end if
end do
end do
call update_3d_halo(hun,hun,au,iimin,jjmin,iimax,jjmax,kmax,HU_TAG)
do j=jjmin-1,jjmax
do i=iimin,iimax
if (av(i,j).gt.0) then
r = cord_relax/dt*HV(i,j)/maxdepth
zz = -HV(i,j)
do k=1,kmax-1
HH=max(ssvn(i,j)+HV(i,j),min_depth)
hvn(i,j,k)=(hvo(i,j,k)*r+HH*0.5*(gga(i,j,k)-gga(i,j,k-1) &
+gga(i,j+1,k)-gga(i,j+1,k-1)))/(r+1.)
zz=zz+hvn(i,j,k)
end do
hvn(i,j,kmax)=ssvn(i,j)-zz
end if
end do
end do
call wait_halo(HU_TAG)
call update_3d_halo(hvn,hvn,av,iimin,jjmin,iimax,jjmax,kmax,HV_TAG)
do j=jjmin-1,jjmax+1
do i=iimin-1,iimax+1
if(az(i,j) .gt. 0) then
r = cord_relax/dt*H(i,j)/maxdepth
zz = -H(i,j)
do k=1,kmax-1
HH=max(ssen(i,j)+H(i,j),min_depth)
hn(i,j,k)=(ho(i,j,k)*r+HH*(gga(i,j,k)-gga(i,j,k-1)))/(r+1.)
zz = zz + hn(i,j,k)
end do
hn(i,j,kmax)=ssen(i,j)-zz
end if
end do
end do
call wait_halo(HV_TAG)
case default
end select
#ifdef DEBUG
write(debug,*) 'Leaving Coordinates()'
write(debug,*)
#endif
return
end subroutine coordinates
!EOC
!-----------------------------------------------------------------------
! Copyright (C) 2001 - Hans Burchard and Karsten Bolding !
!-----------------------------------------------------------------------
!$Id: internal_pressure.F90.SONG,v 1.1 2002-05-02 14:01:01 gotm Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: internal_pressure()
!
! !INTERFACE:
subroutine internal_pressure()
!
! !DESCRIPTION:
!
! !USES:
use domain, only: iimin,iimax,jjmin,jjmax,kmax,az,au,av,H,HU,HV
#ifdef CURVILINEAR
use domain, only: dxu,dyv
#else
use domain, only: dx,dy
#endif
use variables_3d, only: kmin,hn,hun,hvn,idpdx,idpdy,rho
IMPLICIT NONE
!
! !INPUT PARAMETERS:
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !REVISION HISTORY:
! Original author(s): Hans Burchard & Karsten Bolding
!
! $Log: internal_pressure.F90.SONG,v $
! Revision 1.1 2002-05-02 14:01:01 gotm
! Initial revision
!
! Revision 1.15 2001/10/22 08:23:53 bbh
! Removed reference to kplus and kminus when not -DPRESS_GRAD_Z
!
! Revision 1.14 2001/10/22 07:47:59 bbh
! Needed to run dos2unix
!
! Revision 1.13 2001/10/22 07:45:26 bbh
! Fixed a serious bug - now gives the same as old version
!
! Revision 1.10 2001/09/04 07:29:18 bbh
! Internal pressure based on interpolation to z-coordinates
!
! Revision 1.9 2001/09/03 13:03:37 bbh
! Initial pressure gradient can now be subtracted
!
! Revision 1.8 2001/08/31 15:40:37 bbh
! initial pressure can be subtracted now
!
! Revision 1.7 2001/08/01 08:31:22 bbh
! CURVILINEAR now implemented
!
! Revision 1.6 2001/06/22 08:19:10 bbh
! Compiler options such as USE_MASK and OLD_DRY deleted.
! Open and passive boundary for z created.
! Various inconsistencies removed.
! wait_halo added.
! Checked loop boundaries
!
! Revision 1.5 2001/05/25 19:09:48 bbh
! Removed dead code
!
! Revision 1.4 2001/05/20 09:19:09 bbh
! Use who specified twice
!
! Revision 1.3 2001/05/20 07:51:40 bbh
! Internal pressure included
!
! Revision 1.2 2001/05/11 13:47:00 bbh
! Added actual code
!
! Revision 1.1 2001/05/10 11:30:16 bbh
! Added further support for baroclinicity
!
! !LOCAL VARIABLES:
integer :: i,j,k,rc
REALTYPE :: dxm1,dym1,x,y,x1,y1,grdl,grdu,rhol,rhou,prgr,dxz,dyz
REALTYPE,dimension(:,:,:), allocatable :: zz
REALTYPE,save,dimension(:,:,:), allocatable :: idpdx0
REALTYPE,save,dimension(:,:,:), allocatable :: idpdy0
LOGICAL,save :: first=.true.
#ifdef PRESS_GRAD_Z
REALTYPE :: zx(kmax)
REALTYPE :: rhoplus,rhominus
integer :: kplus,kminus
#endif
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
integer, save :: Ncall = 0
Ncall = Ncall+1
write(debug,*) 'internal_pressure() # ',Ncall
#endif
#ifndef CURVILINEAR
dxm1 = _ONE_/DXU
dym1 = _ONE_/DYV
#endif
#ifdef SUBSTR_INI_PRESS
if (first) then
allocate(idpdx0(I3DFIELD),stat=rc) ! Initial x - pressure gradient.
if (rc /= 0) stop 'internal_pressure.F90: Error allocating memory (idpdx0)'
allocate(idpdy0(I3DFIELD),stat=rc) ! Initial y - pressure gradient.
if (rc /= 0) stop 'internal_pressure.F90: Error allocating memory (idpdy0)'
end if
#endif
allocate(zz(I3DFIELD),stat=rc) ! Interface heights.
if (rc /= 0) stop 'internal_pressure.F90: Error allocating memory (zz)'
zz = _ZERO_
idpdx = _ZERO_
idpdy = _ZERO_
! First, the interface heights are calculated in order to get the
! interface slopes further down.
do j=jjmin,jjmax+1
do i=iimin,iimax+1
if (az(i,j) .ge. 1) then
#ifdef PRESS_GRAD_Z
zz(i,j,1)=-H(i,j)+0.5*hn(i,j,1)
do k=2,kmax
zz(i,j,k)=zz(i,j,k-1)+0.5*(hn(i,j,k-1)+hn(i,j,k))
end do
#else
! zz(i,j,1)=-H(i,j)+hn(i,j,1)
! do k=2,kmax
! zz(i,j,k)=zz(i,j,k-1)+hn(i,j,k)
! end do
zz(i,j,1)=-H(i,j)+0.5*hn(i,j,1)
do k=2,kmax
zz(i,j,k)=zz(i,j,k-1)+0.5*(hn(i,j,k-1)+hn(i,j,k))
end do
#endif
end if
end do
end do
! Calculation of layer integrated internal pressure gradient as it
! appears on the right hand side of the u-velocity equation.
do j=jjmin,jjmax
do i=iimin,iimax
#ifdef CURVILINEAR
dxm1=_ONE_/DXU
#endif
if (au(i,j) .ge. 1) then
#ifdef PRESS_GRAD_Z
zx(1)=-HU(i,j)+0.5*hun(i,j,1) ! zx defined on u-points
do k=2,kmax
zx(k)=zx(k-1)+0.5*(hun(i,j,k-1)+hun(i,j,k))
end do
#endif
! grdl=0.5*hun(i,j,kmax)*(rho(i+1,j,kmax)-rho(i,j,kmax))*dxm1 !SONG
grdl=0.5*(rho(i+1,j,kmax)*hn(i+1,j,kmax)-rho(i,j,kmax)*hn(i,j,kmax))*dxm1
rhol=0.5*(rho(i+1,j,kmax)+rho(i,j,kmax))*(zz(i+1,j,kmax)-zz(i,j,kmax))*dxm1 !SONG
prgr=grdl
idpdx(i,j,kmax)=hun(i,j,kmax)*prgr
do k=kmax-1,1,-1
grdu=grdl
#ifdef PRESS_GRAD_Z
do kplus=kmax,1,-1
if (zz(i+1,j,kplus) .le. zx(k)) EXIT
end do
do kminus=kmax,1,-1
if (zz(i,j,kminus) .le. zx(k)) EXIT
end do
if (kplus .eq. kmax) rhoplus=rho(i+1,j,kplus)
if (kplus .lt. kmax .and. kplus .gt. 1) then
rhoplus=((zx(k)-zz(i+1,j,kplus))*rho(i+1,j,kplus+1)+ &
(zz(i+1,j,kplus+1)-zx(k))*rho(i+1,j,kplus))/ &
(0.5*(hn(i+1,j,kplus+1)+hn(i+1,j,kplus)))
end if
if (kminus .eq. kmax) rhominus=rho(i,j,kminus)
if ((kminus .lt. kmax) .and. (kminus .gt. 1)) then
rhominus=((zx(k)-zz(i,j,kminus))*rho(i,j,kminus+1)+ &
(zz(i,j,kminus+1)-zx(k))*rho(i,j,kminus))/ &
(0.5*(hn(i,j,kminus+1)+hn(i,j,kminus)))
end if
if (zx(k) .gt. max(-H(i+1,j),-H(i,j))) then
grdl=0.5*hun(i,j,k)*(rhoplus-rhominus)*dxm1
else
grdl= _ZERO_
end if
prgr=prgr+grdu+grdl
#else
grdl=0.5*hun(i,j,k)*(rho(i+1,j,k)-rho(i,j,k))*dxm1 !SONG
grdl=(0.5*(rho(i+1,j,k)+rho(i+1,j,k+1)) &
*0.5*(hn(i+1,j,k)+hn(i+1,j,k+1)) &
-0.5*(rho(i,j,k)+rho(i,j,k+1)) &
*0.5*(hn(i,j,k)+hn(i,j,k+1)) )*dxm1
rhou=rhol
rhol=0.5*(rho(i+1,j,k)+rho(i,j,k))*(zz(i+1,j,k)-zz(i,j,k))*dxm1
prgr=prgr+grdl-(rhou-rhol) !SONG
#endif
idpdx(i,j,k)=hun(i,j,k)*prgr
end do
end if
end do
end do
! Calculation of layer integrated internal pressure gradient as it
! appears on the right hand side of the v-velocity equation.
do j=jjmin,jjmax
do i=iimin,iimax
#ifdef CURVILINEAR
dym1 = _ONE_/DYV
#endif
if (av(i,j) .ge. 1) then
#ifdef PRESS_GRAD_Z
zx(1)=-HV(i,j)+0.5*hvn(i,j,1) ! zx defined on v-points
do k=2,kmax
zx(k)=zx(k-1)+0.5*(hvn(i,j,k-1)+hvn(i,j,k))
end do
#endif
! grdl=0.5*hvn(i,j,kmax)*(rho(i,j+1,kmax)-rho(i,j,kmax))*dym1
! rhol=0.5*(rho(i,j+1,kmax)+rho(i,j,kmax))
grdl=0.5*(rho(i,j+1,kmax)*hn(i,j+1,kmax)-rho(i,j,kmax)*hn(i,j,kmax))*dym1
rhol=0.5*(rho(i,j+1,kmax)+rho(i,j,kmax))*(zz(i,j+1,kmax)-zz(i,j,kmax))*dxm1 !SONG
prgr=grdl
idpdy(i,j,kmax)=hvn(i,j,kmax)*prgr
do k=kmax-1,1,-1
grdu=grdl
#ifdef PRESS_GRAD_Z
do kplus=kmax,1,-1
if (zz(i,j+1,kplus) .le. zx(k)) EXIT
end do
do kminus=kmax,1,-1
if (zz(i,j,kminus) .le. zx(k)) EXIT
end do
if (kplus .eq. kmax) rhoplus=rho(i,j+1,kplus)
if ((kplus .lt. kmax) .and. (kplus .gt. 1)) then
rhoplus=((zx(k)-zz(i,j+1,kplus))*rho(i,j+1,kplus+1)+ &
(zz(i,j+1,kplus+1)-zx(k))*rho(i,j+1,kplus))/ &
(0.5*(hn(i,j+1,kplus+1)+hn(i,j+1,kplus)))
end if
if (kminus .eq. kmax) rhominus=rho(i,j,kminus)
if ((kminus .lt. kmax) .and. (kminus .gt. 1)) then
rhominus=((zx(k)-zz(i,j,kminus))*rho(i,j,kminus+1)+ &
(zz(i,j,kminus+1)-zx(k))*rho(i,j,kminus))/ &
(0.5*(hn(i,j,kminus+1)+hn(i,j,kminus)))
end if
if (zx(k).gt.max(-H(i,j+1),-H(i,j))) then
grdl=0.5*hvn(i,j,k)*(rhoplus-rhominus)*dym1
else
grdl= _ZERO_
end if
prgr=prgr+grdu+grdl
#else
! grdl=0.5*hvn(i,j,k)*(rho(i,j+1,k)-rho(i,j,k))*dym1
! rhou=rhol
! rhol=0.5*(rho(i,j+1,k)+rho(i,j,k))
! dyz=(zz(i,j+1,k)-zz(i,j,k))*dym1
! prgr=prgr+grdu+grdl-dyz*(rhou-rhol)
grdl=(0.5*(rho(i,j+1,k)+rho(i,j+1,k+1)) &
*0.5*(hn(i,j+1,k)+hn(i,j+1,k+1)) &
-0.5*(rho(i,j,k)+rho(i,j,k+1)) &
*0.5*(hn(i,j,k)+hn(i,j,k+1)) )*dym1
rhou=rhol
rhol=0.5*(rho(i,j+1,k)+rho(i,j,k))*(zz(i,j+1,k)-zz(i,j,k))*dym1
prgr=prgr+grdl-(rhou-rhol) !SONG
#endif
idpdy(i,j,k)=hvn(i,j,k)*prgr
end do
end if
end do
end do
#ifdef SUBSTR_INI_PRESS
if (first) then
first = .false.
do k=0,kmax
do j=jjmin,jjmax
do i=iimin,iimax
idpdx0(i,j,k) = idpdx(i,j,k)
idpdx(i,j,k) = _ZERO_
idpdy0(i,j,k) = idpdy(i,j,k)
idpdy(i,j,k) = _ZERO_
end do
end do
end do
else
do k=0,kmax
do j=jjmin,jjmax
do i=iimin,iimax
idpdx(i,j,k) = idpdx(i,j,k) - idpdx0(i,j,k)
idpdy(i,j,k) = idpdy(i,j,k) - idpdy0(i,j,k)
end do
end do
end do
end if
#endif
#ifdef FORTRAN90
deallocate(zz,stat=rc) ! work array
if (rc /= 0) stop 'internal_pressure: Error de-allocating memory (zz)'
#endif
#ifdef DEBUG
write(debug,*) 'Leaving internal_pressure()'
write(debug,*)
#endif
return
end subroutine internal_pressure
!EOC
!-----------------------------------------------------------------------
! Copyright (C) 2001 - Hans Burchard and Karsten Bolding !
!-----------------------------------------------------------------------
!$Id: internal_pressure.F90.new,v 1.1 2002-05-02 14:01:00 gotm Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: internal_pressure()
!