Updated isostasy code to update Moho position at the end of each time step

src/DOUAR.f90

@@ -1211,6 +1211,14 @@ do while (istep.le.params%nstep)
 
    call interpolate_velocity_on_surface(params,surface,ov)
 
+   !--------------------------------------------------------------------------------------
+   !--------------------------------------------------------------------------------------
+   !     update Moho displacement
+   !--------------------------------------------------------------------------------------
+   !--------------------------------------------------------------------------------------
+   call show_time (total,step,inc,1,'Update Moho displacement$')
+   zi%zisodisp=zi%zisodisp-zi%zisoinc
+
    !--------------------------------------------------------------------------------------
    !--------------------------------------------------------------------------------------
    !     update strain, pressure, temperature  on 3D cloud
@@ -1227,12 +1235,8 @@ do while (istep.le.params%nstep)
    call show_time (total,step,inc,1,'Advect the cloud$')
 !   call move_cloud (params,cl,ov%octree,ov%noctree,ov%unode,ov%vnode,ov%wnode,ov%nnode,ov%icon,ov%nleaves)
 
-   print *,'Cloud will be advected'
-
    call move_cloud (params,cl,ov)
 
-   print *,'Cloud has been moved'
-
    call heap (threadinfo,'olsf%x','main',size(olsf%x),'dp',-1)            ; deallocate (olsf%x) 
    call heap (threadinfo,'olsf%y','main',size(olsf%y),'dp',-1)            ; deallocate (olsf%y) 
    call heap (threadinfo,'olsf%z','main',size(olsf%z),'dp',-1)            ; deallocate (olsf%z)
@@ -1240,25 +1244,19 @@ do while (istep.le.params%nstep)
    call heap (threadinfo,'olsf%icon','main',size(olsf%icon),'int',-1)     ; deallocate (olsf%icon)
    call heap (threadinfo,'olsf%octree','main',size(olsf%octree),'int',-1) ; deallocate (olsf%octree)
 
-   print *,'olsf stuff deallocated'
-
    osolve%u=ov%unode
    osolve%v=ov%vnode
    osolve%w=ov%wnode
    osolve%wpreiso=ov%wnodepreiso
    osolve%temp=ov%temp
 
-   print *,'osolve values transferred ov'
-         
    !--------------------------------------------------------------------------------------
    !--------------------------------------------------------------------------------------
    !     write global output before deallocating arrays
    !--------------------------------------------------------------------------------------
    !--------------------------------------------------------------------------------------
    call show_time (total,step,inc,1,'Write output$')
-   print *,'Output will be generated now'
    call write_global_output (params,istep,iter-1,current_time,osolve,ov,vo,surface,cl,zi,'final')
-   print *,'Output written'
    call mpi_barrier (mpi_comm_world,ierr)
 
    !--------------------------------------------------------------------------------------
@@ -1364,7 +1362,10 @@ deallocate (surface)
 deallocate (surface0)
 deallocate (mat)
 !if (params%isobc) deallocate (zisodisp)
-if (params%isobc) deallocate (zi%zisodisp)
+if (params%isobc) then
+  deallocate (zi%zisodisp)
+  deallocate (zi%zisoinc)
+endif
 deallocate (params%materialn)
 
 call heap_final (threadinfo)

src/define_cloud.f90

@@ -73,7 +73,9 @@ if (params%irestart.eq.0) then
    allocate (cl%tag(1), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%tag in define_cloud$')
    if (params%isobc) then
      allocate (zi%zisodisp(2**params%levelmax_oct+1,2**params%levelmax_oct+1), stat=err); if (err.ne.0) call stop_run ('Error alloc zi%zisodisp in define_cloud$')
+     allocate (zi%zisoinc(2**params%levelmax_oct+1,2**params%levelmax_oct+1), stat=err); if (err.ne.0) call stop_run ('Error alloc zi%zisoinc in define_cloud$')
      zi%zisodisp=0.d0
+     zi%zisoinc=0.d0
    endif
 
 else
@@ -97,7 +99,10 @@ else
    allocate (cl%temp(cl%np),stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%temp in define_cloud$')
    allocate (cl%press(cl%np),stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%press in define_cloud$')
    allocate (cl%tag(cl%np),stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%tag in define_cloud$')
-   if (params%isobc) allocate (zi%zisodisp(2**params%levelmax_oct+1,2**params%levelmax_oct+1), stat=err); if (err.ne.0) call stop_run ('Error alloc zi%zisodisp in define_cloud$')
+   if (params%isobc) then
+     allocate (zi%zisodisp(2**params%levelmax_oct+1,2**params%levelmax_oct+1), stat=err); if (err.ne.0) call stop_run ('Error alloc zi%zisodisp in define_cloud$')
+     allocate (zi%zisoinc(2**params%levelmax_oct+1,2**params%levelmax_oct+1), stat=err); if (err.ne.0) call stop_run ('Error alloc zi%zisoinc in define_cloud$')
+   endif
 
    read (19) (x,               &
                 y,               &
@@ -155,6 +160,7 @@ else
      read (19) dumpi
      read (19) ((zi%zisodisp(i,j),j=1,2**params%levelmax_oct+1),&
                                   i=1,2**params%levelmax_oct+1)
+     zi%zisoinc=0.d0
    endif
 
 ! correct for vertical exaggeration

src/isostasy.f90

@@ -263,7 +263,8 @@ enddo
   endif
 
 ! 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
+  if (params%isobc .and. iter.eq.1) zi%zisoinc=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

src/module_definitions.f90

@@ -330,7 +330,7 @@ module definitions
 ! due to isostasy
 
       type ziso
-        double precision,dimension(:,:),pointer::zisodisp
+        double precision,dimension(:,:),pointer::zisodisp,zisoinc
       end type ziso
 
 end module definitions