Renamed

src/define_bcs/define_bc_segmented_s-line.f90

-!------------------------------------------------------------------------------|
-!------------------------------------------------------------------------------|
-!                                                                              |
-!              ||===\\                                                         | 
-!              ||    \\                                                        |
-!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
-!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
-!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
-!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
-!                                                                              |
-!------------------------------------------------------------------------------|
-!------------------------------------------------------------------------------|
-!                                                                              |
-!              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)
-
-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
-!------------------------------------------------------------------------------|