diff --git a/DOUAR.f90 b/DOUAR.f90
index 68db5bdb7ad0c7ae45788f84d8255ea302d711bd..7b335de6051142c913d103db776c1c548e546bbf 100644
--- a/DOUAR.f90
+++ b/DOUAR.f90
@@ -308,14 +308,22 @@ do while (istep.le.params%nstep)
!---------------------------------------------------------------------------|
do is=1,params%ns
surface0(is)%nsurface=surface(is)%nsurface
- allocate (surface0(is)%r(surface0(is)%nsurface),stat=threadinfo%err) ; call heap (threadinfo,'surface0(is)%r','main',size(surface0(is)%r),'dp',+1)
- allocate (surface0(is)%s(surface0(is)%nsurface),stat=threadinfo%err) ; call heap (threadinfo,'surface0(is)%s','main',size(surface0(is)%s),'dp',+1)
- allocate (surface0(is)%x(surface0(is)%nsurface),stat=threadinfo%err) ; call heap (threadinfo,'surface0(is)%x','main',size(surface0(is)%x),'dp',+1)
- allocate (surface0(is)%y(surface0(is)%nsurface),stat=threadinfo%err) ; call heap (threadinfo,'surface0(is)%y','main',size(surface0(is)%y),'dp',+1)
- allocate (surface0(is)%z(surface0(is)%nsurface),stat=threadinfo%err) ; call heap (threadinfo,'surface0(is)%z','main',size(surface0(is)%z),'dp',+1)
- allocate (surface0(is)%xn(surface0(is)%nsurface),stat=threadinfo%err) ; call heap (threadinfo,'surface0(is)%xn','main',size(surface0(is)%xn),'dp',+1)
- allocate (surface0(is)%yn(surface0(is)%nsurface),stat=threadinfo%err) ; call heap (threadinfo,'surface0(is)%yn','main',size(surface0(is)%yn),'dp',+1)
- allocate (surface0(is)%zn(surface0(is)%nsurface),stat=threadinfo%err) ; call heap (threadinfo,'surface0(is)%zn','main',size(surface0(is)%zn),'dp',+1)
+ allocate (surface0(is)%r(surface0(is)%nsurface),stat=threadinfo%err)
+ call heap (threadinfo,'surface0(is)%r','main',size(surface0(is)%r),'dp',+1)
+ allocate (surface0(is)%s(surface0(is)%nsurface),stat=threadinfo%err)
+ call heap (threadinfo,'surface0(is)%s','main',size(surface0(is)%s),'dp',+1)
+ allocate (surface0(is)%x(surface0(is)%nsurface),stat=threadinfo%err)
+ call heap (threadinfo,'surface0(is)%x','main',size(surface0(is)%x),'dp',+1)
+ allocate (surface0(is)%y(surface0(is)%nsurface),stat=threadinfo%err)
+ call heap (threadinfo,'surface0(is)%y','main',size(surface0(is)%y),'dp',+1)
+ allocate (surface0(is)%z(surface0(is)%nsurface),stat=threadinfo%err)
+ call heap (threadinfo,'surface0(is)%z','main',size(surface0(is)%z),'dp',+1)
+ allocate (surface0(is)%xn(surface0(is)%nsurface),stat=threadinfo%err)
+ call heap (threadinfo,'surface0(is)%xn','main',size(surface0(is)%xn),'dp',+1)
+ allocate (surface0(is)%yn(surface0(is)%nsurface),stat=threadinfo%err)
+ call heap (threadinfo,'surface0(is)%yn','main',size(surface0(is)%yn),'dp',+1)
+ allocate (surface0(is)%zn(surface0(is)%nsurface),stat=threadinfo%err)
+ call heap (threadinfo,'surface0(is)%zn','main',size(surface0(is)%zn),'dp',+1)
surface0(is)%r=surface(is)%r
surface0(is)%s=surface(is)%s
surface0(is)%x=surface(is)%x
@@ -352,7 +360,8 @@ do while (istep.le.params%nstep)
write(*,'(a,L1)') shift//'(1) params%griditer < 0 ',(params%griditer.lt.0)
write(*,'(a,L1)') shift//'(2) current_level==params%levelmax_oct',(current_level==params%levelmax_oct)
write(*,'(a,L1)') shift//'(3) increase_current_level ',increase_current_level
- write(*,'(a,L1)') shift//'converge = (1) & (2) & (3) -> ',(params%griditer .lt. 0 .and. current_level==params%levelmax_oct .and. increase_current_level)
+ write(*,'(a,L1)') shift//'converge = (1) & (2) & (3) -> ',&
+ (params%griditer .lt. 0 .and. current_level==params%levelmax_oct .and. increase_current_level)
end if
@@ -374,7 +383,8 @@ do while (istep.le.params%nstep)
osolve%noctree=params%noctreemax
- allocate (osolve%octree(osolve%noctree),stat=threadinfo%err) ; call heap (threadinfo,'osolve%octree','main',size(osolve%octree),'int',+1)
+ allocate (osolve%octree(osolve%noctree),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%octree','main',size(osolve%octree),'int',+1)
call octree_init (osolve%octree,osolve%noctree)
call octree_create_uniform (osolve%octree,osolve%noctree,params%leveluniform_oct)
@@ -402,11 +412,16 @@ do while (istep.le.params%nstep)
if (iproc.eq.0) write (8,*) osolve%nleaves,' leaves in solve octree '
- allocate (osolve%icon(8,osolve%nleaves),stat=threadinfo%err) ; call heap (threadinfo,'osolve%icon','main',size(osolve%icon),'int',+1)
- allocate (osolve%x(osolve%nnode),stat=threadinfo%err) ; call heap (threadinfo,'osolve%x', 'main',size(osolve%x),'dp',+1)
- allocate (osolve%y(osolve%nnode),stat=threadinfo%err) ; call heap (threadinfo,'osolve%y', 'main',size(osolve%y),'dp',+1)
- allocate (osolve%z(osolve%nnode),stat=threadinfo%err) ; call heap (threadinfo,'osolve%z', 'main',size(osolve%z),'dp',+1)
- allocate (osolve%lsf(osolve%nnode,osolve%nlsf),stat=threadinfo%err) ; call heap (threadinfo,'osolve%lsf', 'main',size(osolve%lsf),'dp',+1)
+ allocate (osolve%icon(8,osolve%nleaves),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%icon','main',size(osolve%icon),'int',+1)
+ allocate (osolve%x(osolve%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%x', 'main',size(osolve%x),'dp',+1)
+ allocate (osolve%y(osolve%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%y', 'main',size(osolve%y),'dp',+1)
+ allocate (osolve%z(osolve%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%z', 'main',size(osolve%z),'dp',+1)
+ allocate (osolve%lsf(osolve%nnode,osolve%nlsf),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%lsf', 'main',size(osolve%lsf),'dp',+1)
osolve%lsf=0.d0
call octree_find_node_connectivity (osolve%octree,osolve%noctree, &
@@ -460,13 +475,20 @@ do while (istep.le.params%nstep)
! allocate fields of osolve
!------------------------------------------------------------------------|
!------------------------------------------------------------------------|
- allocate (osolve%u(osolve%nnode),stat=threadinfo%err) ; call heap (threadinfo,'osolve%u','main',size(osolve%u),'dp',+1)
- allocate (osolve%v(osolve%nnode),stat=threadinfo%err) ; call heap (threadinfo,'osolve%v','main',size(osolve%v),'dp',+1)
- allocate (osolve%w(osolve%nnode),stat=threadinfo%err) ; call heap (threadinfo,'osolve%w','main',size(osolve%w),'dp',+1)
- allocate (osolve%temp(osolve%nnode),stat=threadinfo%err) ; call heap (threadinfo,'osolve%temp','main',size(osolve%temp),'dp',+1)
- allocate (osolve%pressure(osolve%nleaves),stat=threadinfo%err) ; call heap (threadinfo,'osolve%pressure','main',size(osolve%pressure),'dp',+1)
- allocate (osolve%spressure(osolve%nleaves),stat=threadinfo%err) ; call heap (threadinfo,'osolve%spressure','main',size(osolve%spressure),'dp',+1)
- allocate (osolve%strain(osolve%nnode),stat=threadinfo%err) ; call heap (threadinfo,'osolve%strain','main',size(osolve%strain),'dp',+1)
+ allocate (osolve%u(osolve%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%u','main',size(osolve%u),'dp',+1)
+ allocate (osolve%v(osolve%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%v','main',size(osolve%v),'dp',+1)
+ allocate (osolve%w(osolve%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%w','main',size(osolve%w),'dp',+1)
+ allocate (osolve%temp(osolve%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%temp','main',size(osolve%temp),'dp',+1)
+ allocate (osolve%pressure(osolve%nleaves),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%pressure','main',size(osolve%pressure),'dp',+1)
+ allocate (osolve%spressure(osolve%nleaves),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%spressure','main',size(osolve%spressure),'dp',+1)
+ allocate (osolve%strain(osolve%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%strain','main',size(osolve%strain),'dp',+1)
osolve%u = 0.d0
osolve%v = 0.d0
osolve%w = 0.d0
@@ -505,7 +527,8 @@ do while (istep.le.params%nstep)
call show_time (total,step,inc,1,'Find bad faces$')
osolve%nface=osolve%nleaves
- allocate (osolve%iface(9,osolve%nface),stat=threadinfo%err) ; call heap (threadinfo,'osolve%iface','main',size(osolve%iface),'int',+1)
+ allocate (osolve%iface(9,osolve%nface),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%iface','main',size(osolve%iface),'int',+1)
call octree_find_bad_faces (osolve%octree,osolve%noctree, &
osolve%iface,osolve%nface, &
@@ -523,12 +546,18 @@ do while (istep.le.params%nstep)
!------------------------------------------------------------------------|
call show_time (total,step,inc,1,'Find void$')
- allocate (vo%node(osolve%nnode),stat=threadinfo%err) ; call heap (threadinfo,'vo%node','main',size(vo%node),'int',+1)
- allocate (vo%leaf(osolve%nnode),stat=threadinfo%err) ; call heap (threadinfo,'vo%leaf','main',size(vo%leaf),'int',+1)
- allocate (vo%face(osolve%nface),stat=threadinfo%err) ; call heap (threadinfo,'vo%face','main',size(vo%face),'int',+1)
- allocate (vo%rtf(osolve%nnode),stat=threadinfo%err) ; call heap (threadinfo,'vo%rtf','main',size(vo%rtf),'int',+1)
- allocate (vo%ftr(osolve%nnode),stat=threadinfo%err) ; call heap (threadinfo,'vo%ftr','main',size(vo%ftr),'int',+1)
- allocate (vo%influid(osolve%nnode),stat=threadinfo%err) ; call heap (threadinfo,'vo%influid','main',size(vo%influid),'bool',+1)
+ allocate (vo%node(osolve%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'vo%node','main',size(vo%node),'int',+1)
+ allocate (vo%leaf(osolve%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'vo%leaf','main',size(vo%leaf),'int',+1)
+ allocate (vo%face(osolve%nface),stat=threadinfo%err)
+ call heap (threadinfo,'vo%face','main',size(vo%face),'int',+1)
+ allocate (vo%rtf(osolve%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'vo%rtf','main',size(vo%rtf),'int',+1)
+ allocate (vo%ftr(osolve%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'vo%ftr','main',size(vo%ftr),'int',+1)
+ allocate (vo%influid(osolve%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'vo%influid','main',size(vo%influid),'bool',+1)
call find_void_nodes (params,osolve,vo,istep,ov,threadinfo)
@@ -538,8 +567,10 @@ do while (istep.le.params%nstep)
!------------------------------------------------------------------------|
!------------------------------------------------------------------------|
call show_time (total,step,inc,1,'Define boundary conditions$')
- allocate (osolve%kfix(osolve%nnode*3),stat=threadinfo%err) ; call heap (threadinfo,'osolve%kfix','main',size(osolve%kfix),'int',+1)
- allocate (osolve%kfixt(osolve%nnode),stat=threadinfo%err) ; call heap (threadinfo,'osolve%kfixt','main',size(osolve%kfixt),'int',+1)
+ allocate (osolve%kfix(osolve%nnode*3),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%kfix','main',size(osolve%kfix),'int',+1)
+ allocate (osolve%kfixt(osolve%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%kfixt','main',size(osolve%kfixt),'int',+1)
call define_bc (params,osolve,vo)
@@ -588,14 +619,22 @@ do while (istep.le.params%nstep)
!------------------------------------------------------------------------|
!------------------------------------------------------------------------|
- allocate (osolve%eviscosity(osolve%nleaves),stat=threadinfo%err) ; call heap (threadinfo,'osolve%eviscosity','main',size(osolve%eviscosity),'dp',+1)
- allocate (osolve%q(osolve%nleaves),stat=threadinfo%err) ; call heap (threadinfo,'osolve%q','main',size(osolve%q),'dp',+1)
- allocate (osolve%crit(osolve%nleaves),stat=threadinfo%err) ; call heap (threadinfo,'osolve%crit','main',size(osolve%crit),'dp',+1)
- allocate (osolve%e2d(osolve%nleaves),stat=threadinfo%err) ; call heap (threadinfo,'osolve%e2d','main',size(osolve%e2d),'dp',+1)
- allocate (osolve%e3d(osolve%nleaves),stat=threadinfo%err) ; call heap (threadinfo,'osolve%e3d','main',size(osolve%e3d),'dp',+1)
- allocate (osolve%lode(osolve%nleaves),stat=threadinfo%err) ; call heap (threadinfo,'osolve%lode','main',size(osolve%lode),'dp',+1)
- allocate (osolve%is_plastic(osolve%nleaves),stat=threadinfo%err) ; call heap (threadinfo,'osolve%is_plastic','main',size(osolve%is_plastic),'bool',+1)
- allocate (weightel(osolve%nleaves)) ; call heap (threadinfo,'weightel','main',size(weightel),'bool',+1)
+ allocate (osolve%eviscosity(osolve%nleaves),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%eviscosity','main',size(osolve%eviscosity),'dp',+1)
+ allocate (osolve%q(osolve%nleaves),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%q','main',size(osolve%q),'dp',+1)
+ allocate (osolve%crit(osolve%nleaves),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%crit','main',size(osolve%crit),'dp',+1)
+ allocate (osolve%e2d(osolve%nleaves),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%e2d','main',size(osolve%e2d),'dp',+1)
+ allocate (osolve%e3d(osolve%nleaves),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%e3d','main',size(osolve%e3d),'dp',+1)
+ allocate (osolve%lode(osolve%nleaves),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%lode','main',size(osolve%lode),'dp',+1)
+ allocate (osolve%is_plastic(osolve%nleaves),stat=threadinfo%err)
+ call heap (threadinfo,'osolve%is_plastic','main',size(osolve%is_plastic),'bool',+1)
+ allocate (weightel(osolve%nleaves))
+ call heap (threadinfo,'weightel','main',size(weightel),'bool',+1)
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
@@ -731,8 +770,8 @@ do while (istep.le.params%nstep)
call show_time (total,step,inc,1,'slicing the cube$')
call slices(params,osolve,ov,vo,sections,istep,iter,iter_nl)
- call heap (threadinfo,'osolve%eviscosity','main',size(osolve%eviscosity),'dp',-1) ; deallocate (osolve%eviscosity)
- call heap (threadinfo,'osolve%q','main',size(osolve%q),'dp',-1) ; deallocate (osolve%q)
+ call heap (threadinfo,'osolve%eviscosity','main',size(osolve%eviscosity),'dp',-1) ; deallocate (osolve%eviscosity)
+ call heap (threadinfo,'osolve%q','main',size(osolve%q),'dp',-1) ; deallocate (osolve%q)
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! compute isostasy
@@ -784,7 +823,7 @@ do while (istep.le.params%nstep)
if (.not.converge) then
call heap (threadinfo,'osolve%octree','main',size(osolve%octree),'int',-1) ; deallocate (osolve%octree)
- call heap (threadinfo,'osolve%icon','main',size(osolve%icon),'int',-1) ; deallocate (osolve%icon)
+ call heap (threadinfo,'osolve%icon','main',size(osolve%icon),'int',-1) ; deallocate (osolve%icon)
call heap (threadinfo,'osolve%x','main',size(osolve%x),'dp',-1) ; deallocate (osolve%x)
call heap (threadinfo,'osolve%y','main',size(osolve%y),'dp',-1) ; deallocate (osolve%y)
call heap (threadinfo,'osolve%z','main',size(osolve%z),'dp',-1) ; deallocate (osolve%z)
@@ -858,12 +897,18 @@ do while (istep.le.params%nstep)
olsf%nleaves=ov%nleaves
olsf%noctree=ov%noctree
- allocate (olsf%octree(olsf%noctree),stat=threadinfo%err) ; call heap (threadinfo,'olsf%octree','main',size(olsf%octree),'int',+1)
- allocate (olsf%x(olsf%nnode),stat=threadinfo%err) ; call heap (threadinfo,'olsf%x','main',size(olsf%x),'dp',+1)
- allocate (olsf%y(olsf%nnode),stat=threadinfo%err) ; call heap (threadinfo,'olsf%y','main',size(olsf%y),'dp',+1)
- allocate (olsf%z(olsf%nnode),stat=threadinfo%err) ; call heap (threadinfo,'olsf%z','main',size(olsf%z),'dp',+1)
- allocate (olsf%lsf(olsf%nnode),stat=threadinfo%err) ; call heap (threadinfo,'olsf%lsf','main',size(olsf%lsf),'dp',+1)
- allocate (olsf%icon(params%mpe,olsf%nleaves),stat=threadinfo%err) ; call heap (threadinfo,'olsf%icon','main',size(olsf%icon),'int',+1)
+ allocate (olsf%octree(olsf%noctree),stat=threadinfo%err)
+ call heap (threadinfo,'olsf%octree','main',size(olsf%octree),'int',+1)
+ allocate (olsf%x(olsf%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'olsf%x','main',size(olsf%x),'dp',+1)
+ allocate (olsf%y(olsf%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'olsf%y','main',size(olsf%y),'dp',+1)
+ allocate (olsf%z(olsf%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'olsf%z','main',size(olsf%z),'dp',+1)
+ allocate (olsf%lsf(olsf%nnode),stat=threadinfo%err)
+ call heap (threadinfo,'olsf%lsf','main',size(olsf%lsf),'dp',+1)
+ allocate (olsf%icon(params%mpe,olsf%nleaves),stat=threadinfo%err)
+ call heap (threadinfo,'olsf%icon','main',size(olsf%icon),'int',+1)
olsf%octree(1:olsf%noctree)=ov%octree(1:ov%noctree)
olsf%x(1:olsf%nnode)=ov%x(1:ov%nnode)
@@ -1079,12 +1124,13 @@ do while (istep.le.params%nstep)
olsf%nleaves=ov%nleaves
olsf%noctree=ov%noctree
- allocate (olsf%octree(olsf%noctree),stat=threadinfo%err) ; call heap (threadinfo,'olsf%octree','main',size(olsf%octree),'int',+1)
- allocate (olsf%x(olsf%nnode),stat=threadinfo%err) ; call heap (threadinfo,'olsf%x','main',size(olsf%x),'dp',+1)
- allocate (olsf%y(olsf%nnode),stat=threadinfo%err) ; call heap (threadinfo,'olsf%y','main',size(olsf%y),'dp',+1)
- allocate (olsf%z(olsf%nnode),stat=threadinfo%err) ; call heap (threadinfo,'olsf%z','main',size(olsf%z),'dp',+1)
- allocate (olsf%lsf(olsf%nnode),stat=threadinfo%err) ; call heap (threadinfo,'olsf%lsf','main',size(olsf%lsf),'dp',+1)
- allocate (olsf%icon(params%mpe,olsf%nleaves),stat=threadinfo%err) ; call heap (threadinfo,'olsf%icon','main',size(olsf%icon),'int',+1)
+ allocate (olsf%octree(olsf%noctree),stat=threadinfo%err) ; call heap(threadinfo,'olsf%octree','main',size(olsf%octree),'int',+1)
+ allocate (olsf%x(olsf%nnode),stat=threadinfo%err) ; call heap (threadinfo,'olsf%x','main',size(olsf%x),'dp',+1)
+ allocate (olsf%y(olsf%nnode),stat=threadinfo%err) ; call heap (threadinfo,'olsf%y','main',size(olsf%y),'dp',+1)
+ allocate (olsf%z(olsf%nnode),stat=threadinfo%err) ; call heap (threadinfo,'olsf%z','main',size(olsf%z),'dp',+1)
+ allocate (olsf%lsf(olsf%nnode),stat=threadinfo%err) ; call heap (threadinfo,'olsf%lsf','main',size(olsf%lsf),'dp',+1)
+ allocate (olsf%icon(params%mpe,olsf%nleaves),stat=threadinfo%err)
+ call heap (threadinfo,'olsf%icon','main',size(olsf%icon),'int',+1)
olsf%octree(1:olsf%noctree)=ov%octree(1:ov%noctree)
olsf%x(1:olsf%nnode)=ov%x(1:ov%nnode)
@@ -1258,8 +1304,10 @@ call heap (threadinfo,'ov%vnode','main',size(ov%vnode),'dp',-1) ; deallocat
call heap (threadinfo,'ov%wnode','main',size(ov%wnode),'dp',-1) ; deallocate (ov%wnode)
call heap (threadinfo,'ov%icon','main',size(ov%icon),'int',-1) ; deallocate (ov%icon)
call heap (threadinfo,'ov%temp','main',size(ov%temp),'dp',-1) ; deallocate (ov%temp)
-call heap (threadinfo,'ov%temporary_nodal...','main',size(ov%temporary_nodal_pressure),'dp',-1) ; deallocate (ov%temporary_nodal_pressure)
-call heap (threadinfo,'ov%whole_leaf...','main',size(ov%whole_leaf_in_fluid),'bool',-1) ; deallocate (ov%whole_leaf_in_fluid)
+call heap (threadinfo,'ov%temporary_nodal...','main',size(ov%temporary_nodal_pressure),'dp',-1)
+deallocate (ov%temporary_nodal_pressure)
+call heap (threadinfo,'ov%whole_leaf...','main',size(ov%whole_leaf_in_fluid),'bool',-1)
+deallocate (ov%whole_leaf_in_fluid)
do is=1,params%ns