Skip to content

GitLab

G
getm
Commit 56b1152a authored by kbk's avatar kbk
Browse files
Options

introduce buoy() and rename rho() to buoy()i where appropriate

parent c44c0293
Hide whitespace changes
Inline Side-by-side
Showing with 134 additions and 177 deletions
doc/3dIntro.tex
View file @ 56b1152a
%$Id: 3dIntro.tex,v 1.2 2007-01-10 22:01:16 hb Exp $
%$Id: 3dIntro.tex,v 1.3 2007-02-23 12:20:35 kbk Exp $
\section{Introduction to 3d module}
......@@ -329,5 +329,6 @@ will be added as an option in the near future.
For the equation of state, also linearised version are implemented into
GETM, for details, see section \ref{sec-eqstate} on page \pageref{sec-eqstate}.
It should be noted carefully that the GETM variable {\tt rho} actually
stands for the buoyancy $b$ as defined in (\ref{bdef}).
For convinient use in other subroutines the buoyancy $b$ as defined in
(\ref{bdef} is calculated and stored in the GETM variable {tt buoy}.
src/3d/dynamic_allocations_3d.h
View file @ 56b1152a
......@@ -59,9 +59,12 @@
allocate(T(I3DFIELD),stat=rc) ! 3D field for temperature
if (rc /= 0) stop 'init_3d: Error allocating memory (T)'
allocate(rho(I3DFIELD),stat=rc) ! 3D field for density (actual buoyancy)
allocate(rho(I3DFIELD),stat=rc) ! 3D field for density
if (rc /= 0) stop 'init_3d: Error allocating memory (rho)'
allocate(buoy(I3DFIELD),stat=rc) ! 3D field for buoyancy
if (rc /= 0) stop 'init_3d: Error allocating memory (buoy)'
allocate(idpdx(I3DFIELD),stat=rc) ! Internal pressure gradient - x
if (rc /= 0) stop 'init_3d: Error allocating memory (idpdx)'
......
src/3d/dynamic_declarations_3d.h
View file @ 56b1152a
! Remember to update this value if you add more 3D arrays.
#ifdef UV_TVD
integer,parameter :: n3d_fields=30
integer,parameter :: n3d_fields=33
#else
integer,parameter :: n3d_fields=23
integer,parameter :: n3d_fields=26
#endif
! Number of vertical layers in z,u,v columns
......@@ -23,9 +23,9 @@
REALTYPE, dimension(:,:,:), allocatable :: NN
! 3D baroclinic fields
REALTYPE, dimension(:,:,:), allocatable :: S,T,rho
REALTYPE, dimension(:,:,:), allocatable :: S,T,rho,buoy
REALTYPE, dimension(:,:,:), allocatable :: idpdx,idpdy
REALTYPE, dimension(:,:,:), allocatable :: rad
REALTYPE, dimension(:,:,:), allocatable :: rad,light
#endif
! suspended matter
......@@ -34,8 +34,6 @@
REALTYPE, dimension(:,:), allocatable :: spm_pool
#endif
REALTYPE, dimension(:,:,:), allocatable :: light
#ifdef UV_TVD
REALTYPE, dimension(:,:,:), allocatable :: uadv,vadv,wadv
REALTYPE, dimension(:,:,:), allocatable :: huadv,hvadv,hoadv,hnadv
......
src/3d/eqstate.F90
View file @ 56b1152a
!$Id: eqstate.F90,v 1.8 2006-03-23 09:53:43 frv-bjb Exp $
!$Id: eqstate.F90,v 1.9 2007-02-23 12:20:36 kbk Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -104,7 +104,7 @@
!
! !USES:
use domain, only: iimin,iimax,jjmin,jjmax,kmax,az
use variables_3d, only: T,S,rho
use variables_3d, only: T,S,rho,buoy
IMPLICIT NONE
!
! !INPUT PARAMETERS:
......@@ -116,21 +116,6 @@
! !LOCAL VARIABLES:
integer :: i,j,k
REALTYPE :: x
#if 0
! Brydon et. al. - Table 2 - narrow
REALTYPE, parameter :: a1=-1.36471e-1,b1= 5.06423e-1,g1=-5.52640e-4
REALTYPE, parameter :: a2=-4.68181e-2,b2=-3.57109e-3,g2= 4.88584e-6
REALTYPE, parameter :: a3= 8.07004e-1,b3=-8.76148e-4,g3= 9.96027e-7
REALTYPE, parameter :: a4=-7.45353e-3,b4= 5.25243e-5,g4=-7.25139e-8
REALTYPE, parameter :: a5=-2.94418e-3,b5= 1.57976e-5,g5=-3.98736e-9
REALTYPE, parameter :: a6= 3.43570e-5,b6=-3.46686e-7,g6= 4.00631e-10
REALTYPE, parameter :: a7= 3.48658e-5,b7=-1.68764e-7,g7= 8.26368e-11
!
REALTYPE :: p=50.,p2
REALTYPE :: c1,c2,c3,c4,c5,c6,c7
!
REALTYPE :: s1,t1,t2,t3
#endif
REALTYPE :: KK
REALTYPE :: T1,T2,T3,T4,T5,S1,S15,S2,S3,p2
!EOP
......@@ -145,27 +130,17 @@
#define BUOYANCY
select case (eqstate_method)
case (1)
#ifdef DENSITY
forall(i=iimin-1:iimax+1,j=jjmin-1:jjmax+1,az(i,j) .gt. 0) &
rho(i,j,1:kmax) = rho_0 + &
dtr0*(T(i,j,1:kmax)-T0) + dsr0*(S(i,j,1:kmax)-S0)
#endif
#undef DENSITY
#ifdef BUOYANCY
x = -g/rho_0
forall(i=iimin-1:iimax+1,j=jjmin-1:jjmax+1,az(i,j) .gt. 0) &
rho(i,j,1:kmax) = &
x*(dtr0*(T(i,j,1:kmax)-T0) + dsr0*(S(i,j,1:kmax)-S0))
#endif
#ifdef HAIDVOGEL_TEST
forall(i=iimin-1:iimax+1,j=jjmin-1:jjmax+1,az(i,j) .gt. 0) &
rho(i,j,1:kmax) = -g/rho_0*(S(i,j,1:kmax)+1000.-rho_0)
rho(i,j,1:kmax) = 1000. + S(i,j,1:kmax)
#endif
#ifdef CONSTANCE_TEST
x = -g/rho_0
forall(i=iimin-1:iimax+1,j=jjmin-1:jjmax+1,az(i,j) .gt. 0) &
rho(i,j,1:kmax) = x*dtr0*(T(i,j,1:kmax)-T0)
rho(i,j,1:kmax) = 1000. + *dtr0*(T(i,j,1:kmax)-T0)
#endif
case (2)
do k = 1,kmax
......@@ -217,41 +192,19 @@
x=x/(1.-p/KK)
end if
#endif
rho(i,j,k)=-g*(x-rho_0)/rho_0
rho(i,j,k)=x
end if
end do
end do
end do
#if 0
do k = 1,kmax
do j = jjmin-1:jjmax+1
do i = iimin-1:iimax+1
if (az(i,j) .gt. 0) then
p = 50.
p2 = p*p
c1 = a1 + b1*p + g1*p2
c2 = a2 + b2*p + g2*p2
c3 = a3 + b3*p + g3*p2
c4 = a4 + b4*p + g4*p2
c5 = a5 + b5*p + g5*p2
c6 = a6 + b6*p + g6*p2
c7 = a7 + b7*p + g7*p2
s1 = S(i,j,k)
t1 = T(i,j,k)
t2 = t1*t1
t3 = t1*t2
rho(i,j,k)=c1+c2*t1+c3*s1+c4*t2+c5*s1*t1+c6*t3+c7*s1*t2
rho(i,j,k)=-g/rho_0*(rho(i,j,k)-rho_0)
end if
end do
end do
end do
#endif
case default
end select
#undef BUOYANCY
x=-g/rho_0
buoy=x*(rho-rho_0)
#ifdef DEBUG
write(debug,*) 'leaving do_eqstate()'
write(debug,*)
......
src/3d/internal_pressure.F90
View file @ 56b1152a
!$Id: internal_pressure.F90,v 1.16 2006-03-01 15:54:08 kbk Exp $
!$Id: internal_pressure.F90,v 1.17 2007-02-23 12:20:36 kbk Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -64,7 +64,7 @@
#else
use domain, only: dx,dy
#endif
use variables_3d, only: kmin,hn,hun,hvn,idpdx,idpdy,rho
use variables_3d, only: kmin,hn,hun,hvn,idpdx,idpdy,buoy
IMPLICIT NONE
!
! !PUBLIC DATA MEMBERS:
......
src/3d/ip_blumberg_mellor.F90
View file @ 56b1152a
!$Id: ip_blumberg_mellor.F90,v 1.7 2006-03-01 15:54:08 kbk Exp $
!$Id: ip_blumberg_mellor.F90,v 1.8 2007-02-23 12:20:36 kbk Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -59,7 +59,7 @@
! !LOCAL VARIABLES:
integer :: i,j,k
REALTYPE :: dxm1,dym1
REALTYPE :: grdl,grdu,rhol,rhou,prgr,dxz,dyz
REALTYPE :: grdl,grdu,buoyl,buoyu,prgr,dxz,dyz
!EOP
!-----------------------------------------------------------------------
!BOC
......@@ -89,17 +89,17 @@
#if defined(SPHERICAL) || defined(CURVILINEAR)
dxm1=_ONE_/DXU
#endif
grdl=(rho(i+1,j,kmax)-rho(i,j,kmax))*dxm1
rhol=0.5*(rho(i+1,j,kmax)+rho(i,j,kmax))
grdl=(buoy(i+1,j,kmax)-buoy(i,j,kmax))*dxm1
buoyl=0.5*(buoy(i+1,j,kmax)+buoy(i,j,kmax))
prgr=grdl
idpdx(i,j,kmax)=hun(i,j,kmax)*prgr*0.5*hun(i,j,kmax)
do k=kmax-1,1,-1
grdu=grdl
grdl=(rho(i+1,j,k)-rho(i,j,k))*dxm1
rhou=rhol
rhol=0.5*(rho(i+1,j,k)+rho(i,j,k))
grdl=(buoy(i+1,j,k)-buoy(i,j,k))*dxm1
buoyu=buoyl
buoyl=0.5*(buoy(i+1,j,k)+buoy(i,j,k))
dxz=(zz(i+1,j,k)-zz(i,j,k))*dxm1
prgr=prgr+0.5*(grdu+grdl)*0.5*(hun(i,j,k)+hun(i,j,k+1))-dxz*(rhou-rhol)
prgr=prgr+0.5*(grdu+grdl)*0.5*(hun(i,j,k)+hun(i,j,k+1))-dxz*(buoyu-buoyl)
idpdx(i,j,k)=hun(i,j,k)*prgr
end do
end if
......@@ -114,17 +114,17 @@
#if defined(SPHERICAL) || defined(CURVILINEAR)
dym1 = _ONE_/DYV
#endif
grdl=(rho(i,j+1,kmax)-rho(i,j,kmax))*dym1
rhol=0.5*(rho(i,j+1,kmax)+rho(i,j,kmax))
grdl=(buoy(i,j+1,kmax)-buoy(i,j,kmax))*dym1
buoyl=0.5*(buoy(i,j+1,kmax)+buoy(i,j,kmax))
prgr=grdl
idpdy(i,j,kmax)=hvn(i,j,kmax)*prgr*0.5*hvn(i,j,kmax)
do k=kmax-1,1,-1
grdu=grdl
grdl=(rho(i,j+1,k)-rho(i,j,k))*dym1
rhou=rhol
rhol=0.5*(rho(i,j+1,k)+rho(i,j,k))
grdl=(buoy(i,j+1,k)-buoy(i,j,k))*dym1
buoyu=buoyl
buoyl=0.5*(buoy(i,j+1,k)+buoy(i,j,k))
dyz=(zz(i,j+1,k)-zz(i,j,k))*dym1
prgr=prgr+0.5*(grdu+grdl)*0.5*(hvn(i,j,k)+hvn(i,j,k+1))-dyz*(rhou-rhol)
prgr=prgr+0.5*(grdu+grdl)*0.5*(hvn(i,j,k)+hvn(i,j,k+1))-dyz*(buoyu-buoyl)
idpdy(i,j,k)=hvn(i,j,k)*prgr
end do
end if
......
src/3d/ip_blumberg_mellor_lin.F90
View file @ 56b1152a
!$Id: ip_blumberg_mellor_lin.F90,v 1.4 2006-03-01 15:54:08 kbk Exp $
!$Id: ip_blumberg_mellor_lin.F90,v 1.5 2007-02-23 12:20:36 kbk Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -92,15 +92,15 @@
dxm1=_ONE_/DXU
#endif
dxzl=(zz(i+1,j,kmax)-zz(i,j,kmax))*dxm1
dxrl=(rho(i+1,j,kmax)-rho(i,j,kmax))*dxm1
dxrl=(buoy(i+1,j,kmax)-buoy(i,j,kmax))*dxm1
prgr=0.
do k=kmax-1,kumin_pmz(i,j),-1
dxzu=dxzl
dxzl=(zz(i+1,j,k)-zz(i,j,k))*dxm1
dxru=dxrl
dxrl=(rho(i+1,j,k)-rho(i,j,k))*dxm1
dzr2=(rho(i+1,j,k+1)-rho(i+1,j,k))/(zz(i+1,j,k+1)-zz(i+1,j,k))
dzr1=(rho(i ,j,k+1)-rho(i ,j,k))/(zz(i ,j,k+1)-zz(i ,j,k))
dxrl=(buoy(i+1,j,k)-buoy(i,j,k))*dxm1
dzr2=(buoy(i+1,j,k+1)-buoy(i+1,j,k))/(zz(i+1,j,k+1)-zz(i+1,j,k))
dzr1=(buoy(i ,j,k+1)-buoy(i ,j,k))/(zz(i ,j,k+1)-zz(i ,j,k))
aa=0.5*(dxrl+dxru)
bb=0.5*(dxzl+dxzu)
cc=0.5*(dzr2+dzr1)
......@@ -121,16 +121,16 @@
dym1 = _ONE_/DYV
#endif
dyzl=(zz(i,j+1,kmax)-zz(i,j,kmax))*dym1
dyrl=(rho(i,j+1,kmax)-rho(i,j,kmax))*dym1
dyrl=(buoy(i,j+1,kmax)-buoy(i,j,kmax))*dym1
prgr=0.
idpdy(i,j,kmax)=hvn(i,j,kmax)*prgr
do k=kmax-1,kvmin_pmz(i,j),-1
dyzu=dyzl
dyzl=(zz(i,j+1,k)-zz(i,j,k))*dym1
dyru=dyrl
dyrl=(rho(i,j+1,k)-rho(i,j,k))*dym1
dzr2=(rho(i,j+1,k+1)-rho(i,j+1,k))/(zz(i,j+1,k+1)-zz(i,j+1,k))
dzr1=(rho(i,j ,k+1)-rho(i ,j,k))/(zz(i ,j,k+1)-zz(i ,j,k))
dyrl=(buoy(i,j+1,k)-buoy(i,j,k))*dym1
dzr2=(buoy(i,j+1,k+1)-buoy(i,j+1,k))/(zz(i,j+1,k+1)-zz(i,j+1,k))
dzr1=(buoy(i,j ,k+1)-buoy(i ,j,k))/(zz(i ,j,k+1)-zz(i ,j,k))
aa=0.5*(dyrl+dyru)
bb=0.5*(dyzl+dyzu)
cc=0.5*(dzr2+dzr1)
......
src/3d/ip_chu_fan.F90
View file @ 56b1152a
!$Id: ip_chu_fan.F90,v 1.4 2006-03-01 15:54:08 kbk Exp $
!$Id: ip_chu_fan.F90,v 1.5 2007-02-23 12:20:37 kbk Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -24,7 +24,7 @@
! !LOCAL VARIABLES:
integer :: i,j,k
REALTYPE :: dxm1,dym1,x,y,x1,y1,hc
REALTYPE :: grdl,grdu,rhol,rhou,prgr,dxz,dyz
REALTYPE :: grdl,grdu,buoyl,buoyu,prgr,dxz,dyz
!EOP
!-----------------------------------------------------------------------
!BOC
......@@ -54,22 +54,22 @@
#if defined(SPHERICAL) || defined(CURVILINEAR)
dxm1=_ONE_/DXU
#endif
grdl=(rho(i+1,j,kmax)-rho(i,j,kmax))*dxm1
rhol=0.5*(rho(i+1,j,kmax)+rho(i,j,kmax))
grdl=(buoy(i+1,j,kmax)-buoy(i,j,kmax))*dxm1
buoyl=0.5*(buoy(i+1,j,kmax)+buoy(i,j,kmax))
prgr=grdl
idpdx(i,j,kmax)=hun(i,j,kmax)*prgr*0.5*hun(i,j,kmax)
do k=kmax-1,1,-1
grdu=grdl
grdl=(rho(i+1,j,k)-rho(i,j,k))*dxm1
rhou=rhol
rhol=0.5*(rho(i+1,j,k)+rho(i,j,k))
grdl=(buoy(i+1,j,k)-buoy(i,j,k))*dxm1
buoyu=buoyl
buoyl=0.5*(buoy(i+1,j,k)+buoy(i,j,k))
dxz=(zz(i+1,j,k)-zz(i,j,k))*dxm1
hc=(zz(i+1,j,k)-zz(i+1,j,k+1))**2*(rho(i+1,j,k)-rho(i+1,j,k+1))
hc=hc-(zz(i,j,k)-zz(i,j,k+1))**2*(rho(i,j,k)-rho(i,j,k+1))
hc=hc+(zz(i+1,j,k+1)-zz(i,j,k+1))**2*(rho(i+1,j,k+1)-rho(i,j,k+1))
hc=hc-(zz(i+1,j,k)-zz(i,j,k))**2*(rho(i+1,j,k)-rho(i,j,k))
hc=(zz(i+1,j,k)-zz(i+1,j,k+1))**2*(buoy(i+1,j,k)-buoy(i+1,j,k+1))
hc=hc-(zz(i,j,k)-zz(i,j,k+1))**2*(buoy(i,j,k)-buoy(i,j,k+1))
hc=hc+(zz(i+1,j,k+1)-zz(i,j,k+1))**2*(buoy(i+1,j,k+1)-buoy(i,j,k+1))
hc=hc-(zz(i+1,j,k)-zz(i,j,k))**2*(buoy(i+1,j,k)-buoy(i,j,k))
hc=hc*dxm1*0.1666666667/(zz(i,j,k)+zz(i+1,j,k)-zz(i,j,k+1)-zz(i+1,j,k+1))
prgr=prgr+(grdu+grdl+hc)*0.5*(hun(i,j,k)+hun(i,j,k+1))-dxz*(rhou-rhol)
prgr=prgr+(grdu+grdl+hc)*0.5*(hun(i,j,k)+hun(i,j,k+1))-dxz*(buoyu-buoyl)
idpdx(i,j,k)=hun(i,j,k)*prgr
end do
end if
......@@ -84,22 +84,22 @@
#if defined(SPHERICAL) || defined(CURVILINEAR)
dym1 = _ONE_/DYV
#endif
grdl=(rho(i,j+1,kmax)-rho(i,j,kmax))*dym1
rhol=0.5*(rho(i,j+1,kmax)+rho(i,j,kmax))
grdl=(buoy(i,j+1,kmax)-buoy(i,j,kmax))*dym1
buoyl=0.5*(buoy(i,j+1,kmax)+buoy(i,j,kmax))
prgr=grdl
idpdy(i,j,kmax)=hvn(i,j,kmax)*prgr*0.5*hvn(i,j,kmax)
do k=kmax-1,1,-1
grdu=grdl
grdl=(rho(i,j+1,k)-rho(i,j,k))*dym1
rhou=rhol
rhol=0.5*(rho(i,j+1,k)+rho(i,j,k))
grdl=(buoy(i,j+1,k)-buoy(i,j,k))*dym1
buoyu=buoyl
buoyl=0.5*(buoy(i,j+1,k)+buoy(i,j,k))
dyz=(zz(i,j+1,k)-zz(i,j,k))*dym1
hc=(zz(i,j+1,k)-zz(i,j+1,k+1))**2*(rho(i,j+1,k)-rho(i,j+1,k+1))
hc=hc-(zz(i,j,k)-zz(i,j,k+1))**2*(rho(i,j,k)-rho(i,j,k+1))
hc=hc+(zz(i,j+1,k+1)-zz(i,j,k+1))**2*(rho(i,j+1,k+1)-rho(i,j,k+1))
hc=hc-(zz(i,j+1,k)-zz(i,j,k))**2*(rho(i,j+1,k)-rho(i,j,k))
hc=(zz(i,j+1,k)-zz(i,j+1,k+1))**2*(buoy(i,j+1,k)-buoy(i,j+1,k+1))
hc=hc-(zz(i,j,k)-zz(i,j,k+1))**2*(buoy(i,j,k)-buoy(i,j,k+1))
hc=hc+(zz(i,j+1,k+1)-zz(i,j,k+1))**2*(buoy(i,j+1,k+1)-buoy(i,j,k+1))
hc=hc-(zz(i,j+1,k)-zz(i,j,k))**2*(buoy(i,j+1,k)-buoy(i,j,k))
hc=hc*dym1*0.1666666667/(zz(i,j,k)+zz(i,j+1,k)-zz(i,j,k+1)-zz(i,j+1,k+1))
prgr=prgr+(grdu+grdl+hc)*0.5*(hvn(i,j,k)+hvn(i,j,k+1))-dyz*(rhou-rhol)
prgr=prgr+(grdu+grdl+hc)*0.5*(hvn(i,j,k)+hvn(i,j,k+1))-dyz*(buoyu-buoyl)
idpdy(i,j,k)=hvn(i,j,k)*prgr
end do
end if
......
src/3d/ip_song_wright.F90
View file @ 56b1152a
!$Id: ip_song_wright.F90,v 1.4 2006-03-01 15:54:08 kbk Exp $
!$Id: ip_song_wright.F90,v 1.5 2007-02-23 12:20:37 kbk Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -52,7 +52,7 @@
! !LOCAL VARIABLES:
integer :: i,j,k
REALTYPE :: dxm1,dym1
REALTYPE :: grdl,grdu,rhol,rhou,prgr,dxz,dyz
REALTYPE :: grdl,grdu,buoyl,buoyu,prgr,dxz,dyz
!EOP
!-----------------------------------------------------------------------
!BOC
......@@ -82,20 +82,20 @@
#if defined(SPHERICAL) || defined(CURVILINEAR)
dxm1=_ONE_/DXU
#endif
grdl=0.5*hun(i,j,kmax)*(rho(i+1,j,kmax)-rho(i,j,kmax))*dxm1
rhol=0.5*(rho(i+1,j,kmax)+rho(i,j,kmax)) &
grdl=0.5*hun(i,j,kmax)*(buoy(i+1,j,kmax)-buoy(i,j,kmax))*dxm1
buoyl=0.5*(buoy(i+1,j,kmax)+buoy(i,j,kmax)) &
*(zz(i+1,j,kmax)- zz(i,j,kmax))*dxm1
prgr=grdl
idpdx(i,j,kmax)=hun(i,j,kmax)*prgr
do k=kmax-1,1,-1
grdu=grdl
grdl=(0.5*(rho(i+1,j,k)+rho(i+1,j,k+1)) &
grdl=(0.5*(buoy(i+1,j,k)+buoy(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*(buoy(i,j,k)+buoy(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)
buoyu=buoyl
buoyl=0.5*(buoy(i+1,j,k)+buoy(i,j,k))*(zz(i+1,j,k)-zz(i,j,k))*dxm1
prgr=prgr+grdl-(buoyu-buoyl)
idpdx(i,j,k)=hun(i,j,k)*prgr
end do
end if
......@@ -110,20 +110,20 @@
#if defined(SPHERICAL) || defined(CURVILINEAR)
dym1 = _ONE_/DYV
#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*hvn(i,j,kmax)*(buoy(i,j+1,kmax)-buoy(i,j,kmax))*dym1
buoyl=0.5*(buoy(i,j+1,kmax)+buoy(i,j,kmax)) &
*(zz(i,j+1,kmax)- zz(i,j,kmax))*dxm1
prgr=grdl
idpdy(i,j,kmax)=hvn(i,j,kmax)*prgr
do k=kmax-1,1,-1
grdu=grdl
grdl=(0.5*(rho(i,j+1,k)+rho(i,j+1,k+1)) &
grdl=(0.5*(buoy(i,j+1,k)+buoy(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*(buoy(i,j,k)+buoy(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)
buoyu=buoyl
buoyl=0.5*(buoy(i,j+1,k)+buoy(i,j,k))*(zz(i,j+1,k)-zz(i,j,k))*dym1
prgr=prgr+grdl-(buoyu-buoyl)
idpdy(i,j,k)=hvn(i,j,k)*prgr
end do
end if
......
src/3d/ip_z_interpol.F90
View file @ 56b1152a
!$Id: ip_z_interpol.F90,v 1.4 2006-03-01 15:54:08 kbk Exp $
!$Id: ip_z_interpol.F90,v 1.5 2007-02-23 12:20:37 kbk Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -29,10 +29,10 @@
! !LOCAL VARIABLES:
integer :: i,j,k
REALTYPE :: dxm1,dym1
REALTYPE :: grdl,grdu,rhol,rhou,prgr,dxz,dyz
REALTYPE :: grdl,grdu,buoyl,prgr,dxz,dyz
integer :: kplus,kminus
REALTYPE :: zx(kmax)
REALTYPE :: rhoplus,rhominus
REALTYPE :: buoyplus,buoyminus
!EOP
!-----------------------------------------------------------------------
!BOC
......@@ -65,8 +65,8 @@
do k=2,kmax
zx(k)=zx(k-1)+0.5*(hun(i,j,k-1)+hun(i,j,k))
end do
grdl=0.5*hun(i,j,kmax)*(rho(i+1,j,kmax)-rho(i,j,kmax))*dxm1
rhol=0.5*(rho(i+1,j,kmax)+rho(i,j,kmax))
grdl=0.5*hun(i,j,kmax)*(buoy(i+1,j,kmax)-buoy(i,j,kmax))*dxm1
buoyl=0.5*(buoy(i+1,j,kmax)+buoy(i,j,kmax))
prgr=grdl
idpdx(i,j,kmax)=hun(i,j,kmax)*prgr
do k=kmax-1,1,-1
......@@ -77,20 +77,20 @@
do kminus=kmax,1,-1 ! Find neighboring index to west
if (zz(i,j,kminus) .le. zx(k)) EXIT
end do
if (kplus .eq. kmax) rhoplus=rho(i+1,j,kplus)
if (kplus .eq. kmax) buoyplus=buoy(i+1,j,kplus)
if (kplus .lt. kmax .and. kplus .gt. 1) then ! interpolate
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))/ &
buoyplus=((zx(k)-zz(i+1,j,kplus))*buoy(i+1,j,kplus+1)+ &
(zz(i+1,j,kplus+1)-zx(k))*buoy(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 .eq. kmax) buoyminus=buoy(i,j,kminus)
if ((kminus .lt. kmax) .and. (kminus .gt. 1)) then ! interpolate
rhominus=((zx(k)-zz(i,j,kminus))*rho(i,j,kminus+1)+ &
(zz(i,j,kminus+1)-zx(k))*rho(i,j,kminus))/ &
buoyminus=((zx(k)-zz(i,j,kminus))*buoy(i,j,kminus+1)+ &
(zz(i,j,kminus+1)-zx(k))*buoy(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
grdl=0.5*hun(i,j,k)*(buoyplus-buoyminus)*dxm1
else
grdl= _ZERO_
end if
......@@ -113,8 +113,8 @@
do k=2,kmax
zx(k)=zx(k-1)+0.5*(hvn(i,j,k-1)+hvn(i,j,k))
end do
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*hvn(i,j,kmax)*(buoy(i,j+1,kmax)-buoy(i,j,kmax))*dym1
buoyl=0.5*(buoy(i,j+1,kmax)+buoy(i,j,kmax))
prgr=grdl
idpdy(i,j,kmax)=hvn(i,j,kmax)*prgr
do k=kmax-1,1,-1
......@@ -125,20 +125,20 @@
do kminus=kmax,1,-1 ! Find neighboring index to south
if (zz(i,j,kminus) .le. zx(k)) EXIT
end do
if (kplus .eq. kmax) rhoplus=rho(i,j+1,kplus)
if (kplus .eq. kmax) buoyplus=buoy(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))/ &
buoyplus=((zx(k)-zz(i,j+1,kplus))*buoy(i,j+1,kplus+1)+ &
(zz(i,j+1,kplus+1)-zx(k))*buoy(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 .eq. kmax) buoyminus=buoy(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))/ &
buoyminus=((zx(k)-zz(i,j,kminus))*buoy(i,j,kminus+1)+ &
(zz(i,j,kminus+1)-zx(k))*buoy(i,j,kminus))/ &
(0.5*(hn(i,j,kminus+1)+hn(i,j,kminus)))