From fd6da22a36a85e62fce43892701594f1cdf77b4e Mon Sep 17 00:00:00 2001
From: Dave Whipp <dwhipp@dal.ca>
Date: Mon, 21 Dec 2009 15:39:38 +0000
Subject: [PATCH] Added storage of zisodisp array
---
src/isostasy.f90 | 186 ++++++++++++++++++++++++-----------------------
1 file changed, 95 insertions(+), 91 deletions(-)
diff --git a/src/isostasy.f90 b/src/isostasy.f90
index 267f82bd..31c69c75 100644
--- a/src/isostasy.f90
+++ b/src/isostasy.f90
@@ -1,4 +1,4 @@
-subroutine isostasy (params,weightel,ov,surface,mat)
+subroutine isostasy (params,weightel,ov,surface,mat,zisodisp)
! first go at imposing isostasy
! this routine assumes that the original surface are flat
@@ -15,6 +15,7 @@ type(octreev) ov
double precision weightel(ov%nleaves)
type(sheet) surface(params%ns)
type(material) mat(0:params%nmat)
+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
@@ -107,7 +108,7 @@ 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
+!if (iproc.eq.0) print*,'hx,dh,xk,d',hx,dh,xk,d
ii=nflex/2-nb/2-1
jj=ii
@@ -118,7 +119,7 @@ flex=0.d0
enddo
enddo
-if (iproc.eq.0) print*,'poids',minval(flex),maxval(flex)
+!if (iproc.eq.0) print*,'poids',minval(flex),maxval(flex)
do j=1,nflex
call dsinft (flex(1,j),nflex)
@@ -165,7 +166,7 @@ pihx=pi/hx
call dsinft (flex(1,j),nflex)
enddo
-if (iproc.eq.0) print*,'def in m',minval(flex),maxval(flex)
+!if (iproc.eq.0) print*,'def in m',minval(flex),maxval(flex)
ii=nflex/2-nb/2-1
jj=ii
@@ -177,7 +178,7 @@ jj=ii
deallocate (flex,work)
-if (iproc.eq.0) print*,'dimensionless def ',minval(weight),maxval(weight)
+!if (iproc.eq.0) print*,'dimensionless def ',minval(weight),maxval(weight)
endif
@@ -206,59 +207,62 @@ 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
+!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
+!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) zisodisp=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
@@ -267,47 +271,47 @@ np=0
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
+!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)
--
GitLab