Skip to content
Snippets Groups Projects
Select Git revision
  • feature/mumps
  • develop
  • feature/parabolic_subduction
  • feature/hdf5out
  • feature/optimisedMoveCloud
  • feature/thrust_isobc
  • feature/particleAdvectionFast
  • feature/orography_sinclair
  • feature/mixed_cloudmat
  • feature/improved_cloud_material_numbering
  • feature/multibox2
  • feature/EGU2020
  • feature/mumps_final
  • feature/oroprecip_tuebingen
  • MfN_Berlin
  • feature/cloudmat-sphere
  • feature/baselevel
  • feature/LK_manuscript_2018
  • feature/sb04_oroprecip
  • feature/oroprecip2
  • v0.5.1
  • v0.5
  • v0.4.3.0
  • v0.4.2.0
  • v0.4.1.0
  • v0.4.0.2
  • v0.4.0.1
  • v0.4
  • v0.4rc1
  • DOUAR-WSMP_v0.1.1
  • DOUAR-WSMP_v0.1
  • 2009-09-02
  • fromDavidWhipp
33 results
Blame Permalink
define_bc.20090710.f90 9.51 KiB 
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              DEFINE_BC_TIBET   Feb. 2009                                     |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

subroutine define_bc (params,osolve,vo) 

!------------------------------------------------------------------------------|
!((((((((((((((((  Purpose of the routine  )))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! this routine assigns the boundary condition for the Tibet experiment

!------------------------------------------------------------------------------|
!((((((((((((((((  declaration of the subroutine arguments  ))))))))))))))))))))
!------------------------------------------------------------------------------|

use definitions

implicit none

type (parameters) params
type (octreesolve) osolve
type (void) vo

!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|

integer i,iproc,nproc,ierr
double precision eps,lsf0,pi,l,h,x0,x1,x2,phi,spu,spv,crustfact
integer ie,ij,j,jp

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|


INCLUDE 'mpif.h'

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=0.2945d0
h=0.0625d0
x0=0.5d0
phi=45.d0
phi=phi/180.d0*pi
crustfact=0.d0                                                                  
! Factor for what proportion of the crust will be subducted
spu=1.d0                                                                        
! x-velocity of s-point relative to incoming plate
spv=0.d0                                                                        
! y-velocity of s-point relative to incoming plate
l=l+crustfact*h                                                                 
! Apply crustal subduction factor
!x1=x0-l*(1.d0-cos(phi))/sin(phi)
x1=x0+l*(1.d0-cos(phi))/sin(phi)-(1.d0-spu)*params%dt
!x2=x0+l/tan(phi)
x2=x0-l/tan(phi)-(1.d0-spu)*params%dt
select case(trim(params%infile))

case('input.txt','input.small.txt')

  do i=1,osolve%nnode
     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.0.575d0) osolve%u(i)=(1.d0- sin(pi*(osolve%y(i)+.425d0)/2.d0)**100)
        !if (osolve%y(i).ge.0.425d0) osolve%u(i)=(1.d0- sin(pi*(osolve%y(i)+.575d0)/2.d0)**100)
     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)-x1)*(osolve%x(i)-x2).le.0.d0) then
          !if (osolve%y(i).le.0.575d0) osolve%u(i)=cos(phi)*(1.d0-sin(pi*(osolve%y(i)+.425d0)/2.d0)**100)
          !if (osolve%y(i).ge.0.425d0) osolve%u(i)=cos(phi)*(1.d0-sin(pi*(osolve%y(i)+.575d0)/2.d0)**100)
          if (osolve%y(i).le.0.575d0) osolve%w(i)=-sin(phi)*(1.d0-sin(pi*(osolve%y(i)+.425d0)/2.d0)**100)
          !if (osolve%y(i).ge.0.425d0) osolve%w(i)=-sin(phi)*(1.d0-sin(pi*(osolve%y(i)+.575d0)/2.d0)**100)
        elseif (osolve%x(i).lt.x2) then
          if (osolve%y(i).le.0.575d0) osolve%u(i)=(1.d0-sin(pi*(osolve%y(i)+.425d0)/2.d0)**100)
          !if (osolve%y(i).ge.0.425d0) osolve%u(i)=(1.d0-sin(pi*(osolve%y(i)+.575d0)/2.d0)**100)
        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 (spu.ne.0.d0) then                                                     ! 
!If x-comp of s-point velocity is non-zero then
         osolve%u(i)=osolve%u(i)-spu                                           ! 
!subract s-point velocity from entire model domain
     endif                                                                     ! 
!to keep that point stationary
     if (spv.ne.0.d0) then                                                     ! 
!Same condition as above in y-direction
         osolve%v(i)=osolve%v(i)-spv
     endif
  enddo

case ('input.jgr')
  call define_bc_jgr (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo) 

case ('input.sphere')
  call define_bc_sphere (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo) 

case ('input.3Dpunch')
  call define_bc_3Dpunch (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo) 

case ('input.2Dpunch')
  call define_bc_2Dpunch (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo) 

case ('input.2Dpunch_vert')
  call define_bc_2Dpunch_vert (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo) 

case ('input.folding')
  call define_bc_folding (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo) 

case ('input.pipo')
  call define_bc_pipo (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo) 

case ('input.ritske')
  call define_bc_ritske (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo) 

case ('input.ritske_isurf')
  call define_bc_ritske_isurf (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo) 

case ('input.subduction')
  call define_bc_subduction (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo) 

case ('input.riedel')
  call define_bc_riedel (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo) 

case ('input.parallipipede')
  call define_bc_parallipipede (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo) 

case default
  if (iproc.eq.0) print *,params%infile
  call stop_run (' pb with input file in define_bc.f90')

end select

! here we need to impose that the bc satisfy the bad faces constrains

do ie=1,osolve%nface
 do j=1,4
 jp=1+mod(j,4)
 i=osolve%iface(j+4,ie)
 if (osolve%kfix((i-1)*3+1).eq.1) osolve%u(i)=(osolve%u(osolve%iface(j,ie))+osolve%u(osolve%iface(jp,ie)))/2.d0
 if (osolve%kfix((i-1)*3+2).eq.1) osolve%v(i)=(osolve%v(osolve%iface(j,ie))+osolve%v(osolve%iface(jp,ie)))/2.d0
 if (osolve%kfix((i-1)*3+3).eq.1) osolve%w(i)=(osolve%w(osolve%iface(j,ie))+osolve%w(osolve%iface(jp,ie)))/2.d0
 if (osolve%kfixt(i).eq.1) osolve%temp(i)=(osolve%temp(osolve%iface(j,ie))+osolve%temp(osolve%iface(jp,ie)))/2.d0
 enddo
i=osolve%iface(9,ie)
if (osolve%kfix((i-1)*3+1).eq.1) osolve%u(i)=sum(osolve%u(osolve%iface(1:4,ie)))/4.d0
if (osolve%kfix((i-1)*3+2).eq.1) osolve%v(i)=sum(osolve%v(osolve%iface(1:4,ie)))/4.d0
if (osolve%kfix((i-1)*3+3).eq.1) osolve%w(i)=sum(osolve%w(osolve%iface(1:4,ie)))/4.d0
if (osolve%kfixt(i).eq.1) osolve%temp(i)=sum(osolve%temp(osolve%iface(1:4,ie)))/4.d0
enddo

if (params%debug==2) call output_bc (osolve)

return
end
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|