added pre-isostasy velocity storage/output

src/DOUAR.f90

@@ -1314,6 +1314,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%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%wpreiso','main',size(ov%wpreiso),'dp',-1)  ; deallocate (ov%wpreiso)
 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) 

src/define_ov.f90

@@ -94,6 +94,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%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%wpreiso(ov%nnode),stat=threadinfo%err)  ;  call heap (threadinfo,'ov%wpreiso','define_ov',size(ov%wnode),'dp',+1)
    allocate (ov%temp(ov%nnode),stat=threadinfo%err)     ;  call heap (threadinfo,'ov%temp','define_ov',size(ov%temp),'dp',+1)
    allocate (ov%temporary_nodal_pressure(ov%nnode),stat=threadinfo%err)
    call heap (threadinfo,'ov%temporary_nodal_pressure','define_ov',size(ov%temporary_nodal_pressure),'dp',+1)
@@ -102,6 +103,7 @@ if (params%irestart.eq.0) then
    ov%unode=0.d0
    ov%vnode=0.d0
    ov%wnode=0.d0
+   ov%wpreiso=0.d0
 !   ov%pressure=0.d0
    ov%temporary_nodal_pressure=0.d0
    call initialize_temperature (ov,params%ztemp)
@@ -121,9 +123,10 @@ else
    allocate (ov%y(ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%y in define_ov$')
    allocate (ov%z(ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%z in define_ov$')
    allocate (ov%icon(8,ov%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%icon 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%wnode(ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%wnode 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%wnode(ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%wnode in define_ov$')
+   allocate (ov%wpreiso(ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%wpreiso 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%pressure(ov%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%pressure in define_ov$')
    allocate (ov%temporary_nodal_pressure(ov%nnode),stat=err)
@@ -136,6 +139,7 @@ else
                ov%unode(i),     &
                ov%vnode(i),     &
                ov%wnode(i),     &
+               ov%wpreiso(i),     &
                (xlsf,j=1,nlsf), &
                ov%temp(i),      &
                ov%temporary_nodal_pressure(i), &

src/isostasy.f90

@@ -263,6 +263,7 @@ np=0
 ! solution of the Stokes equation before interpolation onto the surfaces and their advection
 
   do i=1,ov%nnode
+  ov%wpreiso(i)=ov%wnode(i)
   ov%wnode(i)=ov%wnode(i)-diso(i)/params%dt
   enddo
 

src/module_definitions.f90

@@ -45,6 +45,7 @@ module definitions
 ! nleaves is number of leaves in the octree
 ! x,y,z are the coordinates of the nodes
 ! unode,vnode,wnode are the velocity conponents at the nodes
+! wpreiso is the velocity z-component prior to isostasy
 ! temp is temperature at the nodes
 ! icon is connectivity array between nodes
 
@@ -53,6 +54,7 @@ module definitions
       integer noctree,nnode,nleaves
       double precision,dimension(:),pointer::x,y,z,temp
       double precision,dimension(:),pointer::unode,vnode,wnode
+      double precision,dimension(:),pointer::wpreiso
       double precision,dimension(:),pointer::pressure
       double precision,dimension(:),pointer::temporary_nodal_pressure
       integer,dimension(:,:),pointer::icon

src/write_global_output.f90

@@ -100,6 +100,7 @@ if (iproc.eq.0) then
                 osolve%u(i),                         &
                 osolve%v(i),                         & 
                 osolve%w(i),                         &
+                ov%wpreiso(i),                       &
                 (osolve%lsf(i,j),j=1,osolve%nlsf),   &
                 osolve%temp(i),                      &
                 ov%temporary_nodal_pressure(i),      &