!
!
! Subroutines: Output_Grid, AutoSave
! 
!
subroutine Output_Grid
!
! output grid infor. in the format of .FEG
! input: Linklists of node, element, fault
!
use dflib
use Global

implicit none
character(len=512) msg0, msg1
integer(4) iRet, ierr
integer(4) i, LRi, no, n1, n2, n3, n4
real lon, lat, elev, q, hc,hm, dips1,dips2, offst, density_anomaly, cooling_curvature
real tempvec(3)

msg0 = ' 'C
msg1 = ' 'C
open(unit = 200, file = feg_sav, status = 'new', action='write', iostat = ierr)
if (ierr /= 0) then
  msg0 = 'Specified file already exists. Do you want to overwrite it?'C
  msg1 = 'Error opening file'C
  iRet = messageboxqq(msg0, msg1, MB$ICONEXCLAMATION.OR.MB$YESNOCANCEL)
  if (iRet == MB$IDYES) then
     open(unit = 200, file = feg_sav, status = 'old', iostat = ierr)
	 rewind(200)
  else
     return
  endif
endif
! write out into file
write(200, '(A80)') new_FEG_title_line
! node
nRealN = numNod
nFakeN = 0
n1000 = MAX(1000000, numNod)
brief =.TRUE.
!
write(200, 2) numnod, nrealn, nfaken, n1000, brief
2  FORMAT(4I8,L8,'  (numNod, nRealN, nFakeN, n1000, brief)')
outputnode: do i = 1, numnod
    tempvec(1) = nodeABG(1,i)
	tempvec(2) = nodeABG(2,i)
	tempvec(3) = nodeABG(3,i)
	call ABG2lonlat(tempvec, lon,lat)
	lon = lon * 57.2957795
	lat = lat * 57.2957795
	elev = eqcm(1,i)
	q = eqcm(2,i)
	hc = eqcm(3,i)
	hm = eqcm(4,i)
    density_anomaly = eqcm(5,i)
    cooling_curvature = eqcm(6,i)
    IF (OrbData5) THEN ! OrbData5 format: 6 real variables per node
	    write(200, 3) i, lon, lat, elev, q, hc, hm, density_anomaly, cooling_curvature
    ELSE ! old OrbData format; only 4 real variables per node
	    write(200, 3) i, lon, lat, elev, q, hc, hm
    END IF
end do outputnode
3  format(I8, 2F11.5, 6ES10.2)
! element
write(200,4) numel
4  format(I10,' (numEl = number of triangular continuum elements)')
if (numel > 0) then
    outputele: do i = 1, numel
       n1 = nodes(1,i)
       n2 = nodes(2,i)
       n3 = nodes(3,i)
       if (appType == 3) then ! Restore3+ format, where 3 per-element data MAY be present:
           if ((element_data(2, i) > 0.0).or.(element_data(3, i) > 0.0)) then ! all 3 per-element data are present
               write(200, '(I8, 3I8, ES9.1, F7.1, ES9.1)') i, n1, n2, n3, element_data(1:3, i) ! write out all 3 (young element_mu_, mu_switch, old element_mu_)
           else if (element_data(1, i) > 0.0) then ! one value is present (time-independent mu_)
               write(200, '(I8, 3I8, ES9.1)') i, n1, n2, n3, element_data(1, i) ! write out mu_ only
           else ! although per-element data are ALLOWED, none are actually present
               write(200, '(I8, 3I8)') i, n1, n2, n3
           end if
       ELSE IF (AppType == 1) THEN ! Shells mode; support v5.0+
           LRi = continuum_LRi(i)
           IF (LRi > 0) THEN
               WRITE (200, "(I8, 3I8, ' LR', I8)") i, n1, n2, n3, LRi
           ELSE
               WRITE (200, "(I8, 3I8)") i, n1, n2, n3
           END IF
       else ! earlier .FEG format, with no per-element data allowed:
           write(200, '(I8, 3I8)') i, n1, n2, n3
       end if
    end do outputele
end if ! numel > 0
! faults
write(200,5) nfl
5 format (I10,' (nFl =  number of linear fault elements)')
if (nfl > 0) then
    outputfault: do i = 1, nfl
       n1 = nodef(1,i)
       n2 = nodef(2,i)
       n3 = nodef(3,i)
       n4 = nodef(4,i)
       dips1 = fdip(1,i)
       dips2 = fdip(2,i)
       offst = offset(i)
       !
       ! set dip from [1 ~ 179] to [-89 ~ 89]
       ! but within OrbWin, dips are in [1 ~ 179], when reading grid, dips in [-89 ~ 89]
       ! are changed back to [1 ~ 179], see ReadGrid.f90 for details
       !
       if (dips1 > 90.0) dips1 = dips1 - 180.0
       if (dips2 > 90.0) dips2 = dips2 - 180.0
       IF (AppType == 1) THEN ! Shells mode; support v5.0+
           LRi = fault_LRi(i)
           IF (LRi > 0) THEN
               WRITE (200, "(I8, 4I8, 2F6.1, ES10.2, ' LR', I8)") i, n1, n2, n3, n4, dips1, dips2, offst, LRi
           ELSE
               WRITE (200, "(I8, 4I8, 2F6.1, ES10.2)") i, n1, n2, n3, n4, dips1, dips2, offst
           END IF
       ELSE ! other value of AppType
           WRITE (200, "(I8, 4I8, 2F6.1, ES10.2)") i, n1, n2, n3, n4, dips1, dips2, offst
       END IF
    end do outputfault
