Subroutine rwdata(irw)
Use comm

Integer, Intent(In) :: irw
Real(Kind=8), Pointer :: bigstor(:,:,:)
Real(Kind=8), Allocatable :: ingrids(:,:)
Integer :: j,k
Integer(Kind=8) :: glength
Character(Len=256) :: oname

Select Case(irw)
Case(0)
  Allocate(ingrids(mmax+2,4),Stat=error)
  If (error/=0) Then
    Print *,"RWDATA: Error allocating memory to ingrids ", error
    Stop
  End If

  Inquire(IOLENGTH=glength) ingrids
  oname = Trim(AdjustL(grid_name))//".bin"
  Print *,' <Reading ',Trim(AdjustL(oname))
  Open (Unit=30,File=Trim(AdjustL(oname)), form="unformatted", access="direct", recl=glength, Status="Old",IOStat=error)
  If (error/=0) Then
    Print *, "RWDATA: Error opening ",Trim(AdjustL(oname)), error
    Stop
  End If
  Read(30,rec=1) ingrids
  Close(30,IOStat=error)
  If (error/=0) Print *, "RWDATA: Error closing ",Trim(AdjustL(oname)), error

    y = ingrids(1:nyo,1)
    z = ingrids(1:nzo,2)
    dy = ingrids(1:nyo,3)
    dz = ingrids(1:nzo,4)

  Inquire(IOLENGTH=glength) ian
  oname = Trim(AdjustL(grid_name2))//".bin"
  Print *,' <Reading ',Trim(AdjustL(oname))
  Open (Unit=30,File=Trim(AdjustL(oname)), form="unformatted", access="direct", recl=glength, Status="Old",IOStat=error)
  If (error/=0) Then
    Print *, "RWDATA: Error opening",Trim(AdjustL(oname)), error
    Stop
  End If
  Read(30,rec=1) ian
  Close(30,IOStat=error)
  If (error/=0) Print *, "RWDATA: Error closing",Trim(AdjustL(oname)), error

Case(1)
  Allocate(bigstor(0:nyo+1,0:nzo+1,7),Stat=error)
  If (error/=0) Then
    Print *,"RWDATA: Error allocating memory to bigstor", error
    Stop
  End If

  If (nyo>ny.OR.nzo>nz) Then
    Write(*,*) "rwdata: nyo or nzo too big", ny, nz, nyo, nzo
    Stop
  End If

  oname = Trim(AdjustL(output_name))//".bin"
  Print *,' <Reading ',Trim(AdjustL(oname))
  Open (Unit=30,File=Trim(AdjustL(oname)), form="unformatted", access="direct", recl=bigstor_length, Status="Old",IOStat=error)
  If (error/=0) Then
    Print *, "RWDATA: Error opening",Trim(AdjustL(oname)), error
    Stop
  End If
  Read(30,rec=1) bigstor
  Close(30,IOStat=error)
  If (error/=0) Print *, "RWDATA: Error closing",Trim(AdjustL(oname)), error

  y(0:nyo+1) = bigstor(0:nyo+1,1,1)
  z(0:nzo+1) = bigstor(1,0:nzo+1,2)
  ForAll(j=0:nyo+1,k=0:nzo+1)
    phi(j,k) = bigstor(j,k,3) 
    rne(j,k) = bigstor(j,k,4) 
    rnp(j,k) = bigstor(j,k,5) 
    rnn(j,k) = bigstor(j,k,6) 
    enth(j,k) = bigstor(j,k,7) 
  End ForAll

  DeAllocate(bigstor,Stat=error)
  If (error/=0) Print *,"RWDATA: Error deallocating memory to bigstor", error

Case(2)
  oname = Trim(AdjustL(output_name))//"_time.bin"
  Inquire(IOLENGTH=bigstor_length) timekeep
  Open(Unit=50,File=Trim(AdjustL(oname)), form="unformatted", access="direct", recl=bigstor_length, Status="replace",IOStat=error)
  If (error/=0) Then
    Print *, "RWDATA: Error opening ", Trim(AdjustL(oname)), error, bigstor_length
  Else
    Write(50,rec=1) timekeep
    Close(50,IOStat=error)
    If (error/=0) Print *, "RWDATA: Error closing ", Trim(AdjustL(oname)), error
  End If
  
  Allocate(bigstor(0:ny+1,0:nz+1,7),Stat=error)
  If (error/=0) Then
    Print *,"RWDATA: Error allocating memory to bigstor", error
    Stop
  End If
  
  ForAll(j=0:ny+1,k=0:nz+1)
    bigstor(j,k,1) = y(j)
    bigstor(j,k,2) = z(k)
    bigstor(j,k,3) = phi(j,k)
    bigstor(j,k,4) = rne(j,k)
    bigstor(j,k,5) = rnp(j,k)
    bigstor(j,k,6) = rnn(j,k)
    bigstor(j,k,7) = enth(j,k)
  End ForAll
  
  oname = Trim(AdjustL(output_name))//".bin"
  Inquire(IOLENGTH=bigstor_length) bigstor
  Open(Unit=10,File=Trim(AdjustL(oname)), form="unformatted", access="direct", recl=bigstor_length, Status="replace",IOStat=error)
  If (error/=0) Then
    Print *, "RWDATA: Error opening ", Trim(AdjustL(oname)), error, bigstor_length
  Else
    Write(10,rec=1) bigstor
    Close(10,IOStat=error)
    If (error/=0) Print *, "RWDATA: Error closing ", Trim(AdjustL(oname)), error
  End If

  nyo=ny
  nzo=nz

  Open(11, File=Trim(AdjustL(input_file)), Action = "Write", Status = "replace", IOStat=error) 
  input_file = Trim(AdjustL(input_file))
  ifread = .True.
  If (error/=0) Then
    Print *, "RWDATA: Error opening new ", Trim(AdjustL(input_file)), error
    Stop
  End If
  Write(11, NML=inputs)                                                              
  Close(11, IOStat=error)                                                
  If (error/=0) Print *, "RWDATA: Error closing new", Trim(AdjustL(input_file)), error

Case Default
  Print *, "RWDAYA: Error, irw must have value between 0 and 3 ", irw
End Select

End Subroutine rwdata
