Skip to content
Snippets Groups Projects
read_restart.f90 11.8 KiB
Newer Older
  • Learn to ignore specific revisions
  • subroutine read_restart(params,istep,current_time,osolve,ov,vo,surface,cl,&
    
    
    use definitions
    !use mpi
    
    implicit none
    
    include 'mpif.h'
    
    ! Passed variables
    type (parameters) :: params
    
    double precision :: current_time
    type (octreesolve) :: osolve
    type (octreev) :: ov
    type (void) :: vo
    type (sheet) :: surface(params%ns)
    type (cloud) :: cl
    
    type (nest_info) :: nest
    
    ! Internal variables
    
    double precision :: junk
    integer :: err,i,ires,is,j,k,stat,zisosize
    character(len=1) :: answer
    
    ! Read required values from the input.txt file
    write (*,'(a,a)') 'Reading input file ',trim(params%infile)
    params%vex=1.d0
    call scanfile (params%infile,'vex',params%vex,ires)
    params%isobc=.false.
    call scanfile (params%infile,'isobc',answer,ires)  
    params%isobc=(trim(answer)=='T') 
    params%nest=.false.
    call scanfile (params%infile,'nest',answer,ires)  
    params%nest=(trim(answer)=='T') 
    
    ! Echo values to screen
    write (*,'(a,a,a)') 'Input values from ',trim(params%infile),':'
    write (*,'(a,e11.4)') 'vex =',params%vex
    write (*,'(a,l1)') 'isobc = ',params%isobc
    write (*,'(a,l1)') 'nest = ',params%nest
    
    ! Open/read restart file
    
    open (9,file=trim(params%restartfile),status='old',form='unformatted')
    
    write (*,'(a,a)') 'Reading restart file ',trim(params%restartfile)
    write (*,'(a)') ''
    
             osolve%nnode,     &
             osolve%nleaves,   &
             osolve%nface,     &
             osolve%nlsf,      &
             cl%np,            &
             current_time
    
    ! Massive allocation
    allocate (osolve%x(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%x in read_restart$')
    allocate (osolve%y(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%y in read_restart$')
    allocate (osolve%z(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%z in read_restart$')
    
    allocate (osolve%u(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%u in read_restart$')
    allocate (osolve%v(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%v in read_restart$')
    allocate (osolve%w(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%w in read_restart$')
    allocate (osolve%wiso(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%wiso in read_restart$')
    
    allocate (osolve%lsf(osolve%nnode,osolve%nlsf),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%lsf in read_restart$')
    allocate (osolve%temp(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%temp in read_restart$')
    allocate (osolve%strain(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%strain in read_restart$')
    allocate (osolve%kfix(osolve%nnode*3),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%kfix in read_restart$')
    allocate (osolve%kfixt(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%kfixt in read_restart$')
    allocate (osolve%icon(8,osolve%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%icon in read_restart$')
    allocate (osolve%pressure(osolve%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%pressure in read_restart$') 
    allocate (osolve%spressure(osolve%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%spressure in read_restart$')
    allocate (osolve%crit(osolve%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%crit in read_restart$')
    allocate (osolve%e2d(osolve%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%e2d in read_restart$')
    allocate (osolve%eviscosity(osolve%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%eviscosity in read_restart$')
    allocate (osolve%is_plastic(osolve%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%is_plastic in read_restart$')
    allocate (osolve%dilatr(osolve%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%dilatr in read_restart$')
    allocate (osolve%matnum(osolve%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%matnum in read_restart$')
    
    allocate (osolve%octree(osolve%noctree),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%octree in read_restart$')
    
    allocate (osolve%iface(9,osolve%nface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%iface in read_restart$')
    
    allocate (vo%node(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc vo%node in read_restart$')
    allocate (vo%leaf(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc vo%leaf in read_restart$')
    allocate (vo%ftr(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc vo%ftr in read_restart$')
    allocate (vo%rtf(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc vo%rtf in read_restart$')
    allocate (vo%influid(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc vo%influid in read_restart$')
    allocate (vo%face(osolve%nface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc vo%face in read_restart$')
    
    allocate (ov%temporary_nodal_pressure(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%temp_nodal_pressure in read_restart$')
    allocate (ov%whole_leaf_in_fluid(osolve%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%whole_leaf_in_fluid in read_restart$')
    
    allocate (cl%x(cl%np), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%x in read_restart$')
    allocate (cl%y(cl%np), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%y in read_restart$')
    allocate (cl%z(cl%np), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%z in read_restart$')
    allocate (cl%x0(cl%np), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%x0 in read_restart$')
    allocate (cl%y0(cl%np), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%x0 in read_restart$')
    allocate (cl%z0(cl%np), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%x0 in read_restart$')
    allocate (cl%strain(cl%np), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%strain in read_restart$')
    allocate (cl%lsf0(cl%np), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%lsf0 in read_restart$')
    allocate (cl%temp(cl%np), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%temp in read_restart$')
    allocate (cl%press(cl%np), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%press in read_restart$')
    allocate (cl%tag(cl%np), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%tag in read_restart$')
    
    
    ! Read info for octree solve nodes (x,y,z,u,v,w,lsf,temp)
    read (9) (osolve%x(i),                         &
              osolve%y(i),                         &  
              osolve%z(i),                         &
              osolve%u(i),                         &
              osolve%v(i),                         & 
              osolve%w(i),                         &
              osolve%wiso(i),                      &
    
              (osolve%lsf(i,j),j=1,osolve%nlsf),   &
    
              osolve%temp(i),                      &
              ov%temporary_nodal_pressure(i),      &
              osolve%strain(i),                    &
              osolve%kfix((i-1)*3+1),              &
              osolve%kfix((i-1)*3+2),              &
              osolve%kfix((i-1)*3+3),              &
              osolve%kfixt(i),                     &
              i=1,osolve%nnode)
    
    
    ! Correct for vertical exaggeration
    osolve%z=osolve%z/params%vex
    
    read (9) ((osolve%icon(j,i),j=1,8),   &
    
               osolve%pressure(i),        &
               osolve%spressure(i),       &
               osolve%crit(i),            &
               osolve%e2d(i),             &
               osolve%eviscosity(i),      & 
               osolve%is_plastic(i),      &
               osolve%dilatr(i),          &
               osolve%matnum(i),          &
    
               ov%whole_leaf_in_fluid(i), &
               i=1,osolve%nleaves)
    
    ! Read octree information
    
    read (9) (osolve%octree(i),i=1,osolve%noctree)
    
    
    ! Read bad face information
    
    read (9) ((osolve%iface(j,i),j=1,9),i=1,osolve%nface)
    
    
    ! Read void information
    read (9) (vo%node(i),    &
              vo%leaf(i),    &
              vo%ftr(i),     &
              vo%rtf(i),     &
              vo%influid(i), &
              i=1,osolve%nnode)
    
    ! Read bad faces information
    read (9) (vo%face(i),i=1,osolve%nface)
    
    ! Read surface information (r,s,x,y,z,xn,yn,zn)
    do is=1,osolve%nlsf
    
      read (9) surface(is)%nsurface,surface(is)%activation_time,surface(is)%nt
    
    
      allocate (surface(is)%r(surface(is)%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc surface(is)%r in read_restart$')
      allocate (surface(is)%s(surface(is)%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc surface(is)%s in read_restart$')
      allocate (surface(is)%x(surface(is)%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc surface(is)%x in read_restart$')
      allocate (surface(is)%y(surface(is)%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc surface(is)%y in read_restart$')
      allocate (surface(is)%z(surface(is)%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc surface(is)%z in read_restart$')
      allocate (surface(is)%xn(surface(is)%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc surface(is)%xn in read_restart$')
      allocate (surface(is)%yn(surface(is)%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc surface(is)%yn in read_restart$')
      allocate (surface(is)%zn(surface(is)%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc surface(is)%zn in read_restart$')
      allocate (surface(is)%u(surface(is)%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc surface(is)%u in read_restart$')
      allocate (surface(is)%v(surface(is)%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc surface(is)%v in read_restart$')
      allocate (surface(is)%w(surface(is)%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc surface(is)%w in read_restart$')
      allocate (surface(is)%icon(3,surface(is)%nt),stat=err) ; if (err.ne.0) call stop_run ('Error alloc surface(is)%icon in read_restart$')
    
      read (9) (surface(is)%r(i),  &
                surface(is)%s(i),  &
                surface(is)%x(i),  &
                surface(is)%y(i),  &
                surface(is)%z(i),  &
                surface(is)%xn(i), &
                surface(is)%yn(i), &
                surface(is)%zn(i), &
                surface(is)%u(i),  &
                surface(is)%v(i),  &
                surface(is)%w(i),  &
                i=1,surface(is)%nsurface)
    
      read (9) ((surface(is)%icon(j,i),j=1,3),i=1,surface(is)%nt)
    
    
      ! correct surface for vertical exageration
      surface(is)%z=surface(is)%z/params%vex
      surface(is)%zn=surface(is)%zn/params%vex
    enddo
    
    ! Read cloud information
    read (9) (cl%x(i),             &
              cl%y(i),             &
              cl%z(i),             &
              cl%x0(i),            &
              cl%y0(i),            &
              cl%z0(i),            &
              cl%strain(i),        &
              cl%lsf0(i),          &
              cl%temp(i),          &
              cl%press(i),         &
              cl%tag(i),           &
              i=1,cl%np)
    
    ! correct for vertical exaggeration
    cl%z=cl%z/params%vex
    cl%z0=cl%z0/params%vex
    
    ! Read isostasy basal displacement array
    
    allocate (bcdef%zisodisp(zisosize+1,zisosize+1), stat=err) ; if (err.ne.0) call stop_run ('Error alloc bcdef%zisodisp in read_restart$')
    read (9) ((bcdef%zisodisp(i,j),j=1,zisosize+1),i=1,zisosize+1)
    bcdef%zisodisp=bcdef%zisodisp-surface(osolve%nlsf)%sp01
    
    
    ! Read nested model info
    if (params%nest) then
      read (9) nest%sselemx,nest%sselemy,nest%sselemz,nest%xminls,nest%yminls,&
               nest%zminls
    else
      nest%sselemx=0.d0
      nest%sselemy=0.d0
      nest%sselemz=0.d0
      nest%xminls=0.d0
      nest%yminls=0.d0
      nest%zminls=0.d0
    endif
    
    
    !write (*,*) '******************************************************************'
    !write (*,*) 'Reading to EOF'
    !write (*,*) '******************************************************************'
    !stat=0
    !do i=1,1000000
    !  read (9,iostat=stat) junk
    !  if (stat /= 0) exit
    !  write (*,*) junk
    !enddo