!
!
! Subroutine Read_Grid
!
!
Subroutine Read_Grid
! 
! read into a mesh grid file. If file does not exist, then display a messagebox
! If read successfully, Gridloaded set to .True. Otherwise, it remains .false.
! Utilize control switch: GrayMenu
!
use dflib
use global

implicit none
character(len=132) input_record
character(len=512) msg0, msg1
integer(4) iret, ierr, ierr2
integer(4) i, j, no, n1, n2, n3, n4
real lon, lat, elev, q, hc,hm, dips1,dips2, offst, density_anomaly, cooling_curvature
real(8) Asum, Bsum, Csum

msg0 = ' 'C
msg1 = ' 'C
open(unit = 100, file = feg_inp, status = 'old', Action = 'read', iostat = ierr, PAD = "YES")
!NOTE: PAD = "YES" allows reading either old OrbData/Shells .feg, or newer OrbData5/Shells .feg
if (ierr /= 0) then
    msg0 = ' error opening input file ' //feg_inp//' 'C
	msg1 = ' error opening file 'C
	iret = messageboxqq(msg0, msg1, MB$ICONEXCLAMATION.OR.MB$OK)
	return
endif
read(100,'(A80)') title
! node
read(100, *, iostat = ierr) numnod, nrealn, nfaken, n1000, brief
if (ierr /= 0) then
    call error1
	return
endif
! allocate memory

1000 if (GrayMenu) then
  mxnode = numnod
else
  if (numnod > mxnode) then
     call error2('node',numnod, mxnode)
	 ! NUMNOD = 0
     ! return
     !
	 ! modified on March 10, 2005, ask for node, element, fault size again
	 !
	 call InitArraySize
	 ! after call, mxnode, mxel, mxfel change
	 if(numnod <= mxnode) then
	      if(allocated(nodeABG)) deallocate(nodeABG)
          if(allocated(eqcm)) deallocate(eqcm)
		  if(allocated(NMemo)) deallocate(NMemo)
     else 
	      goto 1000
     endif
  end if
end if
if (.not.allocated(nodeABG)) allocate(nodeABG(3,mxnode))
if (.not.allocated(eqcm)) allocate(eqcm(6,mxnode)) ! but (5:6,mxnod) may remain filled with zeros, unless OrbData5 = .TRUE.
if (.not.allocated(NMemo)) allocate(NMemo(mxnode))
! node input
Asum = 0.0D0
Bsum = 0.0D0
Csum = 0.0D0
OrbData5 = .FALSE. ! global switch initialized as FALSE, unless test below succeeds:
nodeinput: do i = 1, numnod
  READ (100, "(A)", IOSTAT = ierr) input_record
 !Note: This indirect reading method prevents reading TWO records in the search for density_anomaly!
  READ (input_record, *, IOSTAT = ierr) no, lon, lat, elev, q, hc, hm, density_anomaly, cooling_curvature
  IF (ierr /= 0) THEN ! perhaps the last 2 columns are not there?
    density_anomaly = 0.0
    cooling_curvature = 0.0
    READ (input_record, *, IOSTAT = ierr2) no, lon, lat, elev, q, hc, hm
    IF (ierr2 /= 0) THEN
        call error1
        return
    END IF
  endif
  lon = lon * .017453293
  lat = lat * .017453293
  nodeABG(1,i) = cos(lat)*cos(lon)
  nodeABG(2,i) = cos(lat)*sin(lon)
  nodeABG(3,i) = sin(lat)
  Asum = Asum + nodeABG(1,i)
  Bsum = Bsum + nodeABG(2,i)
  Csum = Csum + nodeABG(3,i)
  eqcm(1,i) = elev
  eqcm(2,i) = q
  eqcm(3,i) = hc
  eqcm(4,i) = hm
  eqcm(5,i) = density_anomaly   ! may be 0.0 in older files
  eqcm(6,i) = cooling_curvature ! may be 0.0 in older files
  OrbData5 = OrbData5 .OR. (density_anomaly /= 0.0) .OR. (cooling_curvature /= 0.0)
end do nodeinput
Asum = Asum/dble(numnod)
Bsum = Bsum/dble(numnod)
Csum = Csum/dble(numnod)
call unitize(Asum, Bsum, Csum)
nettempvec(1) = Asum
nettempvec(2) = Bsum
nettempvec(3) = Csum
! elevment
read(100, *, iostat = ierr) numel
if (ierr /= 0) then
    call error1
	return
endif
! allocate memory
2000 if (GrayMenu) then
  mxel = numel
else
  if (numel > mxel) then
     call error2('element',numel, mxel)
	 ! NUMEL = 0
     ! return
     !
	 ! modified on March 10, 2005, ask for node, element, fault size again
	 !
	 call InitArraySize
	 ! after call, mxnode, mxel, mxfel change
	 if(numel <= mxel) then
	      if(allocated(nodes)) deallocate(nodes)
          if(allocated(EMemo)) deallocate(EMemo)
     else 
	      goto 2000
     endif
  end if
end if
if(.not.allocated(nodes)) allocate(nodes(3,mxel))
if(.not.allocated(EMemo)) allocate(EMemo(mxel))
! element input
eleinput: do i = 1, numel
    read(100, *, iostat = ierr) no, n1, n2, n3
	if (ierr /=0) then
	  call error1
	  return
	endif
	nodes(1,i) = n1
	nodes(2,i) = n2
	nodes(3,i) = n3
end do eleinput
! fault
if(AppType == 1) then
   read(100,*, iostat = ierr) nfl
3000 If (GrayMenu) then
       mxfel = nfl
   else
     if (nfl > mxfel) then
        call error2('fault',nfl, mxfel)
	    ! NFL = 0
        ! return
        !
  	    ! modified on March 10, 2005, ask for node, element, fault size again
	    !
	    call InitArraySize
	    ! after call, mxnode, mxel, mxfel change
	    if(nfl <= mxfel) then
	      if(allocated(nodef)) deallocate(nodef)
          if(allocated(fdip)) deallocate(fdip)
		  if(allocated(offset)) deallocate(offset)
        else 
	      goto 3000
        endif
     end if
   end if
   ! allocate memory
   if (.not.allocated(nodef))  allocate(nodef(4, mxfel))
   if (.not.allocated(fdip))   allocate(fdip(2,mxfel))
   if (.not.allocated(offset)) allocate(offset(mxfel))
   if(ierr /= 0) then
     call error1
     return
   end if
   ! fault input
   faultinput: do i = 1, nfl
     read(100, *, iostat = ierr) no, n1, n2, n3, n4, dips1, dips2, offst
     if (ierr /=0) then
        call error1
	    return
     end if
     nodef(1,i) = n1
     nodef(2,i) = n2
     nodef(3,i) = n3
     nodef(4,i) = n4
     fdip(1,i) = dips1
     fdip(2,i) = dips2
     ! Dips are changed to [1 179] within orbwin if dips < 0
	 do j = 1, 2
        if (fdip(j,i) > 0.0) then
	    else
	       fdip(j,i) = 180.0 + fdip(j,i)
 	    end if
     end do
     offset(i) = offst
   end do faultinput
else ! Not thin-shell type, no fault allowed
endif
close(unit=100)
!
! successfully read (option to display confirmation dialog)
!msg0 = ' Successfully reading input file ' //feg_inp//' 'C
!msg1 = ' Successfully reading file 'C
!iret = messageboxqq(msg0, msg1, MB$ICONEXCLAMATION.OR.MB$OK)
!
! set global control switch
GridLoaded = .True.
!
return
end subroutine Read_Grid