# 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.
#

# 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/"     "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:
  VEL='hes-s2.out'  ISRF=2  MPAR=0  ICB=' '
  grid:

# 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:

# Computing P wave slowness
  CAL='inv.cal' GRD1='hess-vp.out' GRD2='hes-up.out'
  grdcal:
# 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:


# Gridding background P wave velocities (velocities in complex block 2)
# and slownesses
  VEL='hes-vpb.out'  ISRF=0  MPAR=1  ICB=              ICBEXT=2
  grid:
  VEL=               ISRF=   MPAR=   ICB=              ICBEXT=
  CAL='inv.cal' GRD1='hes-vpb.out' GRD2='hes-upb.out'
  grdcal:
# 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:


# 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:
  N1=67   N2=106 D1=600  D2=200  O1=200  O2=0
  VEL=' '            ISRF=0  MPAR=1  ICB='hes-ibb.out' ICBEXT=2
  grid:
  VEL=' '            ISRF=0  MPAR=1  ICB='hes-ibs.out' ICBEXT=1
  grid:
  VEL=' '            ISRF=0  MPAR=1  ICB='hes-ibc.out' ICBEXT=
  grid:
  VEL=               ISRF=   MPAR=   ICB=              ICBEXT=
  CAL='eq.cal' GRD1='hes-ibc.out' GRD2='hes-ibb.out' GRD3='hes-i.tmp'
  grdcal:
  CAL='mul.cal' GRD1='hes-vpc.tmp' GRD2='hes-i.tmp' GRD3='hes-v1.tmp'
  grdcal:
  VEL='hes-vpc.tmp'  ISRF=0  MPAR=1  ICB='hes-ibc.out'
  grid:
  CAL='eq.cal' GRD1='hes-ibc.out' GRD2='hes-ibs.out' GRD3='hes-i.tmp'
  grdcal:
  CAL='mul.cal' GRD1='hes-vpc.tmp' GRD2='hes-i.tmp' GRD3='hes-v2.tmp'
  grdcal:
  CAL='add.cal' GRD1='hes-v1.tmp' GRD2='hes-v2.tmp' GRD3='hes-vpc.out'
  grdcal:

# 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'
  MATFORM='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:
  MODL2=' '  MODSOB=' '  SOBW00=

  MODL2=' '  MODSOB='modsobv.out'
  SOBOLEV='sob22.dat'    SOBW01=1 # minimizing second velocity derivatives
  invsoft:
  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:
  M2IN='m2.out'
  GRD='hes-s2.out'  INDFUN=2               MPAR=0 ERRMUL=14.18  # SQRT(201)
  invpts:
  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:
  GRD=' '  M2IN=' ' GRDICB=

# Matrix operations
  MATIN1='dm1.out'                       MATOUT='dm2.out'
  MATFUN='inv'
  matfun:
  MATIN1='dm2.out'   MATIN2='gm1.out'    MATOUT='dm2gm1.out'
  SYMMETRY=      MATT1=  MATT2=1
  matmul:
  MATIN1='gm1.out'   MATIN2='dm2gm1.out' MATOUT='sm1.out'
  SYMMETRY='sym' MATT1=  MATT2=
  matmul:
# Minimizing interface curvature
  MATIN1='sm1.out'   MATIN2='modsobs.out' MATOUT='sm2a.out'
  COEF1=  COEF2=6250000
  matlin:
# Minimizing second derivatives of P wave velocities
  MATIN1='sm2a.out'  MATIN2='modsobv.out' MATOUT='sm2.out'
  COEF1=  COEF2=625000000
  matlin:
  MATIN1='sm2.out'                       MATOUT='sm3.out'
  matinv:
  MATIN1='dm2.out'   MATIN2='gm2.out'    MATOUT='gm3.out'
  SYMMETRY=' '   MATT1=  MATT2=
  matmul:
  MATIN1='gm1.out'   MATIN2='gm3.out'    MATOUT='gm4.out'
  SYMMETRY=' '   MATT1=  MATT2=
  matmul:
  MATIN1='sm3.out'   MATIN2='gm4.out'    MATOUT='gm5.out'
  SYMMETRY=' '   MATT1=  MATT2=
  matmul:

