copyright (C) 2001 MSC-RPN COMM %%%MC2%%%
*
subroutine rhs3_p0m ( dimphy, wrk , lminx, lmaxx, lminy, lmaxy ) 1
implicit none
*
integer dimphy,lminx,lmaxx,lminy,lmaxy
real*8 wrk(*)
*
*AUTHORs C. Girard & M. Desgagne
*
*OBJECT
*
* *******************************************************************
* * *
* * CALCULATES *
* * right-hand-side LINEAR Pm-terms: *
* * Pu, Pv, Pw, Pb, Pp *
* * FROM *
* * time-levels t-Dt and t variables: *
* * um, vm, wm, bm, pm *
* * u0, v0 , p0 *
* * STORES *
* * results in *
* * up=Pu, vp=Pv, wp=Pw, bp=Pb, pp=Pp *
* * *
* * PERFORMS essentially the following calculations: *
* * *
* * ~ *
* * up = um - dt0 * d(p0)/dX *
* * *
* * ~ *
* * vp = vm - dt0 * d(p0)/dY *
* * *
* * ~ *
* * wp = c03 * wm - dtm * [ d(pm)/dZ - bm / (1+alfa) ] *
* * *
* * __Z *
* * bp = bm - gama_star * pm - dtm*c06*n02g * wm *
* * *
* * ~ __Z *
* * pp = pm/c2_star - dtm*[ d(wm)/dZ-gc02*wm ] - dt0*hdiv(u0,v0) *
* * *
* * *
* * *
* *******************************************************************
*
*
#include "grd.cdk"
#include "consdyn_8.cdk"
#include "yomdyn.cdk"
#include "yomdyn1.cdk"
#include "dynpar.cdk"
#include "dynmem.cdk"
#include "topo.cdk"
#include "phymem.cdk"
#include "physnml.cdk"
#include "physcom.cdk"
#include "nbcpu.cdk"
#include "partopo.cdk"
#include "dtmdtp.cdk"
#include "alfmem.cdk"
*
integer i, j, k, n, id, jd, iff, jf, km1, kp1
real*8 p25, pt5, one, dt0, c01k
real s, d
pointer (pas ,s(*)), (pad ,d(*))
real*8 q3,div
pointer (paq3 , q3(minx:maxx,miny:maxy,*)),
$ (padiv, div(minx:maxx,miny:maxy,*))
parameter (p25=0.25d0, pt5=0.5d0, one=1.0d0)
*----------------------------------------------------------------------
*
paq3 = loc(wrk( 1))
padiv = loc(wrk(dim3d+1))
*
dt0 = dtm+dtp
*
id =1-hx*west
jd =1-hy*south
iff=ldni+(hx-1)*east
jf =ldnj+(hy-1)*north
*
do k = gnk, 1, -1
if (k.eq.gnk) then
do j = jd, jf
do i = id, iff
q3(i,j,gnk) = 0.0
div(i,j,gnk) = 0.0
end do
end do
else
km1 = max( k-1 , 1 )
do j = jd, jf
do i = id, iff
q3(i,j,k) = sbxy(i,j) * pt5 *
$ (-gg1(i+1,j,k) * ( uu0(i+1,j,k) - uu0(i+1,j,km1) )
$ - gg1(i ,j,k) * ( uu0(i ,j,k) - uu0(i ,j,km1) )
$ - gg2(i,j+1,k) * ( vv0(i,j+1,k) - vv0(i,j+1,km1) )
$ - gg2(i,j ,k) * ( vv0(i,j ,k) - vv0(i,j ,km1) ) )
div(i,j,k) = sbxy(i,j)* ((uu0(i+1,j,k) - uu0(i,j,k)) * odx(1) +
$ (vv0(i,j+1,k) - vv0(i,j,k)) * ody(j) )
$ + gg0r(i,j,k) * ( q3(i,j,k+1) + q3(i,j,k)) * pt5
end do
end do
endif
end do
*
do k = gnk, 1, -1
km1 = max( k-1 , 1 )
if (k.lt.gnk) then
kp1 = min( k+1 , gnk-1 )
do j = jd, jf
do i = id+west, iff
uup(i,j,k) = uum(i,j,k) - dt0 * (
$ ( pp0(i,j,k ) - pp0(i-1,j,k) ) * odxu(1) +
$ g0ur(i,j,k) * p25 * (-gg1(i,j,k+1) *
$ ( pp0(i,j,kp1) + pp0(i-1,j,kp1)
$ - pp0(i,j,k ) - pp0(i-1,j,k ) ) -
$ gg1(i,j,k ) *
$ ( pp0(i,j,k ) + pp0(i-1,j,k )
$ - pp0(i,j,km1) - pp0(i-1,j,km1) )) )
end do
end do
do j = jd+south, jf
do i = id, iff
vvp(i,j,k) = vvm(i,j,k) - dt0 * (
$ ( pp0(i,j,k ) - pp0(i,j-1,k) ) * odyv(j) +
$ g0vr(i,j,k) * p25 * (-gg2(i,j,k+1) *
$ ( pp0(i,j,kp1) + pp0(i,j-1,kp1)
$ - pp0(i,j,k ) - pp0(i,j-1,k ) ) -
$ gg2(i,j,k ) *
$ ( pp0(i,j,k ) + pp0(i,j-1,k )
$ - pp0(i,j,km1) - pp0(i,j-1,km1) )) )
end do
end do
end if
if (k.eq.gnk) then
c01k=one
else
c01k=c01
endif
kp1 = min( k+1 , gnk )
do j = jd, jf
do i = id, iff
wwp(i,j,k) = c03*wwm(i,j,k) - dtm *
$ ( (ppm(i,j,k) - ppm(i,j,k-1) ) * g0wr(i,j,k)
$ - bbm(i,j,k) * ap1r(i,j,k) )
bbp(i,j,k) = bbm(i,j,k)
$ - gama_star*pt5*(ppm(i,j,k)+ppm(i,j,k-1))
$ - dtm*c06*n02g(i,j,k)*wwm(i,j,k)
ppp(i,j,k) = c01k * c2r_star * ppm(i,j,k)
$ + dtm*(pt5*gc02(i,j,k)*(wwm(i,j,kp1)+wwm(i,j,k))
$ - gg0r(i,j,k)*(wwm(i,j,kp1)-wwm(i,j,k)))
$ - dt0 * div(i,j,k)
end do
end do
end do
*
*----------------------------------------------------------------------
return
end