Revised isostasy B/C code logic to do final calculation outside of linear iterations

src/DOUAR.f90

@@ -794,11 +794,12 @@ do while (istep.le.params%nstep)
       !     compute isostasy
       !--------------------------------------------------------------------------------------
       !--------------------------------------------------------------------------------------
-
       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$')
-        call isostasy (params,weightel,ov,surface,mat,iter,zi)
+        if (.not.converge) 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$')
+          call isostasy (params,weightel,ov,surface,mat,iter,zi)
+        endif
       endif
 
       !--------------------------------------------------------------------------------------
@@ -1047,7 +1048,7 @@ do while (istep.le.params%nstep)
       osolve%u(i)=ov%unode(i)
       osolve%v(i)=ov%vnode(i)
       osolve%w(i)=ov%wnode(i)
-      osolve%wpreiso(i)=ov%wnodepreiso(i)
+      !osolve%wpreiso(i)=ov%wnodepreiso(i)
    enddo
 
    if (params%calculate_temp) then
@@ -1209,15 +1210,21 @@ do while (istep.le.params%nstep)
 
    enddo
 
-   call interpolate_velocity_on_surface(params,surface,ov)
+   !call interpolate_velocity_on_surface(params,surface,ov)
 
    !--------------------------------------------------------------------------------------
    !--------------------------------------------------------------------------------------
-   !     update Moho displacement
+   !     apply isostasy and update Moho displacement
    !--------------------------------------------------------------------------------------
    !--------------------------------------------------------------------------------------
-   call show_time (total,step,inc,1,'Update Moho displacement$')
-   zi%zisodisp=zi%zisodisp-zi%zisoinc
+   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$')
+     call isostasy (params,weightel,ov,surface,mat,iter,zi)
+     zi%zisodisp=zi%zisodisp-zi%zisoinc
+   endif
+
+   call interpolate_velocity_on_surface(params,surface,ov)
 
    !--------------------------------------------------------------------------------------
    !--------------------------------------------------------------------------------------