# Updating the model
  M1='m1.out'  MODNEW='gm5.out'
  modmod:

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

# Sobolev norm 'hes-sn1.out' of the surfaces curvature in the updated model
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  MATIN1='modsobs.out'   MATIN2='gm5.out'    MATOUT='gm6.out'
  SYMMETRY=' '   MATT1=  MATT2=
  matmul:
  MATIN1='gm5.out'       MATIN2='gm6.out'    MATOUT='gm7.out'
  SYMMETRY='diag' MATT1=1 MATT2=0
  matmul:
  MATIN1='gm7.out'                       MATOUT='hes-sn1.out'
  SPARSE=-1 MATFUN='sqrt'                MATFORM='formatted'
  matfun:
                                         MATFORM='unformatted'

# Sobolev norm 'hes-vn1.out' of the velocity in the updated model
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  MATIN1='modsobv.out'   MATIN2='gm5.out'    MATOUT='gm60.out'
  SYMMETRY=' '   MATT1=  MATT2=
  matmul:
  MATIN1='gm5.out'       MATIN2='gm60.out'    MATOUT='gm70.out'
  SYMMETRY='diag' MATT1=1 MATT2=0
  matmul:
  MATIN1='gm70.out'                      MATOUT='hes-vn1.out'
  SPARSE=-1 MATFUN='sqrt'                MATFORM='formatted'
  matfun:
                                         MATFORM='unformatted'

# 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:
  GRD=' '           GRDICB=
  MATIN1='gm1.out'   MATIN2='gm5.out'    MATOUT='gm8.out'
  SYMMETRY=      MATT1=1 MATT2=
  matmul:
  MATIN1='gm2.out'   MATIN2='gm8.out' MATOUT='gm9.out'
  COEF1=  COEF2=-1
  matlin:
  MATIN1='dm1.out'                       MATOUT='dm2.out'
  MATFUN='inv'
  matfun:
  MATIN1='dm2.out'   MATIN2='gm9.out'    MATOUT='dm2gm9.out'
  SYMMETRY=       MATT1=  MATT2=
  matmul:
  MATIN1='gm9.out'   MATIN2='dm2gm9.out' MATOUT='gm0.out'
  SYMMETRY='diag' MATT1=1 MATT2=
  matmul:
  MATIN1='gm0.out'                       MATOUT='hes-vd1.out'
  SPARSE=-1 MATFUN='sqrt'                MATFORM='formatted'
  matfun:
                                         MATFORM='unformatted'

#
# 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:
          M2IN='m2.out'
  GRD='hes-s2.out'  INDFUN=2      MPAR=0  ERRMUL=11.314  # SQRT(128)
  invpts:
  GRD=' ' M2IN=
  MATIN1='gm1.out'   MATIN2='gm5.out'    MATOUT='gm8.out'
  SYMMETRY=      MATT1=1 MATT2=
  matmul:
  MATIN1='gm2.out'   MATIN2='gm8.out' MATOUT='gm9.out'
  COEF1=  COEF2=-1
  matlin:
  MATIN1='dm1.out'                       MATOUT='dm2.out'
  MATFUN='inv'
  matfun:
  MATIN1='dm2.out'   MATIN2='gm9.out'    MATOUT='dm2gm9.out'
  SYMMETRY=       MATT1=  MATT2=
  matmul:
  MATIN1='gm9.out'   MATIN2='dm2gm9.out' MATOUT='gm0.out'
  SYMMETRY='diag' MATT1=1 MATT2=
  matmul:
  MATIN1='gm0.out'                       MATOUT='hes-sd1.out'
  SPARSE=-1 MATFUN='sqrt'                MATFORM='formatted'
  matfun:
                                         MATFORM='unformatted'

# 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.ou@'... Standard velocity deviation of the model.
# 'hes-sd1.ou@'... Standard surfaces deviation of the model.
# 'hes-vn1.ou@'... Sobolev norm of the velocity in the model.
# 'hes-sn1.ou@'... Sobolev norm of the interfaces in the model.