copyright (C) 2001 MSC-RPN COMM %%%MC2%%%
***s/r bulle_3d
*
subroutine bulle_3d 1,10
implicit none
*
**
#include "nesting.cdk"
#include "dynmem.cdk"
#include "tracers.cdk"
#include "nestpnt.cdk"
#include "yomdyn.cdk"
#include "yomdyn1.cdk"
#include "consdyn_8.cdk"
#include "levels.cdk"
#include "partopo.cdk"
#include "bubble.cdk"
#include "topo.cdk"
#include "h_geop.cdk"
#include "bcsdim.cdk"
#include "bcsmem.cdk"
*
integer i,j,k,n,gnix,gnjy,dim,err,
$ id,indn,indw,inde,ofn,ofw,ofe
real hhm(minx:maxx,miny:maxy,0:gnk),
$ hht(minx:maxx,miny:maxy, gnk)
real tp1,th1,s,d
pointer (patp1, tp1(gni+2*hx,gnj+2*hy,gnk)),
$ (path1, th1(gni+2*hx,gnj+2*hy,gnk)),
$ (pas ,s(*)), (pad ,d(*))
real*8 p00,t00,dtt,hm1,r,rx,rz,one,c00,pistr,hhh,con1
parameter(one=1.d0)
*
*---------------------------------------------------------------------
*
gnix = gni+2*hx
gnjy = gnj+2*hy
if (myproc.eq.0) then
call hpalloc (patp1, gnix*gnjy*gnk , err, 1)
call hpalloc (path1, gnix*gnjy*gnk , err, 1)
endif
*
p00=100000.d0
t00=303.16
c00=grav_8/(cpd_8*t00)
*
if (myproc.eq.0) then
do k=1,gnk
do j=1,gnjy
do i=1,gnix
dtt = 0.
r=sqrt ((float(i)-bb_xcntr)**2.+(float(k)-bb_zcntr)**2.)
if(r.le.bb_radius) dtt = bb_dpth
ccccc square bubble cccccccccccccccccccccccc
c rx=abs(float(i)-bb_xcntr)
c rz=abs(float(k)-bb_zcntr)
c if(rx.le.bb_radius.and.rz.le.bb_radius) dtt = bb_dpth
ccccc square bubble cccccccccccccccccccccccc
tp1 (i,j,k)= dtt
th1 (i,j,k)= dtt
end do
end do
end do
endif
*
call glbdist2 (tp1,bbp ,minx,maxx,miny,maxy,gni+hx,gnj+hy,gnk)
call glbdist2 (th1,trp(1-hx,1-hy,1,1) ,minx,maxx,miny,maxy,
$ gni+hx,gnj+hy,gnk)
*
do j=1-hy,ldnj+hy
do i=1-hx,ldni+hx
hhm(i,j,0)=hh0i(i,j,1)
ppp(i,j,0)=dlog(p00)
end do
end do
*
do k=1,gnk
do j=1-hy,ldnj+hy
do i=1-hx,ldni+hx
hht(i,j,k)=h_geop(ztr(k),i,j)
hhh = hht(i,j,k)
pistr=one-c00*hhh
bbp(i,j,k)=t00*pistr+bbp(i,j,k)
hhm(i,j,k)=h_geop(zm(k),i,j)
hhh = hhm(i,j,k)
pistr=one-c00*hhh
ppp(i,j,k)=dlog(p00)+cpd_8/rgasd_8*dlog(pistr)
end do
end do
end do
*
call qntstar(qstr,nssq,gots,orts,hht,hhm,
$ (maxx-minx+1)*(maxy-miny+1),0,gnk)
*
do j=1-hy,ldnj+hy
do i=1-hx,ldni+hx
ppp(i,j,0) = (ppp(i,j,0) - qstr(1,1,0))*rgasd_8*grtstar
end do
end do
*
con1=grav_8/grtstar
do k=1,gnk
do j=1-hy,ldnj+hy
do i=1-hx,ldni+hx
ppp(i,j,k) = (ppp(i,j,k)-qstr(1,1,k))*rgasd_8*grtstar
bbp(i,j,k) = con1*(bbp(i,j,k)-grtstar)
end do
end do
end do
*
* Initialize moisture and velocities to zero
*
do k=1,gnk
do j=1-hy,ldnj+hy
do i=1-hx,ldni+hx
uup(i,j,k) = 0.
vvp(i,j,k) = 0.
hup(i,j,k) = 0.
wwp(i,j,k) = 0.
end do
end do
end do
*
dim = ndynvar*dim3d+dim2d
pas = pappp
pad = papp0
do i=1,dim
d(i) = s(i)
end do
*
ofn = (bcs_in-bcs_is)/gnk
ofw = (bcs_iw-bcs_in)/gnk
ofe = (bcs_ie-bcs_iw)/gnk
indn = bcs_in + ofn
indw = bcs_iw + ofn + ofw
inde = bcs_ie + ofn + ofw + ofe
call trnes (uup,bcs_uu(bcs_is),bcs_uu(bcs_in),bcs_uu(bcs_iw),
$ bcs_uu(bcs_ie),minx,maxx,miny,maxy,
$ minxs,maxxs,minys,maxys,minxw,maxxw,minyw,maxyw,gnk)
call trnes (vvp,bcs_vv(bcs_is),bcs_vv(bcs_in),bcs_vv(bcs_iw),
$ bcs_vv(bcs_ie),minx,maxx,miny,maxy,
$ minxs,maxxs,minys,maxys,minxw,maxxw,minyw,maxyw,gnk)
call trnes (wwp,bcs_ww(bcs_is),bcs_ww(bcs_in),bcs_ww(bcs_iw),
$ bcs_ww(bcs_ie),minx,maxx,miny,maxy,
$ minxs,maxxs,minys,maxys,minxw,maxxw,minyw,maxyw,gnk)
call trnes (bbp,bcs_bb(bcs_is),bcs_bb(bcs_in),bcs_bb(bcs_iw),
$ bcs_bb(bcs_ie),minx,maxx,miny,maxy,
$ minxs,maxxs,minys,maxys,minxw,maxxw,minyw,maxyw,gnk)
call trnes (ppp,bcs_pp(bcs_is),bcs_pp(indn) ,bcs_pp(indw) ,
$ bcs_pp(inde) ,minx,maxx,miny,maxy,
$ minxs,maxxs,minys,maxys,minxw,maxxw,minyw,maxyw,gnk+1)
call trnes (hup,bcs_hu(bcs_is),bcs_hu(bcs_in),bcs_hu(bcs_iw),
$ bcs_hu(bcs_ie),minx,maxx,miny,maxy,
$ minxs,maxxs,minys,maxys,minxw,maxxw,minyw,maxyw,gnk)
do n = 1, ntr
id = (n-1)*bcs_sz+1
call trnes (trp(1-hx,1-hy,1,n),bcs_tr(id),bcs_tr(id+bcs_in-1),
$ bcs_tr(id+bcs_iw-1),bcs_tr(id+bcs_ie-1),minx,maxx,miny,maxy,
$ minxs,maxxs,minys,maxys,minxw,maxxw,minyw,maxyw,gnk)
end do
if (ctebcs) then
do i=1,(ntr+6)*bcs_sz+bcs_sz/gnk
bcs_ppa(i) = bcs_pp(i)
end do
endif
*
if (gnpilver.gt.0) then
do k=1,gnk
do j=1-hy,ldnj+hy
do i=1-hx,ldni+hx
ppntt(i,j,k) = ppp(i,j,k)
uuntt(i,j,k) = uup(i,j,k)
vvntt(i,j,k) = vvp(i,j,k)
wwntt(i,j,k) = wwp(i,j,k)
bbntt(i,j,k) = bbp(i,j,k)
huntt(i,j,k) = hup(i,j,k)
end do
end do
end do
do n = 1, ntr
do k=1,gnk
do j=1-hy,ldnj+hy
do i=1-hx,ldni+hx
trntt(i,j,k,n) = trp(i,j,k,n)
end do
end do
end do
end do
do j=1-hy,ldnj+hy
do i=1-hx,ldni+hx
ppntt(i,j,0) = ppp(i,j,0)
end do
end do
*
if (ctebcs) then
do k=1,gnk
do j=1-hy,ldnj+hy
do i=1-hx,ldni+hx
ppnta(i,j,k) = ppp(i,j,k)
uunta(i,j,k) = uup(i,j,k)
vvnta(i,j,k) = vvp(i,j,k)
wwnta(i,j,k) = wwp(i,j,k)
bbnta(i,j,k) = bbp(i,j,k)
hunta(i,j,k) = hup(i,j,k)
end do
end do
do n=1,ntr
do j=1-hy,ldnj+hy
do i=1-hx,ldni+hx
trnta(i,j,k,n) = trp(i,j,k,n)
end do
end do
end do
end do
do j=1-hy,ldnj+hy
do i=1-hx,ldni+hx
ppnta(i,j,0) = ppp(i,j,0)
end do
end do
endif
endif
*
if (myproc.eq.0) then
write (6,602)
call hpdeallc (patp1 ,err,1)
call hpdeallc (path1 ,err,1)
endif
*
602 format (/' BUBBLE INITIALIZATION COMPLETED')
*---------------------------------------------------------------------
return
end