copyright (C) 2001 MSC-RPN COMM %%%MC2%%%
*
subroutine rhs3_pm ( dimphy, wrk, lminx, lmaxx, lminy, lmaxy ) 1,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-level t-Dt variables: *
* * um, vm, wm, bm, pm *
* * STORES *
* * results in *
* * up=Pu, vp=Pv, wp=Pw, bp=Pb, pp=Pp *
* * *
* * PERFORMS essentially the following calculations: *
* * *
* * ~ *
* * up = um - dtm * d(pm)/dX *
* * *
* * ~ *
* * vp = vm - dtm * d(pm)/dY *
* * *
* * ~ *
* * wp = c03 * wm - dtm * [ d(pm)/dZ - bm /(1+alfa) ] *
* * *
* * __Z *
* * bp = bm - gama_star * pm - dtm * c06*n02g * wm *
* * *
* * __Z *
* * pp = pm*c01k/c2_star + dtm*gc02* wm - dtm * DIV(um,vm,wm) *
* * *
* * *
* * ******************************************* *
* * * ________Z * *
* * * ~ _X * *
* * * dq/dX = dq/dX + G0ur*G1*dq/dZ * *
* * * * *
* * * ________Z * *
* * * ~ _Y * *
* * * dq/dY = dq/dY + G0vr*G2*dq/dZ * *
* * * * *
* * * ~ * *
* * * dq/dz = G0wr*dq/dZ * *
* * * * *
* * ******************************************* *
* * *
* * *
* * SEE - divcal for calculation of DIV *
* * *
* * - cxxpar for definitions of parameters c01, c03, c06 *
* * and gama_star, c2_star *
* * *
* * -1 *
* * - nacbar for calculation of ap1r=(1+alfa), n02g, gc02 *
* * *
* * *
* *******************************************************************
*
#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, fulltp, 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))
*
id =1-hx*west
jd =1-hy*south
iff=ldni+(hx-1)*east
jf =ldnj+(hy-1)*north
*
call divcal2 ( div,uum,vvm,wwm,sbxy,gg1,gg2,gg0r,g0ur,g0vr,q3,
$ odx,ody,minx,maxx,miny,maxy,gnk,id,jd,iff,jf )
!$omp do
do k = gnk, 1, -1
km1 = max( k-1 , 1 )
if (k.lt.gnk) then
if(topo_folwing) then
kp1 = min( k+1 , gnk-1 )
do j = jd, jf
do i = id+west, iff
uup(i,j,k) = uum(i,j,k) - dtm * (
$ ( ppm(i,j,k ) - ppm(i-1,j,k) ) * odxu(1) +
$ g0ur(i,j,k) * p25 * (-gg1(i,j,k+1) *
$ ( ppm(i,j,kp1) + ppm(i-1,j,kp1)
$ - ppm(i,j,k ) - ppm(i-1,j,k ) ) -
$ gg1(i,j,k ) *
$ ( ppm(i,j,k ) + ppm(i-1,j,k )
$ - ppm(i,j,km1) - ppm(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) - dtm * (
$ ( ppm(i,j,k ) - ppm(i,j-1,k) ) * odyv(j) +
$ g0vr(i,j,k) * p25 * (-gg2(i,j,k+1) *
$ ( ppm(i,j,kp1) + ppm(i,j-1,kp1)
$ - ppm(i,j,k ) - ppm(i,j-1,k ) ) -
$ gg2(i,j,k ) *
$ ( ppm(i,j,k ) + ppm(i,j-1,k )
$ - ppm(i,j,km1) - ppm(i,j-1,km1) )) )
end do
end do
else
do j = jd, jf
do i = id+west, iff
uup(i,j,k) = uum(i,j,k) - dtm *
$ ( ppm(i,j,k ) - ppm(i-1,j,k) ) * odxu(1)
end do
end do
do j = jd+south, jf
do i = id, iff
vvp(i,j,k) = vvm(i,j,k) - dtm *
$ ( ppm(i,j,k ) - ppm(i,j-1,k) ) * odyv(j)
end do
end do
endif
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*gc02(i,j,k)*pt5*(wwm(i,j,k)+wwm(i,j,kp1))
$ - dtm*div(i,j,k)
end do
end do
end do
!$omp enddo
*
if ((dimphy.gt.0).and.(gnpfb.gt.1)) then
!$omp do
do k = 1, gnk-1
do j = jd, jf
do i = id+west, iff
uup (i,j,k) = uup (i,j,k) + utdp2 (i,j,k)
end do
end do
do j = jd+south, jf
do i = id, iff
vvp (i,j,k) = vvp (i,j,k) + vtdp2 (i,j,k)
end do
end do
do j = jd, jf
do i = id, iff
fulltp = bbm(i,j,k)+grav_8
ppp(i,j,k) = ppp(i,j,k) + ttdp2 (i,j,k)/fulltp
bbp(i,j,k) = bbp(i,j,k) + ttdp2 (i,j,k)
$ * (1 + c07*bbm(i,j,k)/fulltp)
hup(i,j,k) = hup(i,j,k) + hutdp2(i,j,k)
end do
end do
if (diffuw) then
do j = jd, jf
do i = id, iff
wwp (i,j,k) = wwp (i,j,k) + swtdp2 (i,j,k)
end do
end do
endif
do n = bghyd, edhyd
do j = jd, jf
do i = id, iff
trp(i,j,k,n) = trp(i,j,k,n) + cltdp2(i,j,k,n)
end do
end do
end do
end do
!$omp enddo
*
pas = pautp1
pad = pautp2
!$omp do
do i=1,dimphy/2
d(i) = s(i)
end do
!$omp enddo
endif
*
*----------------------------------------------------------------------
return
end