# History file 'plu-inv.h' for building and smoothing model Pluto 1.5
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Dimension of array RAM in include file 'ram.inc' should be
# MRAM=15200000 or greater for this history file.

# Input files required
# ~~~~~~~~~~~~~~~~~~~~
 #chk.pl: "data/plu/" "plu-mod.dat"
 #chk.pl: "model/"    "sob22.dat"
 #chk.pl: "forms/"    "inv.cal"
 #chk.pl: "forms/"    "addsob.cal"
 #chk.pl: "forms/"    "sqrt.cal"
 #chk.pl: "forms/"    "echo.pl"
 #chk.pl: "forms/"    "copy.pl"
 #chk.pl: ""          "plu-vp.grd"
# File 'plu-vp.grd' should be extracted from the ZIPped file
# 'plu-vp.zip' located in directory 'bin'


# Preparing data to be inverted
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Data points at interfaces
# ~~~~~~~~~~~~~~~~~~~~~~~~~
  N1=1201 N2=6960 D1=25. D2=25. O1=0. O2=0.
  VALUE=14000. MODE=-1
  GRD='plu-vp.grd'
  PTS='p-d-i1.out' PTS1='p-d-i2.out'
  N1NEW=1000
  N2NEW=873 ND2=3
  grdiso:

  PTS='p-d-i3.out' PTS1='p-d-i4.out'
  NO2=2640  N2NEW=960 ND2=3
  grdiso:

  PTS='p-d-i5.out' PTS1=' '
  NO1=601   N1NEW=600
  NO2=      N2NEW=1   ND2=
  grdiso:
#
# Data for P-wave velocity
# ~~~~~~~~~~~~~~~~~~~~~~~~
  N1NEW=50 N2NEW=87 D1NEW=600. D2NEW=2000. O1NEW=50. O2NEW=1000.
  GRD='plu-vp.grd' GRDNEW='p-d-v.out'
  grdnew:


# Average slowness and velocity for estimations of SOBMUL
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  copy.pl: "plu-vp.grd" "plu-slop.out"
  echo.pl: "@1=1./@1"   ">  cal.tmp"
  CAL='cal.tmp' GRD1='plu-slop.out'
  grdcal:
  N1NEW=1  N2NEW=1          GNORM=2
  GRD='plu-slop.out' GRDNEW='p-avslo.out'      grdnorm:
  GRD='plu-vp.grd'   GRDNEW='p-avvp.out'       grdnorm:
  GRD=               GRDNEW=

# Inversion
# ~~~~~~~~~
# Form of the files with matrices
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  FORMM='unformatted'

# Inversion of the interfaces
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ICLASS=1

# Initial model
  MODEL='plu-mod.dat'

# Calculating matrices for inversion (without densities)
  M1='m1.out'  M2='m2.out' SOBOLEV='sob22.dat'
  MODSOB='modsobs.out' SOBW00=1 # minimizing interface curvature
  invsoft:
  MODSOB=              SOBW00=
  GM1='gm1.out'  GM2='gm2.out'  GM3=' '  DM1='dm1.out'
  PTS='p-d-i1.out' INDFUN=1  ERRMUL=29.55  invpts:  #   873 points
  M2IN='m2.out'                                     #
  PTS='p-d-i2.out' INDFUN=2  ERRMUL=29.55  invpts:  #   873 points
  PTS='p-d-i3.out' INDFUN=3  ERRMUL=28.44  invpts:  #   809 points
  PTS='p-d-i4.out' INDFUN=4  ERRMUL=28.44  invpts:  #   809 points
  PTS='p-d-i5.out' INDFUN=5  ERRMUL=1.     invpts:  #     1 points
  M2IN=     PTS=     INDFUN=   ERRMUL=              #  3365 total

# Matrix operations
  N1=0  N2=1  N3=1  M1='m2.out'
  CAL='inv.cal' GRD1='dm1.out' GRD2='dm2.out'
  grdcal:
  M1='m1.out'  M2='m2.out'  GM1='gm1.out'  DM1='dm2.out'  SM1='sm1.out'
  gmdmgmt:
  N1=0  N2=0  N3=1  M1='m1.out' M2='m1.out'
  CAL='addsob.cal'  GRD1='sm1.out'  GRD2='modsobs.out' GRD3='sm2.out'
  SOBMUL=30000 25000 50000 2500 0
  grdcal:
  SOBMUL=
  M1='m1.out'               SM1='sm2.out'  SM2='sm3.out'
  sminv:
  M1='m2.out'  M2=' '       DM1='dm2.out'  GM1='gm2.out'  GM2='gm3.out'
  dmgm:
  M1='m1.out'  M2='m2.out'  GM1='gm1.out'  GM2='gm3.out'  GM3='gm4.out'
  gmgm:
  M1='m1.out'  M2=' '       SM1='sm3.out'  GM1='gm4.out'  GM2='gm5.out'
  smgm:

# Updating the model
  M1='m1.out'  MODNEW='gm5.out'  MODOUT='plu-mod.tmp'
  modmod:
  M1= M2=

# Indices of complex blocks relevant to the velocities
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Velocities greater than 14700 => ICB=1; other velocities => ICB=2
  N1=50 N2=87 D1=600. D2=2000. O1=50. O2=1000.
  echo.pl: "A=@1/14700."     ">  cal.tmp"
  echo.pl: "A=INT(A)"        ">> cal.tmp"
  echo.pl: "A=2-A"           ">> cal.tmp"
  echo.pl: "@2=NINT(A)"      ">> cal.tmp"
  CAL='cal.tmp' GRD1='p-d-v.out'   GRD2='p-d-icb.out'
  grdcal:


# Inversion of velocities
# ~~~~~~~~~~~~~~~~~~~~~~~
  MODEL='plu-mod.tmp'    MODIN='plu-mod.tmp'
  ICLASS=2

# Calculating matrices for inversion
  M1='m1.out'  M2='m2.out'
  MODSOB='modsobv.out' SOBW01=1 # minimizing second velocity derivatives
  invsoft:
  MODSOB=               SOBW01=
  GM1='gm1.out'  GM2='gm2.out'  GM3=' '  DM1='dm1.out'
  MPAR=1 POWERM=1 ERRMUL=65.955 # 4350 points
  GRD='p-d-v.out' GRDICB='p-d-icb.out'
  invpts:
  GRD= GRDICB= MPAR= POWERM= ERRMUL=

# Matrix operations
  N1=0  N2=1  N3=1  M1='m2.out'
  CAL='inv.cal' GRD1='dm1.out' GRD2='dm2.out'
  grdcal:
  M1='m1.out'  M2='m2.out'  GM1='gm1.out'  DM1='dm2.out'  SM1='sm1.out'
  gmdmgmt:

  N1=0  N2=0  N3=1  M1='m1.out' M2='m1.out'
  CAL='addsob.cal'  GRD1='sm1.out'  GRD2='modsobv.out'  GRD3='sm2.out'
  SOBMUL=200000 125000 1000000 1000 0 75000 150000 50000 100000 500000
  grdcal:
  SOBMUL=
  M1='m1.out'               SM1='sm2.out'  SM2='sm3.out'
  sminv:
  M1='m2.out'  M2=' '       DM1='dm2.out'  GM1='gm2.out'  GM2='gm3.out'
  dmgm:
  M1='m1.out'  M2='m2.out'  GM1='gm1.out'  GM2='gm3.out'  GM3='gm4.out'
  gmgm:
  M1='m1.out'  M2=' '       SM1='sm3.out'  GM1='gm4.out'  GM2='gm5.out'
  smgm:

# Updating the model
  M1='m1.out'  MODNEW='gm5.out'
  MODOUT='plu-mod.out'
  modmod:


# Sobolev norm of the velocity in the updated model
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  M1='m1.out'  M2=' '  SM1='modsobv.out' GM1='gm5.out'  GM2='gm6.out'
  smgm:
  M1=' '  M2='m1.out'  GM1='gm5.out'  GM2='gm6.out'  GM3='gm7.out'
  gmgm:
  N1=0  N2=1  N3=1  M1=' '
  CAL='sqrt.cal' GRD1='gm7.out' GRD2='p-snv.out' FORMMW='formatted'
  grdcal:


# Important output files
# ~~~~~~~~~~~~~~~~~~~~~~
# 'p-d-v.out' ...  Gridded velocity being fit.
# 'p-d-i1.out',...,'p-d-i5.out' ...  Points at interfaces being fit.
# 'plu-mod.out'... Inverted model.
# 'p-snv.out'...   Sobolev norm of the velocity in the model.