isostasy.f90

subroutine isostasy (params,weightel,ov,surface,mat)

! 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)

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

! 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%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