From 84de35e8c2110dce4ede774c6c6f5a440e097485 Mon Sep 17 00:00:00 2001 From: Dave Whipp <dwhipp@dal.ca> Date: Mon, 21 Dec 2009 14:59:08 +0000 Subject: [PATCH] added new isostasy B/C routine --- src/define_isostasy_bc.f90 | 157 +++++++++++++++++++++++++++++++++++++ 1 file changed, 157 insertions(+) create mode 100644 src/define_isostasy_bc.f90 diff --git a/src/define_isostasy_bc.f90 b/src/define_isostasy_bc.f90 new file mode 100644 index 00000000..65b26e4f --- /dev/null +++ b/src/define_isostasy_bc.f90 @@ -0,0 +1,157 @@ +!------------------------------------------------------------------------------| +!------------------------------------------------------------------------------| +! | +! ||===\\ | +! || \\ | +! || || //==\\ || || //==|| ||/==\\ | +! || || || || || || || || || || | +! || // || || || || || || || | +! ||===// \\==// \\==\\ \\==\\ || | +! | +!------------------------------------------------------------------------------| +!------------------------------------------------------------------------------| +! | +! DEFINE_ISOSTASY_BC OCT. 2009 | +! | +!------------------------------------------------------------------------------| +!------------------------------------------------------------------------------| + +subroutine define_isostasy_bc(params,osolve,vo,zisodisp,firstcall,l,x0,spu,spv,& + ldisp) + +!------------------------------------------------------------------------------| +!(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))) +!------------------------------------------------------------------------------| +! New subroutine to apply the modified isostasy velocity boundary conditions +! - This routine uses the basal displacement array defined in isostasy.f90 to +! modify the basal velocity boundary conditions to have the incoming flux +! velocity tangent to the moho. This might not be exactly correct, but works +! for now. +! +! dwhipp (11/09) +! +!------------------------------------------------------------------------------| +!(((((((((((((((( declaration of the subroutine arguments )))))))))))))))))))) +!------------------------------------------------------------------------------| + +use definitions + +implicit none + +type (parameters) params +type (octreesolve) osolve +type (void) vo +double precision zisodisp(2**params%levelmax_oct+1,2**params%levelmax_oct+1) +double precision x0(osolve%nnode),ldisp(osolve%nnode) +double precision l,spu,spv +logical firstcall + +!------------------------------------------------------------------------------| +!(((((((((((((((( declaration of the subroutine internal variables ))))))))))))) +!------------------------------------------------------------------------------| + +integer i,j,iproc,nproc,ierr,nb,xdisp,ydisp,xnow,ynow,ie,jp +double precision eps,pi,zsl,dxy,zdisp,vinit +double precision,dimension(:,:),allocatable::zisoslx,zisosly + +!------------------------------------------------------------------------------| +!------------------------------------------------------------------------------| + +INCLUDE 'mpif.h' + +call mpi_comm_size (mpi_comm_world,nproc,ierr) +call mpi_comm_rank (mpi_comm_world,iproc,ierr) + +! General variable definition +eps=1.d-10 ! Tiny number +pi=atan(1.d0)*4.d0 ! Pi +nb=2**params%levelmax_oct ! Number of nodes along one side of model +dxy=1./real(nb) ! Non-dimensional node spacing + +if (firstcall) then + do i=1,osolve%nnode + ! Shift s-line depth + xdisp=nint(x0(i)*nb)+1 + ydisp=nint(osolve%y(i)*nb)+1 + ldisp(i)=l+zisodisp(xdisp,ydisp) + enddo +else + allocate(zisoslx(nb+1,nb+1))!,zisosly(nb+1,nb+1)) + do j=1,nb+1 + do i=1,nb+1 + if (i==1) then + zisoslx(i,j)=(zisodisp(i+1,j)-zisodisp(i,j))/dxy ! Displacement between first and second nodes along x-axis + elseif (i==nb+1) then + zisoslx(i,j)=(zisodisp(i,j)-zisodisp(i-1,j))/dxy ! Displacement between second-to-last and last nodes along x-axis + else + zisoslx(i,j)=((zisodisp(i,j)-zisodisp(i-1,j))/dxy+& ! Average displacement of nodes on either side of current node, along x-axis + (zisodisp(i+1,j)-zisodisp(i,j))/dxy)/2.d0 + endif + !if (j==1) then + ! zisosly(i,j)=(zisodisp(i,j+1)-zisodisp(i,j))/dxy ! Displacement between first and second nodes along y-axis + !elseif (j==nb+1) then + ! zisosly(i,j)=(zisodisp(i,j)-zisodisp(i,j-1))/dxy ! Displacement between second-to-last and last nodes along y-axis + !else + ! zisosly(i,j)=((zisodisp(i,j)-zisodisp(i,j-1))/dxy+& ! Average displacement of nodes on either side of current node, along y-axis + ! (zisodisp(i,j+1)-zisodisp(i,j))/dxy)/2.d0 + !endif + enddo + enddo + + do i=1,osolve%nnode + xnow=nint(osolve%x(i)*nb)+1 + ynow=nint(osolve%y(i)*nb)+1 + zsl=zisoslx(xnow,ynow) + + ! Modify velocities + if (osolve%z(i).lt.eps) then +! if (osolve%x(i).le.x0) then +! if (spu.le.1.d0-eps) then ! No change in B/Cs + if (osolve%u(i).gt.eps) then + zdisp=osolve%u(i)*sin(atan(zsl)) ! Negative for negative slope + vinit=(osolve%u(i)**2.d0+osolve%w(i)**2.d0)**0.5d0 ! Incoming velocity vector magnitude + osolve%w(i)=osolve%w(i)+zdisp ! New z-velocity + osolve%u(i)=sqrt(vinit**2-osolve%w(i)**2.d0) ! New x-velocity + endif +! endif +! else +! if (spu.ge.eps) then +! zdisp=osolve%u(i)*sin(atan(zsl)) ! Negative for negative slope +! osolve%w(i)=zdisp-osolve%w(i) ! New z-velocity +! osolve%u(i)=sqrt(vin**2-osolve%w(i)**2.d0) ! New x-velocity +! endif +! endif + endif + enddo + deallocate(zisoslx) +endif + +! 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 + +return +end +!------------------------------------------------------------------------------| +!------------------------------------------------------------------------------| -- GitLab