# History file 'hes-inv.h' for smoothing of model HESS
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Dimension of array RAM in include file 'ram.inc' should be
# MRAM=14000000 or greater for this history file.
# Note also that this history file is rather time-consuming.
#
# Preparing data to be inverted
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Input files required
  chk.pl: "data/hes/"  "hes-mod1.dat"
  chk.pl: "data/hes/"  "hes-mod2.dat"
  chk.pl: "data/hes/"  "hess.dat"
  chk.pl: "model/"     "sob22.dat"
  chk.pl: "forms/"     "inv.cal"
  chk.pl: "forms/"     "addsob.cal"
  chk.pl: "forms/"     "sub.cal"
  chk.pl: "forms/"     "sqrt.cal"
  chk.pl: "forms/"     "eq.cal"
  chk.pl: "forms/"     "mul.cal"
  chk.pl: "forms/"     "add.cal"
# Original model
  MODEL='hess.dat'

# Gridding interfaces
  N1=64    N2=1
  D1=200   D2=1
  O1=9800  O2=0
  VEL='hes-s1.out'  ISRF=1  MPAR=0  ICB=' '
  grid: 'hes-inv.h' /
  VEL='hes-s2.out'  ISRF=2  MPAR=0  ICB=' '
  grid: 'hes-inv.h' /

# Gridding P wave velocities
  N1=202   N2=106
  D1=200   D2=200
  O1=0     O2=0
  VEL='hess-vp.out'  ISRF=0  MPAR=1  ICB='hes-ib.out'
  grid: 'hes-inv.h' /

# Computing P wave slowness
  CAL='inv.cal' GRD1='hess-vp.out' GRD2='hes-up.out'
  grdcal: 'hes-inv.h' /
# Average P wave slowness (harmonic average)
  N1NEW=1   N2NEW=1    D1NEW=1   D2NEW=1   O1NEW=0  O2NEW=0  GNORM=0
  GRD='hes-up.out' GRDNEW='hes-avup.out'
  grdnorm: 'hes-inv.h' /


# Gridding background P wave velocities (velocities in complex block 2)
# and slownesses
  VEL='hes-vpb.out'  ISRF=0  MPAR=1  ICB=              ICBEXT=2
  grid: 'hes-inv.h' /
  VEL=               ISRF=   MPAR=   ICB=              ICBEXT=
  CAL='inv.cal' GRD1='hes-vpb.out' GRD2='hes-upb.out'
  grdcal: 'hes-inv.h' /
# Average P wave velocity of the background (harmonic average)
  N1NEW=1   N2NEW=1    D1NEW=1   D2NEW=1   O1NEW=0  O2NEW=0  GNORM=0
  GRD='hes-vpb.out' GRDNEW='hes-avvp.out'
  grdnorm: 'hes-inv.h' /


# Average velocity on a coarse grid
  N1=202     N2=106    D1=200     D2=200     O1=0     O2=0
  N1NEW=67   N2NEW=106 D1NEW=600  D2NEW=200  O1NEW=200  O2NEW=0  GNORM=1
  GRD='hes-vpb.out' GRDNEW='hes-vpc.tmp'
  grdnorm: 'hes-inv.h' /
  N1=67   N2=106 D1=600  D2=200  O1=200  O2=0
  VEL=' '            ISRF=0  MPAR=1  ICB='hes-ibb.out' ICBEXT=2
  grid: 'hes-inv.h' /
  VEL=' '            ISRF=0  MPAR=1  ICB='hes-ibs.out' ICBEXT=1
  grid: 'hes-inv.h' /
  VEL=' '            ISRF=0  MPAR=1  ICB='hes-ibc.out' ICBEXT=
  grid: 'hes-inv.h' /
  VEL=               ISRF=   MPAR=   ICB=              ICBEXT=
  CAL='eq.cal' GRD1='hes-ibc.out' GRD2='hes-ibb.out' GRD3='hes-i.tmp'
  grdcal: 'hes-inv.h' /
  CAL='mul.cal' GRD1='hes-vpc.tmp' GRD2='hes-i.tmp' GRD3='hes-v1.tmp'
  grdcal: 'hes-inv.h' /
  VEL='hes-vpc.tmp'  ISRF=0  MPAR=1  ICB='hes-ibc.out'
  grid: 'hes-inv.h' /
  CAL='eq.cal' GRD1='hes-ibc.out' GRD2='hes-ibs.out' GRD3='hes-i.tmp'
  grdcal: 'hes-inv.h' /
  CAL='mul.cal' GRD1='hes-vpc.tmp' GRD2='hes-i.tmp' GRD3='hes-v2.tmp'
  grdcal: 'hes-inv.h' /
  CAL='add.cal' GRD1='hes-v1.tmp' GRD2='hes-v2.tmp' GRD3='hes-vpc.out'
  grdcal: 'hes-inv.h' /


# Inversion with smothing of surfaces and of the velocities
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Initial and updated models
  # MODEL='hes-mod1.dat'
    MODEL='hes-mod2.dat'
    MODOUT='hes-mod.out'
# Form of the files with matrices
  FORMM='UNFORMATTED'

# Calculating matrices for inversion (without densities)
  M1='m1.out'  M2='m2.out'
  MODL2=' '  MODSOB='modsobs.out'
  SOBOLEV='sob22.dat'  SOBW00=1 # minimizing interface curvature
  invsoft: 'hes-inv.h' /
  MODL2=' '  MODSOB=' '  SOBW00=

  MODL2=' '  MODSOB='modsobv.out'
  SOBOLEV='sob22.dat'    SOBW01=1 # minimizing second velocity derivatives
  invsoft: 'hes-inv.h' /
  MODL2=' '  MODSOB=' '  SOBW01=

  GM1='gm1.out'  GM2='gm2.out'  GM3=' '  DM1='dm1.out'
  N1=64    N2=1
  D1=200   D2=1
  O1=9800  O2=0
  GRD='hes-s1.out'  INDFUN=1               MPAR=0 ERRMUL=14.18  # SQRT(201)
  invpts: 'hes-inv.h' /
  M2IN='m2.out'
  GRD='hes-s2.out'  INDFUN=2               MPAR=0 ERRMUL=14.18  # SQRT(201)
  invpts: 'hes-inv.h' /
  N1=202  N2=106 D1=200  D2=200  O1=0    O2=0  # hes-mod1.out
  N1=67   N2=106 D1=600  D2=200  O1=200  O2=0  # hes-mod2.out
  GRD='hes-vp.out'   GRDICB='hes-ib.out'   MPAR=1 ERRMUL=146.33 #SQRT(21412)
  GRD='hes-vpc.out'  GRDICB='hes-ibc.out'  MPAR=1 ERRMUL=84.27  #SQRT(7102)
  invpts: 'hes-inv.h' /
  GRD=' '  M2IN=' ' GRDICB=

# Matrix operations
  N1=0  N2=1  N3=1  M1='m2.out'
  CAL='inv.cal' GRD1='dm1.out' GRD2='dm2.out'
  grdcal: 'hes-inv.h' /
  M1='m1.out'  M2='m2.out'  GM1='gm1.out'  DM1='dm2.out'  SM1='sm1.out'
  gmdmgmt: 'hes-inv.h' /
# Minimizing interface curvature
  N1=0  N2=0  N3=1  M1='m1.out'  M2='m1.out'
  SOBMUL=2500
  CAL='addsob.cal' GRD1='sm1.out' GRD2='modsobs.out' GRD3='sm2a.out'
  grdcal: 'hes-inv.h' /
# Minimizing second derivatives of P wave velocities
  N1=0  N2=0  N3=1  M1='m1.out'  M2='m1.out'
  SOBMUL=25000 0 1000 10000 17500 25000 35000 0 70000
  CAL='addsob.cal' GRD1='sm2a.out' GRD2='modsobv.out' GRD3='sm2.out'
  grdcal: 'hes-inv.h' /
  M1='m1.out'               SM1='sm2.out'  SM2='sm3.out'
  sminv: 'hes-inv.h' /
  M1='m2.out'  M2=' '       DM1='dm2.out'  GM1='gm2.out'  GM2='gm3.out'
  dmgm: 'hes-inv.h' /
  M1='m1.out'  M2='m2.out'  GM1='gm1.out'  GM2='gm3.out'  GM3='gm4.out'
  gmgm: 'hes-inv.h' /
  M1='m1.out'  M2=' '       SM1='sm3.out'  GM1='gm4.out'  GM2='gm5.out'
  smgm: 'hes-inv.h' /

# Updating the model
  M1='m1.out'  MODNEW='gm5.out'
  modmod: 'hes-inv.h' /

# Calculating components of the residual objective function
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

# Sobolev norm 'hes-sn1.out' of the surfaces curvature in the updated model
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  M1='m1.out'  M2=' '  SM1='modsobs.out'  GM1='gm5.out'  GM2='gm6.out'
  smgm: 'hes-inv.h' /
  M1=' '  M2='m1.out'  GM1='gm5.out'  GM2='gm6.out'  GM3='gm7.out'
  gmgm: 'hes-inv.h' /
  N1=0  N2=1  N3=1  M1=' '
  CAL='sqrt.cal' GRD1='gm7.out' GRD2='hes-sn1.out'
  FORMMW='FORMATTED'
  grdcal: 'hes-inv.h' /
  FORMMW=

# Sobolev norm 'hes-vn1.out' of the velocity in the updated model
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  M1='m1.out'  M2=' '  SM1='modsobv.out' GM1='gm5.out'  GM2='gm60.out'
  smgm: 'hes-inv.h' /
  M1=' '  M2='m1.out'  GM1='gm5.out'  GM2='gm60.out'  GM3='gm70.out'
  gmgm: 'hes-inv.h' /
  N1=0  N2=1  N3=1  M1=' '
  CAL='sqrt.cal' GRD1='gm70.out' GRD2='hes-vn1.out'
  FORMMW='FORMATTED'
  grdcal: 'hes-inv.h' /
  FORMMW=


# Standard velocity deviation 'hes-vd1.out' of the updated model
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  M1='m1.out'  M2='m2.out'  M3=
  GM1='gm1.out'  GM2='gm2.out'  GM3=' '  DM1='dm1.out'
  N1=202  N2=106 D1=200  D2=200  O1=0    O2=0  # hes-mod1.out
  N1=67   N2=106 D1=600  D2=200  O1=200  O2=0  # hes-mod2.out
  GRD='hes-vp.out'   GRDICB='hes-ib.out'   MPAR=1 ERRMUL=146.33 #SQRT(21412)
  GRD='hes-vpc.out'  GRDICB='hes-ibc.out'  MPAR=1 ERRMUL=84.27  #SQRT(7102)
  invpts: 'hes-inv.h' /
  GRD=' '           GRDICB=
  M1=' '  M2='m1.out'  M3='m2.out'
  GM1='gm5.out'  GM2='gm1.out'  GM3='gm8.out'
  gmgm: 'hes-inv.h' /
  M3=' '
  N1=0  N2=1  N3=1  M1='m2.out'
  CAL='sub.cal' GRD1='gm2.out' GRD2='gm8.out' GRD3='gm9.out'
  grdcal: 'hes-inv.h' /
  N1=0  N2=1  N3=1  M1='m2.out'
  CAL='inv.cal' GRD1='dm1.out' GRD2='dm2.out'
  grdcal: 'hes-inv.h' /
  M1=' '  M2='m2.out'  GM1='gm9.out'  DM1='dm2.out'  SM1='gm0.out'
  gmdmgmt: 'hes-inv.h' /
  N1=0  N2=1  N3=1  M1=' '
  CAL='sqrt.cal' GRD1='gm0.out' GRD2='hes-vd1.out' GRD3=
  FORMMW='FORMATTED'
  grdcal: 'hes-inv.h' /
  FORMMW=
#
# Standard surface deviation 'hes-sd1.out' of the updated model
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  M1='m1.out'  M2='m2.out'  M3=
  GM1='gm1.out'  GM2='gm2.out'  GM3=' '  DM1='dm1.out'
  N1=64    N2=1
  D1=200   D2=1
  O1=9800  O2=0
  GRD='hes-s1.out'  INDFUN=1      MPAR=0  ERRMUL=11.314  # SQRT(128)
  invpts: 'hes-inv.h' /
          M2IN='m2.out'
  GRD='hes-s2.out'  INDFUN=2      MPAR=0  ERRMUL=11.314  # SQRT(128)
  invpts: 'hes-inv.h' /
  GRD=' ' M2IN=
  M1=' '  M2='m1.out'  M3='m2.out'
  GM1='gm5.out'  GM2='gm1.out'  GM3='gm8.out'
  gmgm: 'hes-inv.h' /
  M3=' '
  N1=0  N2=1  N3=1  M1='m2.out'
  CAL='sub.cal' GRD1='gm2.out' GRD2='gm8.out' GRD3='gm9.out'
  grdcal: 'hes-inv.h' /
  N1=0  N2=1  N3=1  M1='m2.out'
  CAL='inv.cal' GRD1='dm1.out' GRD2='dm2.out'
  grdcal: 'hes-inv.h' /
  M1=' '  M2='m2.out'  GM1='gm9.out'  DM1='dm2.out'  SM1='gm0.out'
  gmdmgmt: 'hes-inv.h' /
  N1=0  N2=1  N3=1  M1=' '
  CAL='sqrt.cal' GRD1='gm0.out' GRD2='hes-sd1.out' GRD3=
  FORMMW='FORMATTED'
  grdcal: 'hes-inv.h' /
  FORMMW=

# Important output files
# ~~~~~~~~~~~~~~~~~~~~~~
# 'hess-vp.out'... Gridded velocity being fit.
# 'hes-s1.out'...  Gridded first surface being fit.
# 'hes-s2.out'...  Gridded second surface being fit.
# 'hes-mod.out'... Smoothed model.
# 'hes-vd1.out'... Standard velocity deviation of the model.
# 'hes-sd1.out'... Standard surfaces deviation of the model.
# 'hes-vn1.out'... Sobolev norm of the velocity in the model.
# 'hes-sn1.out'... Sobolev norm of the interfaces in the model.