end if ! nfl > 0
close(unit = 200)
! success output
msg0 = 'Successfully saved grid as ' //feg_sav//' 'C
msg1 = 'Successfully saved'C
iRet = messageboxqq(msg0, msg1, MB$ICONEXCLAMATION.OR.MB$OK)
return
end subroutine Output_Grid
!
!
!
!
subroutine AutoSave(filesave)
use dflib
use Global
implicit none
character*(*), intent(in) :: filesave
integer(4) iRet, ierr
integer(4) i, no, n1, n2, n3, n4
real lon, lat, elev, q, hc,hm, dips1,dips2, offst, density_anomaly, cooling_curvature
real tempvec(3)

open(unit = 400, file = TRIM(filesave), status = 'unknown', action='write', iostat = ierr)
! write out into file
write(400, "(A)") TRIM(filesave)
! node
brief =.false.
! consistent with old Orbweave.exe
nrealn = numnod
! if (nrealn == 0) nrealn = numnod
if (n1000 == 0) n1000 = 100000
!
write(400, 2) numnod, nrealn, nfaken, n1000, brief
2  format(4I8,L8,' (numNod, nRealN, nFakeN, n1000, brief)')
outputnode: do i = 1, numnod
    tempvec(1) = nodeABG(1,i)
	tempvec(2) = nodeABG(2,i)
	tempvec(3) = nodeABG(3,i)
	call ABG2lonlat(tempvec, lon,lat)
	lon = lon * 57.2957795
	lat = lat * 57.2957795
	elev = eqcm(1,i)
	q = eqcm(2,i)
	hc = eqcm(3,i)
	hm = eqcm(4,i)
    density_anomaly = eqcm(5,i)
    cooling_curvature = eqcm(6,i)
    IF (OrbData5) THEN ! OrbData5 format: 6 real values per node
	    write(400, 3) i, lon, lat, elev, q, hc, hm, density_anomaly, cooling_curvature
    ELSE ! old OrbData format: only 4 real values per node
	    write(400, 3) i, lon, lat, elev, q, hc, hm
    END IF
end do outputnode
3  format(I8, 2F11.5, 6ES10.2)
! element
write(400,4) numel
4  format(I10,' (numEl = number of triangular continuum elements)')
if(numel > 0) then
    outputele: do i = 1, numel
       n1 = nodes(1,i)
       n2 = nodes(2,i)
       n3 = nodes(3,i)
       if (appType == 3) then ! Restore3+ format, where 3 per-element data MAY be present:
           if ((element_data(2, i) > 0.0).or.(element_data(3, i) > 0.0)) then ! all 3 per-element data are present
               write(400, '(I8, 3I8, ES9.1, F7.1, ES9.1)') i, n1, n2, n3, element_data(1:3, i) ! write out all 3 (young element_mu_, mu_switch, old element_mu_)
           else if (element_data(1, i) > 0.0) then ! one value is present (time-independent mu_)
               write(400, '(I8, 3I8, ES9.1)') i, n1, n2, n3, element_data(1, i) ! write out mu_ only
           else ! although per-element data are ALLOWED, none are actually present
               write(400, '(I8, 3I8)') i, n1, n2, n3
           end if
       else ! earlier .FEG format, with no per-element data allowed:
           write(400, '(I8, 3I8)') i, n1, n2, n3
       end if
    end do outputele
end if
! fault
write(400,5) nfl
5 format (I10,' (nFl = number of linear fault elements)')
if(nfl > 0) then
    outputfault: do i = 1, nfl
       n1 = nodef(1,i)
       n2 = nodef(2,i)
       n3 = nodef(3,i)
       n4 = nodef(4,i)
       dips1 = fdip(1,i)
       dips2 = fdip(2,i)
       if (dips1 > 90.0) dips1 = dips1 - 180.0
       if (dips2 > 90.0) dips2 = dips2 - 180.0
       offst = offset(i)
       write(400, '(I8, 4I8, 2F6.1, ES10.2)') i, n1, n2, n3, n4, dips1, dips2, offst
    end do outputfault
end if
close(unit = 400)
! success output
msg = ' Successfully saved backup grid in ' // TRIM(filesave)
call setstatusbar(1, TRIM(msg) // CHAR(0))
return
end subroutine AutoSave