added storage/output of pre-isostasy velocity field; output of viscosity,...

added storage/output of pre-isostasy velocity field; output of viscosity, is_plastic, elemental pressure

src/DOUAR.f90

@@ -176,7 +176,9 @@ if (params%irestart==0) then
 else
    current_level=params%levelmax_oct
 end if
- 
+
+if (params%isobc) allocate(zisodisp(2**params%levelmax_oct+1,2**params%levelmax_oct+1),stat=err) ; if (err.ne.0) call stop_run ('Error alloc zisodisp in main$')
+
 call show_time (total,step,inc,1,'Problem Setup$')
 
 !------------------------------------------------------------------------------|
@@ -782,7 +784,7 @@ 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%eviscosity','main',size(osolve%eviscosity),'dp',-1) ; deallocate (osolve%eviscosity)
       call heap (threadinfo,'osolve%q','main',size(osolve%q),'dp',-1)                   ; deallocate (osolve%q)                                                    
       !--------------------------------------------------------------------------------------
       !--------------------------------------------------------------------------------------
@@ -791,7 +793,8 @@ do while (istep.le.params%nstep)
       !--------------------------------------------------------------------------------------
 
       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$')
       if (params%isostasy) call isostasy (params,weightel,ov,surface,mat,zisodisp)
 
       !--------------------------------------------------------------------------------------
@@ -853,8 +856,10 @@ do while (istep.le.params%nstep)
          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%e3d','main',size(osolve%e3d),'dp',-1)                 ; deallocate (osolve%e3d)  
-         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%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)
       end if
 
       !--------------------------------------------------------------------------------------
@@ -1269,6 +1274,7 @@ do while (istep.le.params%nstep)
    call heap (threadinfo,'osolve%e2d','main',size(osolve%e2d),'dp',-1)                  ; deallocate (osolve%e2d)
    call heap (threadinfo,'osolve%e3d','main',size(osolve%e3d),'dp',-1)                  ; deallocate (osolve%e3d)
    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%is_plastic','main',size(osolve%is_plastic),'bool',-1)  ; deallocate (osolve%is_plastic)
    call heap (threadinfo,'vo%node','main',size(vo%node),'int',-1)                       ; deallocate (vo%node)
    call heap (threadinfo,'vo%leaf','main',size(vo%leaf),'int',-1)                       ; deallocate (vo%leaf)
@@ -1276,6 +1282,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%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,'ov%wpreiso','main',size(ov%wpreiso),'dp',-1)                  ; deallocate (ov%wpreiso)
 
    call show_time (total,step,inc,1,'End of time step$')