subroutine isostasy (params,weightel,ov,surface,mat,iter,zi)
! first go at imposing isostasy
! this routine assumes that the original surface are flat
! it takes the zref height of each surface and the density of the corresponding material
! to compute the reference isostatic column
! note also that this will not work when the courant condition is used
use definitions
implicit none
type(parameters) params
type(octreev) ov
double precision weightel(ov%nleaves)
type(sheet) surface(params%ns)
type(material) mat(0:params%nmat)
integer iter
type(ziso) zi
!double precision zisodisp(2**params%levelmax_oct+1,2**params%levelmax_oct+1)
double precision xc,yc,zc,eps
integer i,j,ii,jj,kk,ie,dlev,is,iproc,nproc,ierr,ic
double precision x0,y0,z0,dxyz,dd
integer leaf,level,loc
double precision x0b,y0b,z0b,dxyzb
integer leafb,levelb,locb
double precision,dimension(:),allocatable::diso
double precision weightref,disoel,w
integer,dimension(:),allocatable::np
character shift*72
double precision,dimension(:,:),allocatable::weight,flex,work
double precision pi,hx,dh,pihx,xk,fi,fj,tij,d
integer nflex
integer im,jm,ip,jp,nb
integer iunit
INCLUDE 'mpif.h'
call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)
eps=1.d-8
shift=' '
nb=2**params%levelmax_oct
allocate (weight(nb,nb))
! first calculates the weight of each column, based on the weight of each elements
! this is performed on a 2**levelmax X 2**levelmax grid at the base of the model
! the result is stored in weight(i,j)
dd=1.d0/2.d0**(params%levelmax_oct+1)
weight=0.d0
do ie=1,ov%nleaves
xc=0.d0;yc=0.d0;zc=0.d0
do i=1,8
xc=xc+ov%x(ov%icon(i,ie))
yc=yc+ov%y(ov%icon(i,ie))
zc=zc+ov%z(ov%icon(i,ie))
enddo
xc=xc/8.d0
yc=yc/8.d0
zc=zc/8.d0
call octree_find_leaf (ov%octree,ov%noctree,xc,yc,zc,leaf,level,loc,x0,y0,z0,dxyz)
dlev=params%levelmax_oct-level
call find_integer_coordinates (x0+dd,y0+dd,z0+dd,ii,jj,kk,params%levelmax_oct)
ii=min(ii,nb-2**dlev)
jj=min(jj,nb-2**dlev)
do i=1,2**dlev
do j=1,2**dlev
weight(ii+i,jj+j)=weight(ii+i,jj+j)+weightel(leaf)/(4.d0**dlev)*(4.d0**params%levelmax_oct)
enddo
enddo
! call octree_find_leaf (ov%octree,ov%noctree,xc,yc,eps,leafb,levelb,locb,x0b,y0b,z0b,dxyzb)
! dlev=max(0,levelb-level)
! w(leafb)=w(leafb)+weightel(leaf)/(4.d0**dlev)*(4.d0**levelb)
enddo
! calculates the reference weight from the initial surface geometry and the density of the material
! near the base
weightref=0.d0
do is=1,params%ns-1
weightref=weightref+mat(surface(is)%material)%density*(surface(is)%sp01-surface(is+1)%sp01)
enddo
weightref=weightref+mat(surface(params%ns)%material)%density*surface(params%ns)%sp01
if (.not.params%flexure) then
! apply local isostasy
! the deflection is stored in array weight
do j=1,nb
do i=1,nb
weight(i,j)=(weight(i,j)-weightref)/mat(surface(params%ns)%material)%density
enddo
enddo
else
! flexural isostasy
nflex=8*nb
allocate(work(nflex,nflex),flex(nflex,nflex))
hx=params%length_scale*8.d3
xk=-mat(surface(params%ns)%material)%density*params%density_difference*9.81d0
d=1.d11/12.d0/(1.d0-0.25d0**2)*params%elastic_plate_thickness**3
dh=hx/(nflex-1)
!if (iproc.eq.0) print*,'hx,dh,xk,d',hx,dh,xk,d
ii=nflex/2-nb/2-1
jj=ii
flex=0.d0
do j=1,nb
do i=1,nb
flex(ii+i,jj+j)=(weight(i,j)-weightref)*dh*dh*params%density_difference*9.81d0*params%length_scale*1.d3
enddo
enddo
!if (iproc.eq.0) print*,'poids',minval(flex),maxval(flex)
do j=1,nflex
call dsinft (flex(1,j),nflex)
enddo
do j=1,nflex
do i=1,nflex
work(j,i)=flex(i,j)
enddo
enddo
do i=1,nflex
call dsinft (work(1,i),nflex)
enddo
do j=1,nflex
do i=1,nflex
work(j,i)=work(j,i)*4.d0/hx/hx
enddo
enddo
pi=atan(1.d0)*4.d0
pihx=pi/hx
do j=1,nflex
fj=(j*pihx)**2
do i=1,nflex
fi=(i*pihx)**2
tij=d/xk*(fi**2+2.d0*fi*fj+fj**2)+1.d0
work(j,i)=work(j,i)/xk/tij
enddo
enddo
do i=1,nflex
call dsinft (work(1,i),nflex)
enddo
do j=1,nflex
do i=1,nflex
flex(i,j)=work(j,i)
enddo
enddo
do j=1,nflex
call dsinft (flex(1,j),nflex)
enddo
!if (iproc.eq.0) print*,'def in m',minval(flex),maxval(flex)
ii=nflex/2-nb/2-1
jj=ii
do j=1,nb
do i=1,nb
weight(i,j)=-flex(ii+i,jj+j)/params%length_scale/1.d3
enddo
enddo
deallocate (flex,work)
!if (iproc.eq.0) print*,'dimensionless def ',minval(weight),maxval(weight)
endif
! redistributes the isostatic response onto the elements
! and distributes it onto the nodes into array diso(nnode)
allocate (np(ov%nnode))
allocate (diso(ov%nnode),work(nb+1,nb+1))
do j=1,nb+1
do i=1,nb+1
im=max(i-1,1)
jm=max(j-1,1)
ip=min(i,nb)
jp=min(j,nb)
work(i,j)=(weight(ip,jp)+weight(im,jp)+weight(ip,jm)+weight(im,jm))/4.d0
enddo
enddo
do i=1,ov%nnode
call find_integer_coordinates(ov%x(i)+1.d-20,ov%y(i)+1.d-20,ov%z(i)+1.d-20,ii,jj,kk,params%levelmax_oct)
ii=ii+1
jj=jj+1
kk=kk+1
if ((ii-1)*(ii-nb-1).gt.0 .or. (jj-1)*(jj-nb-1).gt.0) stop 'error in isostasy'
diso(i)=work(ii,jj)
enddo
!goto 3333
! old version : to be neglected new is much faster
!diso=0.d0
!np=0
! do ie=1,ov%nleaves
! xc=0.d0;yc=0.d0;zc=0.d0
! do i=1,8
! xc=xc+ov%x(ov%icon(i,ie))
! yc=yc+ov%y(ov%icon(i,ie))
! zc=zc+ov%z(ov%icon(i,ie))
! enddo
! xc=xc/8.d0
! yc=yc/8.d0
! zc=zc/8.d0
! call octree_find_leaf (ov%octree,ov%noctree,xc,yc,zc,leaf,level,loc,x0,y0,z0,dxyz)
!! call octree_find_leaf (ov%octree,ov%noctree,xc,yc,eps,leafb,levelb,locb,x0b,y0b,z0b,dxyzb)
!! disoel=(w(leafb)-weightref)/(mat(surface(params%ns)%material)%density)
! w=0.d0
! dlev=params%levelmax_oct-level
! call find_integer_coordinates (x0+dd,y0+dd,z0+dd,ii,jj,kk,params%levelmax_oct)
! ii=min(ii,2**(params%levelmax_oct)-2**dlev)
! jj=min(jj,2**(params%levelmax_oct)-2**dlev)
! do i=1,2**dlev
! do j=1,2**dlev
! w=w+weight(ii+i,jj+j)
! enddo
! enddo
! disoel=w/(2**dlev*2.d0**dlev)
!! disoel=(w-weightref)/(mat(surface(params%ns)%material)%density)
! do i=1,8
! ic=ov%icon(i,ie)
! np(ic)=np(ic)+1
! diso(ic)=diso(ic)+disoel
! enddo
! enddo
!
! do i=1,ov%nnode
! if (np(i).ne.0) then
! diso(i)=diso(i)/np(i)
! else
! diso(i)=0
! print*,'error np',i,ov%x(i),ov%y(i),ov%z(i)
! endif
! enddo
!
!3333 continue
if (params%debug>=1.and.iproc.eq.0) then
write(*,'(a,f10.7,1x,f10.7)') shift//'Isostatic velocity range = ',minval(diso)/params%dt,maxval(diso)/params%dt
endif
! Store the total displacement of the basal node plane
if (params%isobc .and. iter.eq.1) zi%zisodisp=zi%zisodisp-work
! transforms the displacements into velocities and add this contribution to the velocity
! solution of the Stokes equation before interpolation onto the surfaces and their advection
do i=1,ov%nnode
ov%wnodepreiso(i)=ov%wnode(i)
ov%wnode(i)=ov%wnode(i)-diso(i)/params%dt
enddo
!goto 1111
!
!! display results into a vtk file for debugging
!
!iunit=30
!open(unit=iunit,file='isostasy.vtk')
!write(iunit,'(a)')'# vtk DataFile Version 3.0'
!write(iunit,'(a)')'velocities'
!write(iunit,'(a)')'ASCII'
!write(iunit,'(a)')'DATASET UNSTRUCTURED_GRID'
!write(iunit,'(a7,i10,a6)')'POINTS ',ov%nnode,' float'
!
!do i=1,ov%nnode
!write(iunit,'(3f16.11)')ov%x(i),ov%y(i),ov%z(i)
!enddo
!
!write(iunit,'(A6, 2I10)') 'CELLS ',ov%nleaves,(1+8)*ov%nleaves
!do ie=1,ov%nleaves
!write(iunit,'(9I10)')8,ov%icon(1:8,ie)-1
!enddo
!
!write(iunit,'(A11, I10)') 'CELL_TYPES ',ov%nleaves
!do ie=1,ov%nleaves
! write(iunit,'(I2)')11 ! octree (8 nodes)
!end do
!
!write(iunit,'(a11,i10)')'POINT_DATA ',ov%nnode
!
!write(iunit,'(a)')'VECTORS iso float'
!do i=1,ov%nnode
! write(iunit,'(3e11.4)')0.,0.,-diso(i)/params%dt
!enddo
!
!write(iunit,'(a)')'VECTORS velo float'
!do i=1,ov%nnode
! write(iunit,'(3e11.4)')0.,0.,ov%wnode(i)
!enddo
!
!close (iunit)
!
!1111 continue
deallocate (np)
deallocate (diso)
deallocate (weight)
deallocate (work)
return
end subroutine isostasy
! SINFT
! taken from numerical recipes (see book for further information)
! subroutines called:
! NONE
SUBROUTINE DSINFT(Y,N)
common /vocal/ ivocal
REAL*8 WR,WI,WPR,WPI,WTEMP,THETA
double precision Y(N)
THETA=3.14159265358979D0/DBLE(N)
WR=1.0D0
WI=0.0D0
WPR=-2.0D0*DSIN(0.5D0*THETA)**2
WPI=DSIN(THETA)
Y(1)=0.0
M=N/2
DO 11 J=1,M
WTEMP=WR
WR=WR*WPR-WI*WPI+WR
WI=WI*WPR+WTEMP*WPI+WI
Y1=WI*(Y(J+1)+Y(N-J+1))
Y2=0.5*(Y(J+1)-Y(N-J+1))
Y(J+1)=Y1+Y2
Y(N-J+1)=Y1-Y2
11 CONTINUE
CALL DREALFT(Y,M,+1)
SUM=0.0
Y(1)=0.5*Y(1)
Y(2)=0.0
DO 12 J=1,N-1,2
SUM=SUM+Y(J)
Y(J)=Y(J+1)
Y(J+1)=SUM
12 CONTINUE
RETURN
END
! REALFT
! taken from numerical recipes (see book for further information)
! subroutines called:
! NONE
SUBROUTINE DREALFT(DATA,N,ISIGN)
common /vocal/ ivocal
REAL*8 WR,WI,WPR,WPI,WTEMP,THETA
double precision DATA(*)
THETA=6.28318530717959D0/2.0D0/DBLE(N)
C1=0.5
IF (ISIGN.EQ.1) THEN
C2=-0.5
CALL DFOUR1(DATA,N,+1)
ELSE
C2=0.5
THETA=-THETA
ENDIF
WPR=-2.0D0*DSIN(0.5D0*THETA)**2
WPI=DSIN(THETA)
WR=1.0D0+WPR
WI=WPI
N2P3=2*N+3
DO 11 I=2,N/2+1
I1=2*I-1
I2=I1+1
I3=N2P3-I2
I4=I3+1
WRS=SNGL(WR)
WIS=SNGL(WI)
H1R=C1*(DATA(I1)+DATA(I3))
H1I=C1*(DATA(I2)-DATA(I4))
H2R=-C2*(DATA(I2)+DATA(I4))
H2I=C2*(DATA(I1)-DATA(I3))
DATA(I1)=H1R+WRS*H2R-WIS*H2I
DATA(I2)=H1I+WRS*H2I+WIS*H2R
DATA(I3)=H1R-WRS*H2R+WIS*H2I
DATA(I4)=-H1I+WRS*H2I+WIS*H2R
WTEMP=WR
WR=WR*WPR-WI*WPI+WR
WI=WI*WPR+WTEMP*WPI+WI
11 CONTINUE
IF (ISIGN.EQ.1) THEN
H1R=DATA(1)
DATA(1)=H1R+DATA(2)
DATA(2)=H1R-DATA(2)
ELSE
H1R=DATA(1)
DATA(1)=C1*(H1R+DATA(2))
DATA(2)=C1*(H1R-DATA(2))
CALL DFOUR1(DATA,N,-1)
ENDIF
RETURN
END
! FOUR1
! this routine is taken directly out of numerical recipes
! see the book for further information
! subroutines called:
! NONE
SUBROUTINE DFOUR1(DATA,NN,ISIGN)
common /vocal/ ivocal
REAL*8 WR,WI,WPR,WPI,WTEMP,THETA
double precision DATA(*)
N=2*NN
J=1
DO 11 I=1,N,2
IF(J.GT.I)THEN
TEMPR=DATA(J)
TEMPI=DATA(J+1)
DATA(J)=DATA(I)
DATA(J+1)=DATA(I+1)
DATA(I)=TEMPR
DATA(I+1)=TEMPI
ENDIF
M=N/2
1 IF ((M.GE.2).AND.(J.GT.M)) THEN
J=J-M
M=M/2
GO TO 1
ENDIF
J=J+M
11 CONTINUE
MMAX=2
2 IF (N.GT.MMAX) THEN
ISTEP=2*MMAX
THETA=6.28318530717959D0/(ISIGN*MMAX)
WPR=-2.D0*DSIN(0.5D0*THETA)**2
WPI=DSIN(THETA)
WR=1.D0
WI=0.D0
DO 13 M=1,MMAX,2
DO 12 I=M,N,ISTEP
J=I+MMAX
TEMPR=SNGL(WR)*DATA(J)-SNGL(WI)*DATA(J+1)
TEMPI=SNGL(WR)*DATA(J+1)+SNGL(WI)*DATA(J)
DATA(J)=DATA(I)-TEMPR
DATA(J+1)=DATA(I+1)-TEMPI
DATA(I)=DATA(I)+TEMPR
DATA(I+1)=DATA(I+1)+TEMPI
12 CONTINUE
WTEMP=WR
WR=WR*WPR-WI*WPI+WR
WI=WI*WPR+WTEMP*WPI+WI
13 CONTINUE
MMAX=ISTEP
GO TO 2
ENDIF
RETURN
END