Skip to content
Snippets Groups Projects
Commit 120d3a4b authored by Dave Whipp's avatar Dave Whipp
Browse files

Modified allocation/passing of pre-isostasy vertical velocity array (wpreiso)

parent 2936d3ae
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/DOUAR.f90
+ 15
− 12
View file @ 120d3a4b
...@@ -478,7 +478,7 @@ do while (istep.le.params%nstep) ...@@ -478,7 +478,7 @@ do while (istep.le.params%nstep)
write(*,'(a)') shift//'After surfaces embedding:' write(*,'(a)') shift//'After surfaces embedding:'
write(*,'(a,i8,a)') shift//'nleaves_old= ',nleaves_old_mem2,' leaves' write(*,'(a,i8,a)') shift//'nleaves_old= ',nleaves_old_mem2,' leaves'
write(*,'(a,i8,a)') shift//'nleaves_new= ',nleaves_new_mem2,' leaves' write(*,'(a,i8,a)') shift//'nleaves_new= ',nleaves_new_mem2,' leaves'
write(*,'(a,l1)') shift//'C2: authorise increase of refine level = ',increase_current_level write(*,'(a,l1)') shift//'C2: authorise increase of refine level = ',increase_current_level
end if end if
nleaves_old_mem1=nleaves_new_mem1 nleaves_old_mem1=nleaves_new_mem1
...@@ -678,6 +678,7 @@ do while (istep.le.params%nstep) ...@@ -678,6 +678,7 @@ do while (istep.le.params%nstep)
osolve%u(i)=ov%unode(i) osolve%u(i)=ov%unode(i)
osolve%v(i)=ov%vnode(i) osolve%v(i)=ov%vnode(i)
osolve%w(i)=ov%wnode(i) osolve%w(i)=ov%wnode(i)
osolve%wpreiso(i)=ov%wnodepreiso(i)
enddo enddo
end if end if
! call define_bc (params,osolve,vo) ! call define_bc (params,osolve,vo)
...@@ -794,10 +795,11 @@ do while (istep.le.params%nstep) ...@@ -794,10 +795,11 @@ do while (istep.le.params%nstep)
!-------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------
call show_time (total,step,inc,1,'Compute isostasy and adjust vertical velocity$') if (params%isostasy) then
call show_time (total,step,inc,1,'Compute isostasy and adjust vertical velocity$')
allocate(ov%wpreiso(ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%wpreiso in main$') !allocate(ov%wpreiso(ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%wpreiso in main$')
if (params%isostasy) call isostasy (params,weightel,ov,surface,mat,iter,zi) call isostasy (params,weightel,ov,surface,mat,iter,zi)
endif
!-------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------
...@@ -853,7 +855,8 @@ do while (istep.le.params%nstep) ...@@ -853,7 +855,8 @@ do while (istep.le.params%nstep)
call heap (threadinfo,'osolve%kfixt','main',size(osolve%kfixt),'int',-1) ; deallocate (osolve%kfixt) call heap (threadinfo,'osolve%kfixt','main',size(osolve%kfixt),'int',-1) ; deallocate (osolve%kfixt)
call heap (threadinfo,'osolve%u','main',size(osolve%u),'dp',-1) ; deallocate (osolve%u) call heap (threadinfo,'osolve%u','main',size(osolve%u),'dp',-1) ; deallocate (osolve%u)
call heap (threadinfo,'osolve%v','main',size(osolve%v),'dp',-1) ; deallocate (osolve%v) call heap (threadinfo,'osolve%v','main',size(osolve%v),'dp',-1) ; deallocate (osolve%v)
call heap (threadinfo,'osolve%w','main',size(osolve%w),'dp',-1) ; deallocate (osolve%w) call heap (threadinfo,'osolve%w','main',size(osolve%w),'dp',-1) ; deallocate (osolve%w)
call heap (threadinfo,'osolve%wpreiso','main',size(osolve%wpreiso),'dp',-1) ; deallocate (osolve%wpreiso)
call heap (threadinfo,'osolve%temp','main',size(osolve%temp),'dp',-1) ; deallocate (osolve%temp) call heap (threadinfo,'osolve%temp','main',size(osolve%temp),'dp',-1) ; deallocate (osolve%temp)
call heap (threadinfo,'osolve%crip','main',size(osolve%crit),'dp',-1) ; deallocate (osolve%crit) call heap (threadinfo,'osolve%crip','main',size(osolve%crit),'dp',-1) ; deallocate (osolve%crit)
call heap (threadinfo,'osolve%e2d','main',size(osolve%e2d),'dp',-1) ; deallocate (osolve%e2d) call heap (threadinfo,'osolve%e2d','main',size(osolve%e2d),'dp',-1) ; deallocate (osolve%e2d)
...@@ -861,7 +864,7 @@ do while (istep.le.params%nstep) ...@@ -861,7 +864,7 @@ do while (istep.le.params%nstep)
call heap (threadinfo,'osolve%lode','main',size(osolve%lode),'dp',-1) ; deallocate (osolve%lode) call heap (threadinfo,'osolve%lode','main',size(osolve%lode),'dp',-1) ; deallocate (osolve%lode)
call heap (threadinfo,'osolve%eviscosity','main',size(osolve%eviscosity),'dp',-1) ; deallocate (osolve%eviscosity) call heap (threadinfo,'osolve%eviscosity','main',size(osolve%eviscosity),'dp',-1) ; deallocate (osolve%eviscosity)
call heap (threadinfo,'osolve%is_plastic','main',size(osolve%is_plastic),'bool',-1) ; deallocate (osolve%is_plastic) call heap (threadinfo,'osolve%is_plastic','main',size(osolve%is_plastic),'bool',-1) ; deallocate (osolve%is_plastic)
call heap (threadinfo,'ov%wpreiso','main',size(ov%wpreiso),'dp',-1) ; deallocate (ov%wpreiso) !call heap (threadinfo,'ov%wpreiso','main',size(ov%wpreiso),'dp',-1) ; deallocate (ov%wpreiso)
end if end if
!-------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------
...@@ -1044,6 +1047,7 @@ do while (istep.le.params%nstep) ...@@ -1044,6 +1047,7 @@ do while (istep.le.params%nstep)
osolve%u(i)=ov%unode(i) osolve%u(i)=ov%unode(i)
osolve%v(i)=ov%vnode(i) osolve%v(i)=ov%vnode(i)
osolve%w(i)=ov%wnode(i) osolve%w(i)=ov%wnode(i)
osolve%wpreiso(i)=ov%wnodepreiso(i)
enddo enddo
if (params%calculate_temp) then if (params%calculate_temp) then
...@@ -1241,6 +1245,7 @@ do while (istep.le.params%nstep) ...@@ -1241,6 +1245,7 @@ do while (istep.le.params%nstep)
osolve%u=ov%unode osolve%u=ov%unode
osolve%v=ov%vnode osolve%v=ov%vnode
osolve%w=ov%wnode osolve%w=ov%wnode
osolve%wpreiso=ov%wnodepreiso
osolve%temp=ov%temp osolve%temp=ov%temp
print *,'osolve values transferred ov' print *,'osolve values transferred ov'
...@@ -1282,6 +1287,7 @@ do while (istep.le.params%nstep) ...@@ -1282,6 +1287,7 @@ do while (istep.le.params%nstep)
call heap (threadinfo,'osolve%u','main',size(osolve%u),'dp',-1) ; deallocate (osolve%u) call heap (threadinfo,'osolve%u','main',size(osolve%u),'dp',-1) ; deallocate (osolve%u)
call heap (threadinfo,'osolve%v','main',size(osolve%v),'dp',-1) ; deallocate (osolve%v) call heap (threadinfo,'osolve%v','main',size(osolve%v),'dp',-1) ; deallocate (osolve%v)
call heap (threadinfo,'osolve%w','main',size(osolve%w),'dp',-1) ; deallocate (osolve%w) call heap (threadinfo,'osolve%w','main',size(osolve%w),'dp',-1) ; deallocate (osolve%w)
call heap (threadinfo,'osolve%wpreiso','main',size(osolve%wpreiso),'dp',-1) ; deallocate (osolve%wpreiso)
call heap (threadinfo,'osolve%temp','main',size(osolve%temp),'dp',-1) ; deallocate (osolve%temp) call heap (threadinfo,'osolve%temp','main',size(osolve%temp),'dp',-1) ; deallocate (osolve%temp)
call heap (threadinfo,'osolve%crit','main',size(osolve%crit),'dp',-1) ; deallocate (osolve%crit) call heap (threadinfo,'osolve%crit','main',size(osolve%crit),'dp',-1) ; deallocate (osolve%crit)
call heap (threadinfo,'osolve%e2d','main',size(osolve%e2d),'dp',-1) ; deallocate (osolve%e2d) call heap (threadinfo,'osolve%e2d','main',size(osolve%e2d),'dp',-1) ; deallocate (osolve%e2d)
...@@ -1295,7 +1301,7 @@ do while (istep.le.params%nstep) ...@@ -1295,7 +1301,7 @@ do while (istep.le.params%nstep)
call heap (threadinfo,'vo%rtf','main',size(vo%rtf),'int',-1) ; deallocate (vo%rtf) call heap (threadinfo,'vo%rtf','main',size(vo%rtf),'int',-1) ; deallocate (vo%rtf)
call heap (threadinfo,'vo%ftr','main',size(vo%ftr),'int',-1) ; deallocate (vo%ftr) call heap (threadinfo,'vo%ftr','main',size(vo%ftr),'int',-1) ; deallocate (vo%ftr)
call heap (threadinfo,'vo%influid','main',size(vo%influid),'bool',-1) ; deallocate (vo%influid) call heap (threadinfo,'vo%influid','main',size(vo%influid),'bool',-1) ; deallocate (vo%influid)
call heap (threadinfo,'ov%wpreiso','main',size(ov%wpreiso),'dp',-1) ; deallocate (ov%wpreiso) !call heap (threadinfo,'ov%wpreiso','main',size(ov%wpreiso),'dp',-1) ; deallocate (ov%wpreiso)
call show_time (total,step,inc,1,'End of time step$') call show_time (total,step,inc,1,'End of time step$')
...@@ -1334,7 +1340,7 @@ call heap (threadinfo,'ov%z','main',size(ov%z),'dp',-1) ; deallocat ...@@ -1334,7 +1340,7 @@ call heap (threadinfo,'ov%z','main',size(ov%z),'dp',-1) ; deallocat
call heap (threadinfo,'ov%unode','main',size(ov%unode),'dp',-1) ; deallocate (ov%unode) call heap (threadinfo,'ov%unode','main',size(ov%unode),'dp',-1) ; deallocate (ov%unode)
call heap (threadinfo,'ov%vnode','main',size(ov%vnode),'dp',-1) ; deallocate (ov%vnode) call heap (threadinfo,'ov%vnode','main',size(ov%vnode),'dp',-1) ; deallocate (ov%vnode)
call heap (threadinfo,'ov%wnode','main',size(ov%wnode),'dp',-1) ; deallocate (ov%wnode) call heap (threadinfo,'ov%wnode','main',size(ov%wnode),'dp',-1) ; deallocate (ov%wnode)
call heap (threadinfo,'ov%wpreiso','main',size(ov%wpreiso),'dp',-1) ; deallocate (ov%wpreiso) call heap (threadinfo,'ov%wnodepreiso','main',size(ov%wnodepreiso),'dp',-1) ; deallocate (ov%wnodepreiso)
call heap (threadinfo,'ov%icon','main',size(ov%icon),'int',-1) ; deallocate (ov%icon) 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%temp','main',size(ov%temp),'dp',-1) ; deallocate (ov%temp)
call heap (threadinfo,'ov%temporary_nodal...','main',size(ov%temporary_nodal_pressure),'dp',-1) call heap (threadinfo,'ov%temporary_nodal...','main',size(ov%temporary_nodal_pressure),'dp',-1)
...@@ -1342,7 +1348,6 @@ call heap (threadinfo,'ov%temporary_nodal...','main',size(ov%temporary_nodal_pre ...@@ -1342,7 +1348,6 @@ call heap (threadinfo,'ov%temporary_nodal...','main',size(ov%temporary_nodal_pre
call heap (threadinfo,'ov%whole_leaf...','main',size(ov%whole_leaf_in_fluid),'bool',-1) call heap (threadinfo,'ov%whole_leaf...','main',size(ov%whole_leaf_in_fluid),'bool',-1)
deallocate (ov%whole_leaf_in_fluid) deallocate (ov%whole_leaf_in_fluid)
do is=1,params%ns do is=1,params%ns
call heap (threadinfo,'surface(is)%x','main',size(surface(is)%x),'dp',-1) ; deallocate (surface(is)%x) call heap (threadinfo,'surface(is)%x','main',size(surface(is)%x),'dp',-1) ; deallocate (surface(is)%x)
call heap (threadinfo,'surface(is)%y','main',size(surface(is)%y),'dp',-1) ; deallocate (surface(is)%y) call heap (threadinfo,'surface(is)%y','main',size(surface(is)%y),'dp',-1) ; deallocate (surface(is)%y)
...@@ -1355,8 +1360,6 @@ do is=1,params%ns ...@@ -1355,8 +1360,6 @@ do is=1,params%ns
call heap (threadinfo,'surface(is)%icon','main',size(surface(is)%icon),'int',-1) ; deallocate (surface(is)%icon) call heap (threadinfo,'surface(is)%icon','main',size(surface(is)%icon),'int',-1) ; deallocate (surface(is)%icon)
end do end do
deallocate (surface) deallocate (surface)
deallocate (surface0) deallocate (surface0)
deallocate (mat) deallocate (mat)
......
src/define_ov.f90
+ 4
− 1
View file @ 120d3a4b
...@@ -95,6 +95,7 @@ if (params%irestart.eq.0) then ...@@ -95,6 +95,7 @@ if (params%irestart.eq.0) then
allocate (ov%unode(ov%nnode),stat=threadinfo%err) ; call heap (threadinfo,'ov%unode','define_ov',size(ov%unode),'dp',+1) allocate (ov%unode(ov%nnode),stat=threadinfo%err) ; call heap (threadinfo,'ov%unode','define_ov',size(ov%unode),'dp',+1)
allocate (ov%vnode(ov%nnode),stat=threadinfo%err) ; call heap (threadinfo,'ov%vnode','define_ov',size(ov%vnode),'dp',+1) allocate (ov%vnode(ov%nnode),stat=threadinfo%err) ; call heap (threadinfo,'ov%vnode','define_ov',size(ov%vnode),'dp',+1)
allocate (ov%wnode(ov%nnode),stat=threadinfo%err) ; call heap (threadinfo,'ov%wnode','define_ov',size(ov%wnode),'dp',+1) allocate (ov%wnode(ov%nnode),stat=threadinfo%err) ; call heap (threadinfo,'ov%wnode','define_ov',size(ov%wnode),'dp',+1)
allocate (ov%wnodepreiso(ov%nnode),stat=threadinfo%err) ; call heap (threadinfo,'ov%wnodepreiso','define_ov',size(ov%wnodepreiso),'dp',+1)
allocate (ov%temp(ov%nnode),stat=threadinfo%err) ; call heap (threadinfo,'ov%temp','define_ov',size(ov%temp),'dp',+1) allocate (ov%temp(ov%nnode),stat=threadinfo%err) ; call heap (threadinfo,'ov%temp','define_ov',size(ov%temp),'dp',+1)
! Line below added by dwhipp - 12/09 ! Line below added by dwhipp - 12/09
allocate (ov%pressure(ov%nleaves),stat=threadinfo%err) ; call heap (threadinfo,'ov%pressure','define_ov',size(ov%pressure),'dp',+1) allocate (ov%pressure(ov%nleaves),stat=threadinfo%err) ; call heap (threadinfo,'ov%pressure','define_ov',size(ov%pressure),'dp',+1)
...@@ -107,6 +108,7 @@ if (params%irestart.eq.0) then ...@@ -107,6 +108,7 @@ if (params%irestart.eq.0) then
ov%unode=0.d0 ov%unode=0.d0
ov%vnode=0.d0 ov%vnode=0.d0
ov%wnode=0.d0 ov%wnode=0.d0
ov%wnodepreiso=0.d0
! Line below uncommented by dwhipp - 12/09 ! Line below uncommented by dwhipp - 12/09
ov%pressure=0.d0 ov%pressure=0.d0
! Line below added by dwhipp - 12/09 ! Line below added by dwhipp - 12/09
...@@ -132,6 +134,7 @@ else ...@@ -132,6 +134,7 @@ else
allocate (ov%unode(ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%unode in define_ov$') allocate (ov%unode(ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%unode in define_ov$')
allocate (ov%vnode(ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%vnode in define_ov$') allocate (ov%vnode(ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%vnode in define_ov$')
allocate (ov%wnode(ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%wnode in define_ov$') allocate (ov%wnode(ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%wnode in define_ov$')
allocate (ov%wnodepreiso(ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%wnodepreiso in define_ov$')
allocate (ov%temp(ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%temp in define_ov$') allocate (ov%temp(ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%temp in define_ov$')
! Line below uncommented by dwhipp - 12/09 ! Line below uncommented by dwhipp - 12/09
allocate (ov%pressure(ov%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%pressure in define_ov$') allocate (ov%pressure(ov%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%pressure in define_ov$')
...@@ -147,7 +150,7 @@ else ...@@ -147,7 +150,7 @@ else
ov%unode(i), & ov%unode(i), &
ov%vnode(i), & ov%vnode(i), &
ov%wnode(i), & ov%wnode(i), &
wpreiso, & ov%wnodepreiso(i), &
(xlsf,j=1,nlsf), & (xlsf,j=1,nlsf), &
ov%temp(i), & ov%temp(i), &
ov%temporary_nodal_pressure(i), & ov%temporary_nodal_pressure(i), &
......
src/interpolate_ov_on_osolve.f90
+ 28
− 8
View file @ 120d3a4b
...@@ -50,7 +50,7 @@ integer iter ...@@ -50,7 +50,7 @@ integer iter
integer err,ie_ov,ie_level,ie_loc,nproc,ierr,iproc,k,ie,ie_osolve,i,j integer err,ie_ov,ie_level,ie_loc,nproc,ierr,iproc,k,ie,ie_osolve,i,j
double precision :: x,y,z,x0_leaf,y0_leaf,z0_leaf,dxyz_leaf double precision :: x,y,z,x0_leaf,y0_leaf,z0_leaf,dxyz_leaf
double precision,dimension(:),allocatable::pressure,pressurep double precision,dimension(:),allocatable::pressure,pressurep
double precision,dimension(:),allocatable::u,v,w,temp double precision,dimension(:),allocatable::u,v,w,wpreiso,temp
double precision,dimension(:),allocatable::countnode double precision,dimension(:),allocatable::countnode
double precision,dimension(:),allocatable::ov_nodepressure double precision,dimension(:),allocatable::ov_nodepressure
character*72 :: shift character*72 :: shift
...@@ -96,24 +96,37 @@ allocate (pressurep(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Erro ...@@ -96,24 +96,37 @@ allocate (pressurep(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Erro
allocate (u(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc u in main$') allocate (u(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc u in main$')
allocate (v(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc v in main$') allocate (v(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc v in main$')
allocate (w(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc w in main$') allocate (w(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc w in main$')
allocate (wpreiso(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc wpreiso in main$')
allocate (temp(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc temp in main$') allocate (temp(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc temp in main$')
u=0.d0 u=0.d0
v=0.d0 v=0.d0
w=0.d0 w=0.d0
wpreiso=0.d0
temp=0.d0 temp=0.d0
pressure=0.d0 pressure=0.d0
! interpolate velocity solution onto solve octree ! interpolate velocity solution onto solve octree
do k=1+iproc,osolve%nnode,nproc do k=1+iproc,osolve%nnode,nproc
call octree_interpolate_five & ! call octree_interpolate_five &
(5,ov%octree,ov%noctree,ov%icon, & ! (5,ov%octree,ov%noctree,ov%icon, &
! ov%nleaves,ov%nnode, &
! osolve%x(k),osolve%y(k),osolve%z(k), &
! ov%unode,u(k), &
! ov%vnode,v(k), &
! ov%wnode,w(k), &
! ov%temp,temp(k), &
! ov%temporary_nodal_pressure,pressure(k))
call octree_interpolate_six &
(6,ov%octree,ov%noctree,ov%icon, &
ov%nleaves,ov%nnode, & ov%nleaves,ov%nnode, &
osolve%x(k),osolve%y(k),osolve%z(k), & osolve%x(k),osolve%y(k),osolve%z(k), &
ov%unode,u(k), & ov%unode,u(k), &
ov%vnode,v(k), & ov%vnode,v(k), &
ov%wnode,w(k), & ov%wnode,w(k), &
ov%wnodepreiso,wpreiso(k), &
ov%temp,temp(k), & ov%temp,temp(k), &
ov%temporary_nodal_pressure,pressure(k)) ov%temporary_nodal_pressure,pressure(k))
enddo enddo
...@@ -134,6 +147,7 @@ endif ...@@ -134,6 +147,7 @@ endif
osolve%u=0.d0 osolve%u=0.d0
osolve%v=0.d0 osolve%v=0.d0
osolve%w=0.d0 osolve%w=0.d0
osolve%wpreiso=0.d0
osolve%temp=0.d0 osolve%temp=0.d0
pressurep=0.d0 pressurep=0.d0
...@@ -145,11 +159,12 @@ pressurep=0.d0 ...@@ -145,11 +159,12 @@ pressurep=0.d0
! if (ieee_is_nan(pressure(i))) print *,'field -> NaN in pressure$' ! if (ieee_is_nan(pressure(i))) print *,'field -> NaN in pressure$'
!end do !end do
call mpi_allreduce (u, osolve%u, osolve%nnode, mpi_double_precision,mpi_sum,mpi_comm_world,ierr) call mpi_allreduce (u, osolve%u, osolve%nnode, mpi_double_precision,mpi_sum,mpi_comm_world,ierr)
call mpi_allreduce (v, osolve%v, osolve%nnode, mpi_double_precision,mpi_sum,mpi_comm_world,ierr) call mpi_allreduce (v, osolve%v, osolve%nnode, mpi_double_precision,mpi_sum,mpi_comm_world,ierr)
call mpi_allreduce (w, osolve%w, osolve%nnode, mpi_double_precision,mpi_sum,mpi_comm_world,ierr) call mpi_allreduce (w, osolve%w, osolve%nnode, mpi_double_precision,mpi_sum,mpi_comm_world,ierr)
call mpi_allreduce (temp, osolve%temp, osolve%nnode, mpi_double_precision,mpi_sum,mpi_comm_world,ierr) call mpi_allreduce (wpreiso, osolve%wpreiso, osolve%nnode, mpi_double_precision,mpi_sum,mpi_comm_world,ierr)
call mpi_allreduce (pressure, pressurep, osolve%nnode, mpi_double_precision,mpi_sum,mpi_comm_world,ierr) call mpi_allreduce (temp, osolve%temp, osolve%nnode, mpi_double_precision,mpi_sum,mpi_comm_world,ierr)
call mpi_allreduce (pressure, pressurep, osolve%nnode, mpi_double_precision,mpi_sum,mpi_comm_world,ierr)
!if (iproc.eq.0) print *,'int.ov.on.osolve: min/max pressurep',minval(pressurep),maxval(pressurep) !if (iproc.eq.0) print *,'int.ov.on.osolve: min/max pressurep',minval(pressurep),maxval(pressurep)
...@@ -171,6 +186,7 @@ deallocate (pressurep) ...@@ -171,6 +186,7 @@ deallocate (pressurep)
deallocate (u) deallocate (u)
deallocate (v) deallocate (v)
deallocate (w) deallocate (w)
deallocate (wpreiso)
deallocate (temp) deallocate (temp)
! prepare the velo octree to receive the solution and carry ! prepare the velo octree to receive the solution and carry
...@@ -183,6 +199,7 @@ call heap (threadinfo,'ov%z','interpolate_ov...',size(ov%z),'dp',-1) ...@@ -183,6 +199,7 @@ call heap (threadinfo,'ov%z','interpolate_ov...',size(ov%z),'dp',-1)
call heap (threadinfo,'ov%unode','interpolate_ov...',size(ov%unode),'dp',-1) ; deallocate (ov%unode) call heap (threadinfo,'ov%unode','interpolate_ov...',size(ov%unode),'dp',-1) ; deallocate (ov%unode)
call heap (threadinfo,'ov%vnode','interpolate_ov...',size(ov%vnode),'dp',-1) ; deallocate (ov%vnode) call heap (threadinfo,'ov%vnode','interpolate_ov...',size(ov%vnode),'dp',-1) ; deallocate (ov%vnode)
call heap (threadinfo,'ov%wnode','interpolate_ov...',size(ov%wnode),'dp',-1) ; deallocate (ov%wnode) call heap (threadinfo,'ov%wnode','interpolate_ov...',size(ov%wnode),'dp',-1) ; deallocate (ov%wnode)
call heap (threadinfo,'ov%wnodepreiso','interpolate_ov...',size(ov%wnodepreiso),'dp',-1) ; deallocate (ov%wnodepreiso)
call heap (threadinfo,'ov%temp','interpolate_ov...',size(ov%temp),'dp',-1) ; deallocate (ov%temp) call heap (threadinfo,'ov%temp','interpolate_ov...',size(ov%temp),'dp',-1) ; deallocate (ov%temp)
call heap (threadinfo,'ov%icon','interpolate_ov...',size(ov%icon),'int',-1) ; deallocate (ov%icon) call heap (threadinfo,'ov%icon','interpolate_ov...',size(ov%icon),'int',-1) ; deallocate (ov%icon)
call heap (threadinfo,'ov%temporary...','interpolate_ov...',size(ov%temporary_nodal_pressure),'dp',-1) call heap (threadinfo,'ov%temporary...','interpolate_ov...',size(ov%temporary_nodal_pressure),'dp',-1)
...@@ -207,6 +224,8 @@ ov%nnode=osolve%nnode ...@@ -207,6 +224,8 @@ ov%nnode=osolve%nnode
call heap (threadinfo,'ov%vnode','interpolate_ov...',size(ov%vnode),'dp',+1) call heap (threadinfo,'ov%vnode','interpolate_ov...',size(ov%vnode),'dp',+1)
allocate (ov%wnode(ov%nnode),stat=threadinfo%err) allocate (ov%wnode(ov%nnode),stat=threadinfo%err)
call heap (threadinfo,'ov%wnode','interpolate_ov...',size(ov%wnode),'dp',+1) call heap (threadinfo,'ov%wnode','interpolate_ov...',size(ov%wnode),'dp',+1)
allocate (ov%wnodepreiso(ov%nnodepreiso),stat=threadinfo%err)
call heap (threadinfo,'ov%wnodepreiso','interpolate_ov...',size(ov%wnodepreiso),'dp',+1)
allocate (ov%temp(ov%nnode),stat=threadinfo%err) allocate (ov%temp(ov%nnode),stat=threadinfo%err)
call heap (threadinfo,'ov%temp','interpolate_ov...',size(ov%temp),'dp',+1) call heap (threadinfo,'ov%temp','interpolate_ov...',size(ov%temp),'dp',+1)
allocate (ov%temporary_nodal_pressure(ov%nnode),stat=threadinfo%err) allocate (ov%temporary_nodal_pressure(ov%nnode),stat=threadinfo%err)
...@@ -243,6 +262,7 @@ ov%temp=osolve%temp ...@@ -243,6 +262,7 @@ ov%temp=osolve%temp
ov%unode=0.d0 ov%unode=0.d0
ov%vnode=0.d0 ov%vnode=0.d0
ov%wnode=0.d0 ov%wnode=0.d0
ov%wnodepreiso=0.d0
!ov%pressure=0.d0 !ov%pressure=0.d0
ov%temporary_nodal_pressure=0.d0 ov%temporary_nodal_pressure=0.d0
......
src/isostasy.f90
+ 2
− 2
View file @ 120d3a4b
...@@ -263,13 +263,13 @@ enddo ...@@ -263,13 +263,13 @@ enddo
endif endif
! Store the total displacement of the basal node plane ! Store the total displacement of the basal node plane
if (params%isobc .and. iter.eq.1) zi%zisodisp=zi%zisodisp-work if (params%isobc .and. iter.eq.1) zi%zisodisp=zi%zisodisp-work
! transforms the displacements into velocities and add this contribution to the velocity ! 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 ! solution of the Stokes equation before interpolation onto the surfaces and their advection
do i=1,ov%nnode do i=1,ov%nnode
ov%wpreiso(i)=ov%wnode(i) ov%wnodepreiso(i)=ov%wnode(i)
ov%wnode(i)=ov%wnode(i)-diso(i)/params%dt ov%wnode(i)=ov%wnode(i)-diso(i)/params%dt
enddo enddo
......
src/module_definitions.f90
+ 2
− 3
View file @ 120d3a4b
...@@ -53,8 +53,7 @@ module definitions ...@@ -53,8 +53,7 @@ module definitions
integer,dimension(:),pointer::octree integer,dimension(:),pointer::octree
integer noctree,nnode,nleaves integer noctree,nnode,nleaves
double precision,dimension(:),pointer::x,y,z,temp double precision,dimension(:),pointer::x,y,z,temp
double precision,dimension(:),pointer::unode,vnode,wnode double precision,dimension(:),pointer::unode,vnode,wnode,wnodepreiso
double precision,dimension(:),pointer::wpreiso
double precision,dimension(:),pointer::pressure,spressure double precision,dimension(:),pointer::pressure,spressure
double precision,dimension(:),pointer::temporary_nodal_pressure double precision,dimension(:),pointer::temporary_nodal_pressure
integer,dimension(:,:),pointer::icon integer,dimension(:,:),pointer::icon
...@@ -94,7 +93,7 @@ module definitions ...@@ -94,7 +93,7 @@ module definitions
double precision,dimension(:),pointer::e2d double precision,dimension(:),pointer::e2d
double precision,dimension(:),pointer::e3d double precision,dimension(:),pointer::e3d
double precision,dimension(:),pointer::lode double precision,dimension(:),pointer::lode
double precision,dimension(:),pointer::u,v,w,temp double precision,dimension(:),pointer::u,v,w,wpreiso,temp
double precision,dimension(:,:),pointer::lsf double precision,dimension(:,:),pointer::lsf
double precision,dimension(:),pointer::eviscosity,q double precision,dimension(:),pointer::eviscosity,q
integer,dimension(:,:),pointer::icon,iface integer,dimension(:,:),pointer::icon,iface
......
src/write_global_output.f90
+ 2
− 2
View file @ 120d3a4b
...@@ -98,11 +98,11 @@ if (iproc.eq.0) then ...@@ -98,11 +98,11 @@ if (iproc.eq.0) then
! write info on octree solve nodes (x,y,z,u,v,w,lsf,temp) ! write info on octree solve nodes (x,y,z,u,v,w,lsf,temp)
write (9) (osolve%x(i), & write (9) (osolve%x(i), &
osolve%y(i), & osolve%y(i), &
osolve%z(i)*params%vex, & osolve%z(i)*params%vex, &
osolve%u(i), & osolve%u(i), &
osolve%v(i), & osolve%v(i), &
osolve%w(i), & osolve%w(i), &
ov%wpreiso(i), & osolve%wpreiso(i), &
(osolve%lsf(i,j),j=1,osolve%nlsf), & (osolve%lsf(i,j),j=1,osolve%nlsf), &
osolve%temp(i), & osolve%temp(i), &
ov%temporary_nodal_pressure(i), & ov%temporary_nodal_pressure(i), &
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment