Newer
Older
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! 8888888b. .d88888b. 888 888 d8888 8888888b. |
! 888 "Y88b d88P" "Y88b 888 888 d88888 888 Y88b |
! 888 888 888 888 888 888 d88P888 888 888 |
! 888 888 888 888 888 888 d88P 888 888 d88P |
! 888 888 888 888 888 888 d88P 888 8888888P" |
! 888 888 888 888 888 888 d88P 888 888 T88b |
! 888 .d88P Y88b. .d88P Y88b. .d88P d8888888888 888 T88b |
! 8888888P" "Y88888P" "Y88888P" d88P 888 888 T88b |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! DEFINE_BC_NEST May 2011 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine define_bc_nest (nnode,kfix,kfixt,x,y,z,u,v,w,temp,vo,params,nest)
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! this routine assigns the boundary condition for the boundaries of the nest in
! nested DOUAR simulations
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
use definitions
implicit none
integer :: nnode
integer :: kfix(nnode*3)
integer :: kfixt(nnode)
double precision :: x(nnode),y(nnode),z(nnode)
double precision :: u(nnode),v(nnode),w(nnode)
double precision :: temp(nnode)
type (void) vo
type (nest_info) nest
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
integer :: vermajor,verminor,verstat,verrev
Dave Whipp
committed
logical :: isostasyout,isobcout,nestout,compactionout
double precision :: eps,xls,yls,zls,current_time
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
eps=1.d-10
ovls%noctree=params%noctreemax
!if (params%debug.gt.1) call heap (threadinfo,'ov%octree','main',size(ov%octree),'int',+1)
! Read in coarse model octree; needs to be stored in params @ end of step
open (9,file=trim(nest%lsoutfile),status='old',form='unformatted')
! Read output version number
read (9) vermajor,verminor,verstat,verrev
! Read optional output flags
read (9) isostasyout,isobcout,nestout,compactionout
Dave Whipp
committed
read (9) ovls%octree(3),ovls%nnode,ovls%nleaves,nface,nlsf,np,current_time
allocate (ovls%x(ovls%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ovls%x in define_bc_nest$')
allocate (ovls%y(ovls%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ovls%y in define_bc_nest$')
allocate (ovls%z(ovls%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ovls%z in define_bc_nest$')
allocate (ovls%icon(8,ovls%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ovls%icon in define_bc_nest$')
allocate (ovls%unode(ovls%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ovls%unode in define_bc_nest$')
allocate (ovls%vnode(ovls%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ovls%vnode in define_bc_nest$')
allocate (ovls%wnode(ovls%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ovls%wnode in define_bc_nest$')
allocate (ovls%temp(ovls%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ovls%temp in define_bc_nest$')
allocate (ovls%pressure(ovls%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ovls%pressure in define_bc_nest$')
allocate (ovls%spressure(ovls%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ovls%spressure in define_bc_nest$')
allocate (ovls%temporary_nodal_pressure(ovls%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ovls%temp_nodal_pressure in define_bc_nest$')
allocate (ovls%whole_leaf_in_fluid(ovls%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ovls%whole_leaf_in_fluid in define_bc_nest$')
Dave Whipp
committed
read (9) ovls%x
read (9) ovls%y
read (9) ovls%z
read (9) ovls%unode
read (9) ovls%vnode
read (9) ovls%wnode
read (9)
read (9) ovls%temp
read (9) ovls%temporary_nodal_pressure
read (9)
read (9)
read (9)
if (isostasyout) then
!if (params%isostasy) then
! allocate (ovls%wnodeiso(ovls%nnode),stat=err)
! if (err.ne.0) call stop_run ('Error alloc ovls%wnodeiso in define_bc_nest$')
! read (9) (ovls%wnodeiso(i),i=1,ovls%nnode)
!else
Dave Whipp
committed
read (9)
!endif
endif
if (compactionout) then
!if (params%compaction) then
! allocate (ovls%wnodecompact(ovls%nnode),stat=err)
! if (err.ne.0) call stop_run ('Error alloc ovls%wnodecompact in define_bc_nest$')
! read (9) (ovls%wnodecompact(i),i=1,ovls%nnode)
!else
Dave Whipp
committed
read (9)
!endif
endif
Dave Whipp
committed
read (9) ovls%icon
read (9) ovls%pressure
read (9) ovls%spressure
read (9)
read (9)
read (9)
read (9)
read (9)
read (9)
read (9) ovls%whole_leaf_in_fluid
if (compactionout) then
!if (params%compaction) then
! allocate (ovls%compaction_density(ovls%nleaves),stat=err)
! if (err.ne.0) call stop_run ('Error alloc ovls%compaction_density in define_bc_nest$')
! read (9) (ovls%compaction_density(i),i=1,ovls%nleaves)
!else
Dave Whipp
committed
read (9)
!endif
endif
Dave Whipp
committed
read (9) ovls%octree
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
close (9)
! Loop through all nodes in nest and interpolate B/C velocities from LS model
do i=1,nnode
xls=x(i)*nest%sselemx+nest%xminls
yls=y(i)*nest%sselemy+nest%yminls
zls=z(i)*nest%sselemz+nest%zminls
if (x(i).lt.eps) then
kfix((i-1)*3+1)=1
kfix((i-1)*3+2)=1
kfix((i-1)*3+3)=1
call octree_interpolate_three (3,ovls%octree,ovls%noctree,ovls%icon, &
ovls%nleaves,ovls%nnode,xls,yls,zls, &
ovls%unode,u(i),ovls%vnode,v(i),ovls%wnode, &
w(i))
elseif (x(i).gt.1.d0-eps) then
kfix((i-1)*3+1)=1
kfix((i-1)*3+2)=1
kfix((i-1)*3+3)=1
call octree_interpolate_three (3,ovls%octree,ovls%noctree,ovls%icon, &
ovls%nleaves,ovls%nnode,xls,yls,zls, &
ovls%unode,u(i),ovls%vnode,v(i),ovls%wnode, &
w(i))
elseif (y(i).lt.eps) then
kfix((i-1)*3+1)=1
kfix((i-1)*3+2)=1
kfix((i-1)*3+3)=1
call octree_interpolate_three (3,ovls%octree,ovls%noctree,ovls%icon, &
ovls%nleaves,ovls%nnode,xls,yls,zls, &
ovls%unode,u(i),ovls%vnode,v(i),ovls%wnode, &
w(i))
elseif (y(i).gt.1.d0-eps) then
kfix((i-1)*3+1)=1
kfix((i-1)*3+2)=1
kfix((i-1)*3+3)=1
call octree_interpolate_three (3,ovls%octree,ovls%noctree,ovls%icon, &
ovls%nleaves,ovls%nnode,xls,yls,zls, &
ovls%unode,u(i),ovls%vnode,v(i),ovls%wnode, &
w(i))
elseif (z(i).lt.eps) then
kfix((i-1)*3+1)=1
kfix((i-1)*3+2)=1
kfix((i-1)*3+3)=1
call octree_interpolate_three (3,ovls%octree,ovls%noctree,ovls%icon, &
ovls%nleaves,ovls%nnode,xls,yls,zls, &
ovls%unode,u(i),ovls%vnode,v(i),ovls%wnode, &
w(i))
kfixt(i)=1
temp(i)=1.d0
elseif (z(i).gt.1.d0-eps) then
kfix((i-1)*3+1)=1
kfix((i-1)*3+2)=1
kfix((i-1)*3+3)=1
call octree_interpolate_three (3,ovls%octree,ovls%noctree,ovls%icon, &
ovls%nleaves,ovls%nnode,xls,yls,zls, &
ovls%unode,u(i),ovls%vnode,v(i),ovls%wnode, &
w(i))
kfixt(i)=1
temp(i)=0.d0
endif
if (.not.vo%influid(i)) then
kfixt(i)=1
temp(i)=0.d0
endif
end do
deallocate (ovls%x)
deallocate (ovls%y)
deallocate (ovls%z)
deallocate (ovls%icon)
deallocate (ovls%unode)
deallocate (ovls%vnode)
deallocate (ovls%wnode)
deallocate (ovls%temp)
deallocate (ovls%pressure)
deallocate (ovls%spressure)
deallocate (ovls%temporary_nodal_pressure)
deallocate (ovls%whole_leaf_in_fluid)
!if (isostasyout .and. params%isostasy) deallocate (ovls%wnodeiso)
!if (compactionout .and. params%compaction) deallocate(ovls%wnodecompact,ovls%compaction_density)
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------