From 837e8a584b36f061867cb91577df97605edda2ec Mon Sep 17 00:00:00 2001
From: Dave Whipp <dwhipp@dal.ca>
Date: Tue, 6 Sep 2011 16:15:03 +0000
Subject: [PATCH] Added modified segmented s-line geometry for use with DOUAR
v0.2+
---
src/define_bc_segmented_s-line.f90 | 212 +++++++++++++++++++++++++++++
1 file changed, 212 insertions(+)
create mode 100644 src/define_bc_segmented_s-line.f90
diff --git a/src/define_bc_segmented_s-line.f90 b/src/define_bc_segmented_s-line.f90
new file mode 100644
index 00000000..3eb364f8
--- /dev/null
+++ b/src/define_bc_segmented_s-line.f90
@@ -0,0 +1,212 @@
+!------------------------------------------------------------------------------|
+!------------------------------------------------------------------------------|
+! |
+! ||===\\ |
+! || \\ |
+! || || //==\\ || || //==|| ||/==\\ |
+! || || || || || || || || || || |
+! || // || || || || || || || |
+! ||===// \\==// \\==\\ \\==\\ || |
+! |
+!------------------------------------------------------------------------------|
+!------------------------------------------------------------------------------|
+! |
+! DEFINE_BC_SEGMENTED_S-LINE Feb. 2009 |
+! |
+!------------------------------------------------------------------------------|
+!------------------------------------------------------------------------------|
+
+subroutine define_bc_segmented_s-line(params,osolve,vo,bcdef,nest)
+
+!------------------------------------------------------------------------------|
+!(((((((((((((((( Purpose of the routine )))))))))))))))))))))))))))))))))))))
+!------------------------------------------------------------------------------|
+! This routine assigns the velocity boundary conditions for the segmented s-line
+! geometry
+
+!------------------------------------------------------------------------------|
+!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
+!------------------------------------------------------------------------------|
+
+use definitions
+!use mpi
+
+implicit none
+
+include 'mpif.h'
+
+type (parameters) params
+type (octreesolve) osolve
+type (void) vo
+type (bc_definition) bcdef
+type (nest_info) nest
+
+!------------------------------------------------------------------------------|
+!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
+!------------------------------------------------------------------------------|
+
+integer i,iproc,nproc,ierr
+double precision eps,lsf0,pi,lorig,h,x1,x2,phi,yend,cper,cscl,xstart,ystart,
+double precision theta,l,vin,vzfluxscl,cntvel,dxy,xend,xsym,ymax,xwidth,ywidth
+double precision,dimension(:),allocatable :: x0,ldisp
+integer ie,ij,j,jp,nelemx,nelemz
+logical firstcall
+
+!------------------------------------------------------------------------------|
+!------------------------------------------------------------------------------|
+
+call mpi_comm_size (mpi_comm_world,nproc,ierr)
+call mpi_comm_rank (mpi_comm_world,iproc,ierr)
+
+eps=1.d-10
+osolve%kfix=0
+osolve%kfixt=0
+pi=atan(1.d0)*4.d0
+
+l=bcdef%bc_parameter(1)
+phi=bcdef%bc_parameter(2)
+theta=bcdef%bc_parameter(3)
+xsym=bcdef%bc_parameter(4)
+ystart=bcdef%bc_parameter(5)
+yend=bcdef%bc_parameter(6)
+ymax=bcdef%bc_parameter(7)
+vin=bcdef%bc_parameter(8)
+nelemx=int(bcdef%bc_parameter(9))
+nelemz=int(bcdef%bc_parameter(10))
+
+ywidth=yend-ystart
+phi=phi/180.d0*pi
+theta=theta/180.d0*pi
+xwidth=ywidth*tan(theta)
+xstart=xsym-(xwidth/2.d0)
+xend=xsym+(xwidth/2.d0)
+vzfluxscl=0.5d0*(1.d0+(tan(phi)/tan((pi-phi)/2.d0)))
+nb=2**params%levelmax_oct
+dxy=1.d0/2**(params%levelmax_oct+1.d0)
+
+select case(trim(params%infile))
+case('input.txt','input.small.txt')
+ allocate(x0(osolve%nnode))
+ do i=1,osolve%nnode
+ if (osolve%y(i).le.ystart) then
+ x0(i)=xstart
+ else if (osolve%y(i).le.yend) then
+ x0(i)=(osolve%y(i)-ystart)*tan(theta)+xstart
+ else
+ x0(i)=xend
+ endif
+ enddo
+
+ if (params%isobc) then
+ allocate (ldisp(osolve%nnode))
+ do i=1,osolve%nnode
+ xdisp=nint(x0(i)*nb)+1
+ ydisp=nint(osolve%y(i)*nb)+1
+ ldisp(i)=l+bcdef%zisodisp(xdisp,ydisp)
+ enddo
+
+ do i=1,osolve%nnode
+ cntvel=0.d0
+ if (params%isobc) l=ldisp(i)
+ x1=x0(i)-l/tan(phi)
+ x2=x0(i)+l/tan(phi)
+ if (((x2-x0(i))-(x0(i)-x1)).gt.eps) vzfluxscl=1.d0 ! Set vzfluxscl equal to one if B/Cs are not symmetric
+
+ if (osolve%x(i).lt.eps) then
+ osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=0.d0
+ osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
+ osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
+ if (osolve%y(i).le.(ymax-(real(nelemx)*dxy))) then
+ osolve%u(i)=1.d0
+ elseif (osolve%y(i).le.(ymax+(real(nelemx)*dxy))) then
+ osolve%u(i)=1.d0-(osolve%y(i)-(ymax-(real(nelemx)*dxy)))/(2.d0*real(nelemx)*dxy)
+ endif
+ endif
+ if (osolve%x(i).gt.1.d0-eps) then
+ osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=0.d0
+ osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
+ osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
+ endif
+ if (osolve%y(i).lt.eps) then
+ !osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=0.d0
+ osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
+ !osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
+ endif
+ if (osolve%y(i).gt.1.d0-eps) then
+ !osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=0.d0
+ osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
+ !osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
+ endif
+ if (osolve%z(i).lt.eps) then
+ osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=0.d0
+ osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
+ osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
+
+ if (osolve%x(i).le.(x1-real(nelemz)*dxy)) then
+ if (osolve%y(i).le.(ymax-(real(nelemx)*dxy))) then
+ osolve%u(i)=1.d0
+ elseif (osolve%y(i).le.(ymax+(real(nelemx)*dxy))) then
+ osolve%u(i)=1.d0-(osolve%y(i)-(ymax-(real(nelemx)*dxy)))/(2.d0*real(nelemx)*dxy)
+ endif
+
+ elseif (osolve%x(i).le.(x1+real(nelemz)*dxy)) then
+ if (osolve%y(i).le.(ymax-(real(nelemx)*dxy))) then
+ osolve%u(i)=1.d0
+ elseif (osolve%y(i).le.(ymax+(real(nelemx)*dxy))) then
+ osolve%u(i)=1.d0-(osolve%y(i)-(ymax-(real(nelemx)*dxy)))/(2.d0*real(nelemx)*dxy)
+ endif
+ osolve%w(i)=osolve%u(i)*(osolve%x(i)-(x1-(real(nelemz)*dxy)))/(2.d0*real(nelemz)*dxy)*(vzfluxscl*-sin(phi))
+ osolve%u(i)=osolve%u(i)*(((1.d0-(osolve%x(i)-(x1-(real(nelemz)*dxy)))/(2.d0*real(nelemz)*dxy))*(1.d0-vzfluxscl*cos(phi)))+vzfluxscl*cos(phi))
+
+ elseif (osolve%x(i).le.(x2-real(nelemz)*dxy)) then
+ if (osolve%y(i).le.(ymax-(real(nelemx)*dxy))) then
+ osolve%u(i)=1.d0
+ elseif (osolve%y(i).le.(ymax+(real(nelemx)*dxy))) then
+ osolve%u(i)=1.d0-(osolve%y(i)-(ymax-(real(nelemx)*dxy)))/(2.d0*real(nelemx)*dxy)
+ endif
+ osolve%w(i)=osolve%u(i)*vzfluxscl*-sin(phi)
+ osolve%u(i)=osolve%u(i)*vzfluxscl*cos(phi)
+
+ elseif (osolve%x(i).le.(x2+real(nelemz)*dxy)) then
+ if (osolve%y(i).le.(ymax-(real(nelemx)*dxy))) then
+ osolve%u(i)=1.d0
+ elseif (osolve%y(i).le.(ymax+(real(nelemx)*dxy))) then
+ osolve%u(i)=1.d0-(osolve%y(i)-(ymax-(real(nelemx)*dxy)))/(2.d0*real(nelemx)*dxy)
+ endif
+ osolve%w(i)=osolve%u(i)*((1.d0-(osolve%x(i)-(x2-(real(nelemz)*dxy)))/(2.d0*real(nelemz)*dxy))*(vzfluxscl*-sin(phi)))
+ osolve%u(i)=osolve%u(i)*((1.d0-(osolve%x(i)-(x2-(real(nelemz)*dxy)))/(2.d0*real(nelemz)*dxy))*(vzfluxscl*cos(phi)))
+ endif
+
+ osolve%kfixt(i)=1
+ osolve%temp(i)=1.d0
+ endif
+ if (osolve%z(i).gt.1.d0-eps) then
+ osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=0.d0
+ osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
+ osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
+ osolve%kfixt(i)=1
+ osolve%temp(i)=0.d0
+ endif
+ if (.not.vo%influid(i)) then
+ osolve%kfixt(i)=1
+ osolve%temp(i)=0.d0
+ endif
+ if (abs(bcdef%utrans).gt.eps) then
+ osolve%u(i)=osolve%u(i)+bcdef%utrans
+ endif
+ if (abs(bcdef%vtrans.gt.eps) then
+ osolve%v(i)=osolve%v(i)+bcdef%vtrans
+ endif
+ osolve%u(i)=osolve%u(i)*vin
+ osolve%v(i)=osolve%v(i)*vin
+ osolve%w(i)=osolve%w(i)*vin
+ enddo
+
+ if (params%isobc) then
+ call define_isostasy_bc(params,osolve,bcdef)
+ deallocate(ldisp)
+ endif
+ deallocate(x0)
+
+end
+!------------------------------------------------------------------------------|
\ No newline at end of file
--
GitLab