copyright (C) 2001 MSC-RPN COMM %%%MC2%%%
***s/r rhs3_r0
*
subroutine rhs3_r0 ( w1 ) 1,1
implicit none
*
real*8 w1(*)
*
*AUTHORs C. Girard & M. Desgagne
*
*OBJECT
*
* *******************************************************************
* * *
* * CALCULATES *
* * right-hand-side CORIOLIS+KE+NON-LINEAR *
* * *
* * explicit terms: *
* * Ru, Rv, Rw, Rb, Rp *
* * FROM *
* * centered time-level t variables: *
* * u0, v0, w0, b0, p0, hu0, tr0 *
* * STORES *
* * results in *
* * ur, vr, wr, br, pr *
* * *
* * PERFORMS essentially the following calculations: *
* * *
* * *
* * bv0 = bv0(b0,hu0,tr0) *
* * *
* * __XY ____Y _XY ___XZ ~ *
* * ur = + v0 *fcor - KE*d(S )/dX - c04*bv0 * d(p0)/dX *
* * *
* * *
* * __XY ____X _XY ___YZ ~ *
* * vr = - u0 *fcor - KE*d(S )/dY - c04*bv0 * d(p0)/dY *
* * *
* * ~ *
* * wr = b0*alfa/(1+alfa) + c02*(bv0-b0) - c04*bv0 * d(p0)/dZ *
* * *
* * _____________Z *
* * br = - c06*n02g * w0 * alfa - c07 * b0 * DIV(u0,v0,w0) *
* * *
* * __Z ____Z *
* * pr = gc02 * w0 * alfa *
* * *
* * *
* * ~ ~ ~ *
* * SEE - rhs3_pm for description of d/dX, d/dY, d/dZ *
* * *
* * - divcal for calculation of DIV *
* * *
* * - cxxpar for definitions of parameters c02,c04,c06,c07 *
* * *
* * - nacbar for calculation of alfa, n02g, gc02 *
* * *
* * *
* *******************************************************************
*
#include "consdyn_8.cdk"
#include "yomdyn.cdk"
#include "yomdyn1.cdk"
#include "dynpar.cdk"
#include "dynmem.cdk"
#include "topo.cdk"
#include "wrnmem.cdk"
#include "physcom.cdk"
#include "nbcpu.cdk"
#include "alfmem.cdk"
*
integer i, j, k, n, kkk,kp1, km1, id, jd, iff, jf
real*8 div,bv0,g1d,g2d,t1
pointer (padiv, div(minx:maxx,miny:maxy,*)),
$ (pabv0, bv0(minx:maxx,miny:maxy,*)),
$ (pag1d, g1d(minx:maxx,miny:maxy,*)),
$ (pag2d, g2d(minx:maxx,miny:maxy,*)),
$ (pat1, t1 (minx:maxx,miny:maxy,*))
real*8 p25, pt5, one, two, coe1,
$ d1, vbarxy, ubarxy, rtbarxz, rtbaryz, kinetic
parameter (p25=0.25d0, pt5=0.5d0, one=1.0d0, two=2.0d0)
*
*----------------------------------------------------------------------
*
padiv = loc(w1( 1))
pabv0 = loc(w1( dim3d+1))
pag1d = loc(w1(2*dim3d+1))
pag2d = loc(w1(3*dim3d+1))
pat1 = loc(w1(4*dim3d+1))
*
!$omp do
do k = 1, gnk
*
id = (1-hx)*west
jd = (1-hy)*south
iff = ldni+(hx-1)*east
jf = ldnj+(hy-1)*north
*
do j = jd,jf
do i = id,iff
t1(i,j,k) = 0.0
end do
end do
*
if (k.gt.1.and.gnload.eq.1) then
do n = bghyd, edhyd
do j = jd,jf
do i = id,iff
t1(i,j,k) =t1(i,j,k) + tr0(i,j,k,n)
end do
end do
end do
endif
do j = jd,jf
do i = id,iff
d1 = bb0(i,j,k) + grav_8
bv0(i,j,k) = bb0(i,j,k) + d1 *
$ ( hu0(i,j,k)*(rgasv_8/rgasd_8-one)-t1(i,j,k))
end do
end do
*
id = 1-hx*west
jd = 1-hy*south
iff = ldni+(hx-1)*east
jf = ldnj+(hy-1)*north
*
coe1 = 0.5
if(k.eq.1.or.k.eq.gnk) coe1 = 1.
km1 = max( k-1, 1 )
*
do j = jd, jf
do i = id+west, iff
g1d(i,j,k) = -gg1(i,j,k) * coe1 *
$ ( pp0(i,j,k ) + pp0(i-1,j,k )
$ - pp0(i,j,k-1) - pp0(i-1,j,k-1) )
end do
end do
do j = jd+south, jf
do i = id, iff
g2d(i,j,k) = -gg2(i,j,k) * coe1 *
$ ( pp0(i,j,k ) + pp0(i,j-1,k )
$ - pp0(i,j,k-1) - pp0(i,j-1,k-1) )
end do
end do
*
end do
!$omp enddo
*
call divcal2 ( div,uu0,vv0,ww0,sbxy,gg1,gg2,gg0r,g0ur,g0vr,t1,
$ odx,ody,minx,maxx,miny,maxy,gnk,id,jd,iff,jf )
*
!$omp do
do k = 1, gnk
km1 = max( k-1, 1 )
*
if (k.lt.gnk) then
*
do j = jd, jf
do i = id+west, iff
vbarxy = p25 * (vv0(i-1,j+1,k)+vv0(i-1,j,k ) +
$ vv0(i ,j+1,k)+vv0(i ,j,k ) )
rtbarxz = p25 * (bv0(i ,j,k+1)+bv0(i-1,j,k+1) +
$ bv0(i ,j,k )+bv0(i-1,j,k ) ) * c04
kinetic = ( uu0(i,j,k)*uu0(i,j,k) + vbarxy*vbarxy ) * pt5
uur(i,j,k) =
$ + vbarxy * ( fcor(i,j+1) + fcor(i,j) ) * pt5
$ - kinetic * ( sbxy(i,j) - sbxy(i-1,j) ) * odxu(1)
$ - rtbarxz * ( (pp0(i,j,k)-pp0(i-1,j,k)) * odxu(1)
$ + g0ur(i,j,k) * ( g1d(i,j,k+1)+g1d(i,j,k)) * pt5 )
end do
end do
*
do j = jd+south, jf
do i = id, iff
ubarxy = p25 * (uu0(i+1,j-1,k)+uu0(i,j-1,k ) +
$ uu0(i+1,j ,k)+uu0(i,j ,k ) )
rtbaryz = p25 * (bv0(i ,j,k+1)+bv0(i,j-1,k+1) +
$ bv0(i ,j,k )+bv0(i,j-1,k ) ) * c04
kinetic = ( vv0(i,j,k)*vv0(i,j,k) + ubarxy*ubarxy ) * pt5
vvr(i,j,k) =
$ - ubarxy * ( fcor(i+1,j) + fcor(i,j) ) * pt5
$ - kinetic * ( sbxy(i,j) - sbxy(i,j-1) ) * odyv(j)
$ - rtbaryz * ( (pp0(i,j,k)-pp0(i,j-1,k)) * odyv(j)
$ + g0vr(i,j,k) * ( g2d(i,j,k+1)+g2d(i,j,k)) * pt5 )
end do
end do
*
endif
*
kp1 = min( k+1, gnk )
kkk = min( k, gnk-1 )
*
do j = jd, jf
do i = id,iff
wwr(i,j,k) = - c04 * bv0(i,j,k) *
$ (pp0(i,j,k)-pp0(i,j,k-1))*g0wr(i,j,k)
$ + bb0(i,j,k) * alfa(i,j,k) * ap1r(i,j,k)
$ + c02 * ( bv0(i,j,k) - bb0(i,j,k ) )
bbr(i,j,k) = - c07*bb0(i,j,k)*(div(i,j,kkk)+div(i,j,km1))*pt5
$ - c06 * n02g(i,j,k) * ww0(i,j,k) * alfa(i,j,k)
ppr(i,j,k) = p25*gc02(i,j,k)*( ww0(i,j,k)+ ww0(i,j,kp1))
$ *(alfa(i,j,k)+alfa(i,j,kp1))
end do
end do
*
if (k.eq.1) then
do j = jd, jf
do i = id+west, iff
uur(i,j,1) = uur(i,j,1) - pt5
$ * g0ur(i,j,1 ) * g1d(i,j,1 )
end do
end do
do j = jd+south, jf
do i = id, iff
vvr(i,j,1) = vvr(i,j,1) - pt5
$ * g0vr(i,j,1 ) * g2d(i,j,1 )
end do
end do
endif
if (k.eq.gnk-1) then
do j = jd, jf
do i = id+west, iff
uur(i,j,gnk-1) = uur(i,j,gnk-1) - pt5
$ * g0ur(i,j,gnk) * g1d(i,j,gnk)
end do
end do
do j = jd+south, jf
do i = id, iff
vvr(i,j,gnk-1) = vvr(i,j,gnk-1) - pt5
$ * g0vr(i,j,gnk) * g2d(i,j,gnk)
end do
end do
endif
*
end do
!$omp enddo
*----------------------------------------------------------------------
*
return
end