!==============================================================================
!  Reversible Jump MCMC Sampling with parallel tempering for SWD modeling
!   (Last modified: Surya Pachhai
!                   29 Nov. 2015, Bullard Lab., University of Cambridge)
!                   17 August 2016, ICTP
!  This version of code allows to interact different cpu s. all cpus are 
!  assigned different temperature. 
!  NOTE: This code  inverts for number of layers, thicknesses, S velocity. 
!  Vp/Vs ratio is used to estimate Vp, but kept fixed as there
!  is no sensitivity on Vp/Vs 
!------------------------------------------------------------------------------
!I would highly appreciate if you could reference the following papers
!1) Pachhai., S., M. Li, M.S. Thorne, J. Dettmer, and H. Tkalcic (2022), Internal structure of ultralow-velicity zones consistent with origin from a basal magma ocean, Nature Geoscience, 15(1), 79-84.
!2) Manu-Marfo, D., Aoudia, A., Pachhai, S., and R. Kherchouche (2019), 3D shear wave velocity model of the crust and uppermost mantle beneath the Tyrrhenian basin and margins, Scientific Reports, 9:3609, https://doi.org/10.1038/s41598-019-40510-z. 
!3) Pachhai, S., J. Dettmer, and H. Tkalčić (2015), Ultra-low velocity zones beneath the Philippine and Tasman Seas revealed by trans-dimensional Bayesian waveform inversion, Geophysical Journal International, 203 (2), 1302-1318, doi:10.1093/gji/ggv368. 
!4) Pachhai, S., H. Tkalčić, and J. Dettmer (2014), Bayesian inference for ultralow velocity zones in the Earth’s lowermost mantle: complex ULVZ beneath the east of Philippines, Journal of Geophysical Research, 11(119), 8346-8365, doi: 10.1002/2014JB011067.  
!5) Thapa, H.R., S. Pachhai, A. Aoudia, D. Manu-Marfo, K. Priestley, and S. Mitra (2023), Crustal structure and discontinuities beneath the Nepal Himalaya using seismic ambient noise and teleseismic P wave coda autocorrelation, Journal of Geophysical Research,
!==============================================================================
!To do : still need to use sdave for rf, pv, and gv
MODULE DATA_TYPE
IMPLICIT NONE
INTEGER(KIND=4), PARAMETER               :: IB=4, RP=KIND(0.D0), SP=KIND(0.0)
REAL(KIND=SP),   PARAMETER               :: PI2  = 3.141592653589793238462643383279502884197_SP
END MODULE DATA_TYPE
!==============================================================================
MODULE RJMCMC_COM
USE MPI
USE DATA_TYPE
IMPLICIT NONE
!==============================================================================
!!
!! General switches
!!
!==============================================================================
INTEGER(KIND=IB), PARAMETER              :: ICOV         = 1      ! 0 = Sample implicit over sigma
                                                                  ! 1 = Sample over sigma
INTEGER(KIND=IB), PARAMETER              :: IPV          = 1      ! 1 = Phase velocity inversion
INTEGER(KIND=IB), PARAMETER              :: IGV          = 1      ! 1 = group velocity inversion
INTEGER(KIND=IB), PARAMETER              :: IHRMN        = 1      ! 1 to call Harman's code
INTEGER(KIND=IB), PARAMETER              :: IMAP         = 0      ! WRITE REPLICA AND EXIT
INTEGER(KIND=IB), PARAMETER              :: ISETSEED     = 0      ! Fix the random seed 
INTEGER(KIND=IB), PARAMETER              :: ISMPPRIOR    = 0      ! CHECK DEATH/BIRTH
!==============================================================================
!!
!! Sampling move related switches
!!
!==============================================================================

INTEGER(KIND=IB), PARAMETER              :: IEXCHANGE    = 1      ! 1 = turn on exchange moves (parallel tempering)
!=============================================================================
!!
!!  AUV DATA RJMCMC trial ping
!!
!==============================================================================
INTEGER(KIND=IB), PARAMETER              :: IAR          = 0      ! 1 = Use Autoregressive error model
INTEGER(KIND=IB), PARAMETER              :: NPTS         = 2048   ! Maximum number of parameters in  RF
INTEGER(KIND=IB)                         :: NDAT                  ! Actual No. data
INTEGER(KIND=IB), PARAMETER              :: NLMX         = 35     ! Max number of layers #changed for 35
INTEGER(KIND=IB), PARAMETER              :: NPL          = 2      ! Number parameters per layer (thickness, and S velocity, and Vp/Vs ratio)
INTEGER(KIND=IB), PARAMETER              :: IMISC        = 0      ! Switch to turn on/off misc sampling
INTEGER(KIND=IB), PARAMETER              :: NMISC        = 1      ! Number misc parameters (e.g., alpha for Gaussian filter)
CHARACTER(len=64)                        :: filebasefile           = 'filebase.txt'
INTEGER(KIND=IB), PARAMETER              :: NWF          = 1      ! Single station waveform

!==============================================================================
!!
!!  Prior variables and good seeding model
!!
!==============================================================================
REAL(KIND=SP),DIMENSION(NPL)             :: minlim     = 0.       ! Minimum limit for parameters
REAL(KIND=SP),DIMENSION(NPL)             :: maxlim     = 0.       ! Maximum limit for parameters
INTEGER(KIND=IB)                         :: kmin       = 1        ! Min number of interfaces
INTEGER(KIND=IB)                         :: kmax       = NLMX     ! Max number of interfaces
REAL(KIND=SP),PARAMETER                  :: hmin       = 0.5_SP   ! Min allowed layer thickness in km
REAL(KIND=SP),PARAMETER                  :: Zmin       = 1.0      ! Minimum height from the CMB in km
REAL(KIND=SP),PARAMETER                  :: fact       = 1.00_SP  ! Factor for rotated space perturbation
REAL(KIND=SP),DIMENSION(NPL)             :: maxpert    = 0.       ! Factor for perturbation width defined later
REAL(KIND=SP),DIMENSION(NPL)             :: pertsd     = 0.       ! Perturbation width used for proposal
REAL(KIND=SP),DIMENSION(NPL)             :: pertsdsc   = (/1050.0, 50.0/)
!-----------------------------------------------------------------------
!increase in these parameters increases the acceptance rate as the width
!decreases
!-----------------------------------------------------------------------
!==============================================================================
!!
!!  Autoregressive model prior variables:
!!
!==============================================================================

INTEGER(KIND=IB),PARAMETER               :: NARFP     = 1          ! Order of each Autoregressive process
REAL(KIND=SP)                            :: armx      = 0.2_SP     ! Max AR and ARI model range (R units)
REAL(KIND=SP),DIMENSION(NARFP*NWF)       :: minlimar  = 0.0_SP    ! Minimum prior bounds for ar parameters
REAL(KIND=SP),DIMENSION(NARFP*NWF)       :: maxlimar  = 0.999_SP   ! Maximum prior bounds for ar parameters
REAL(KIND=SP),DIMENSION(NARFP*NWF)       :: maxpertar = 0._SP
REAL(KIND=SP),DIMENSION(NARFP*NWF)       :: pertarsd  = 30._SP     ! 
REAL(KIND=SP),DIMENSION(NARFP*NWF)       :: pertarsdsc= 18._SP
!==============================================================================
!!
!!  MISC parameter prior bounds etc (in log space)
!!
!==============================================================================

REAL(KIND=SP), DIMENSION(NARFP*NWF)      :: minlimmisc = 0._SP
REAL(KIND=SP), DIMENSION(NARFP*NWF)      :: maxlimmisc = 0._SP
REAL(KIND=SP), DIMENSION(NARFP*NWF)      :: maxpertmisc= 0._SP
REAL(KIND=SP), DIMENSION(NARFP*NWF)      :: pertmiscsd = 0._SP
REAL(KIND=SP), DIMENSION(NARFP*NWF)      :: pertmiscsdsc= 40._SP
!==============================================================================
!!
!!  Standard deviation prior variables:
!!
!==============================================================================

REAL(KIND=SP), DIMENSION(NWF)            :: minlimsd  = 0._SP
REAL(KIND=SP), DIMENSION(NWF)            :: maxlimsd  = 0._SP
REAL(KIND=SP), DIMENSION(NWF)            :: maxpertsd = 0._SP
REAL(KIND=SP), DIMENSION(NWF)            :: pertsdsd  = 0._SP
REAL(KIND=SP), DIMENSION(NWF)            :: pertsdsdsc= 10._SP
REAL(KIND=SP), PARAMETER                 :: hmx       = 300.0_SP    !! Max thickness of the layer below the surface # modified from 200 by Daniel on 20 July 2017
REAL(KIND=SP), PARAMETER                 :: vpvs      = 1.732_SP 
CHARACTER(len=64)                        :: filebase
INTEGER(KIND=IB)                         :: filebaselen
CHARACTER(LEN=100)                       :: infile_pv
CHARACTER(LEN=100)                       :: infile_gv
CHARACTER(LEN=64)                        :: logfile
CHARACTER(LEN=64)                        :: seedfile
CHARACTER(LEN=64)                        :: mapfile
CHARACTER(LEN=64)                        :: covfile
CHARACTER(LEN=64)                        :: obsfile_pv
CHARACTER(LEN=64)                        :: obsfile_gv
CHARACTER(LEN=64)                        :: repfile_pv
CHARACTER(LEN=64)                        :: repfile_gv
CHARACTER(LEN=64)                        :: sdfile
CHARACTER(LEN=64)                        :: samplefile
CHARACTER(LEN=64)                        :: lincovfile
CHARACTER(LEN=64)                        :: stepsizefile
CHARACTER(LEN=64)                        :: refmodel
INTEGER(KIND=IB),PARAMETER               :: ustat   = 57         !unit for status file
INTEGER(KIND=IB),PARAMETER               :: ulog    = 58         !unit for log file
INTEGER(KIND=IB),PARAMETER               :: usample = 59         !unit for sample file
INTEGER(KIND=IB),PARAMETER               :: ustep   = 60         !unit for stepsize file
!==============================================================================
!!
!! Parallel Tempering parameters
!!
!==============================================================================

REAL(KIND=SP),PARAMETER                  :: dTlog      = 1.25_SP ! Temperature increment (increase in dlog 
                                                                 ! Increases the swap probability and decrease 
                                                                 ! In dTlog decreases the swap probability. 
                                                                 ! This also depends on the number of parallel chains too
INTEGER(KIND=IB)                         :: NPTCHAINS1 = 1       ! # Chains T=1
INTEGER(KIND=IB)                         :: NT                   ! # Tempering levels (temperatures)
INTEGER(KIND=IB)                         :: NPTCHAINS            ! # Parallel tempering chains
INTEGER(KIND=IB)                         :: ncswap     = 0_IB    ! Temp swap accepted counter
INTEGER(KIND=IB)                         :: ncswapprop = 0_IB    ! Temp swap proposed counter
INTEGER(KIND=IB),ALLOCATABLE,DIMENSION(:):: NCHAINT
REAL(KIND=SP),   ALLOCATABLE,DIMENSION(:):: beta_pt              ! Temperature array parallel tempering
INTEGER(KIND=IB)                         :: ibirth = 0,ideath = 0

!-----------------------------------------------------------------------------
!For reference Earth model
!-----------------------------------------------------------------------------
INTEGER(KIND=IB)                         :: nreflyr
INTEGER(KIND=IB), PARAMETER              :: nlaymax = 100
REAL(KIND=SP), DIMENSION(nlaymax)        :: vsrf,vpvsrf,dthrf 
REAL(KIND=SP), ALLOCATABLE,DIMENSION(:,:):: vel_ak135

!==============================================================================
!!
!!  Sampling specific parameters
!!
!==============================================================================
INTEGER(KIND=IB)                        :: NFPMX       = (NLMX * NPL) + NPL-1
INTEGER(KIND=IB)                        :: NSDEVM      = ((NLMX*NPL)+NPL-1)*NLMX
INTEGER(KIND=IB)                        :: ioutside    = 0
INTEGER(KIND=IB)                        :: i_bd                   ! Birth-Death track (0=MCMC, 1=birth, 2=death)
INTEGER(KIND=IB)                        :: i_sdpert    = 0        ! if sigma is perturbed, don't compute forward model
INTEGER(KIND=IB)                        :: ibuf = 1, ibuf2 = 1
!==============================================================================   
!!
!!  Convergence parameters
!!
!==============================================================================

INTEGER(KIND=IB)                        :: iconv       = 0        ! Convergence switch slaves
INTEGER(KIND=IB)                        :: iconv2      = 0        ! Convergence switch master
INTEGER(KIND=IB)                        :: iconv3      = 0        ! Convergence switch master
INTEGER(KIND=IB)                        :: iarfail_pv  = 0        ! Tracks failure of AR model when predicted AR series too large
INTEGER(KIND=IB)                        :: iarfail_gv  = 0
!==============================================================================
!!
!! RJMCMC parameters
!!
!==============================================================================

INTEGER(KIND=IB),PARAMETER              :: NCHAIN     = 1E9_IB     ! Ignore (# iterations (max # MCMC steps))
INTEGER(KIND=IB),PARAMETER              :: ICHAINTHIN = 1E0_IB     ! Chain thinning interval
INTEGER(KIND=IB),PARAMETER              :: NKEEP      = 1E1_IB     ! Number models to keep before writing
INTEGER(KIND=IB),PARAMETER              :: NDM        = 100        ! No. steps in lin rot est
!==============================================================================
!!
!! Annealing burn-in parameters (sets beta schedule)
!!
!==============================================================================

REAL(KIND=SP),ALLOCATABLE,DIMENSION(:)  :: beta                    ! Inverse Temprerature
INTEGER(KIND=IB),PARAMETER              :: NTEMP = 2e3             ! Number of beta values
!==============================================================================
!!
!!  Structures for objects and data 
!!
!==============================================================================

INTEGER(KIND=IB)                        :: imcmc1 = 1              ! Counter for models at T=1 (needs to survive checkpointing!)
INTEGER(KIND=IB)                        :: imcmc2 = 1              ! Counter for mcmc steps to scale diminishing adaptation (needs to survive checkpointing!)
INTEGER                                 :: NFIELD = 31             ! The number of fields in objstruc
INTEGER                                 :: objtype1                ! Name of objtype for MPI sending
INTEGER                                 :: objtype2                ! Name of objtype for MPI sending
INTEGER                                 :: objtype3                ! Name of objtype for MPI sending
TYPE :: objstruc
  INTEGER                                  :: k                  ! Layer dimension (number of layers)
  INTEGER                                  :: NFP                ! Number forward parameters
  INTEGER                                  :: ireject_bd = 0     ! Reject birth/death
  INTEGER                                  :: iaccept_bd = 0     ! Accept birth/death
  INTEGER                                  :: iaccept    = 0     ! MHS accept
  INTEGER                                  :: ireject    = 0     ! MHS reject
  INTEGER                                  :: NDAT_pv = 0        ! No. of data points for rf
  INTEGER                                  :: NDAT_gv = 0        ! No. of data points for rf
  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:)   :: par                ! Forward parameters for layers above half space
  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:)   :: z                  ! Forward parameters (interface positions)
  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:)   :: h                  ! Depth parameters
  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:)   :: arpar_pv           ! AR for PV forward parameters
  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:)   :: arpar_gv           ! AR for GV forward parameters
  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:)   :: misc               ! Miscellaneous parameter (for e.g., gauss filter width)
  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:)   :: sdpar_pv           ! Std dev for PV forward parameters
  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:)   :: sdpar_gv           ! Std dev for GV forward parameters
  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:)   :: sdave              ! Std dev of 1000/.. of sdpar
  REAL(KIND=SP)                            :: beta               ! Tempered level
  REAL(KIND=SP)                            :: logL               ! Log likelihood
  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:,:) :: Dobs_pv            ! Observed phase velocity
  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:,:) :: Dobs_gv            ! Observed group velocity
  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:,:) :: t_pv               ! Observed/samples times phase velocity
  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:,:) :: t_gv               ! Observed/samples times group velocity
  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:,:) :: Drep_pv            ! Replica data for phase velocity
  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:,:) :: Drep_gv            ! Replica data for group velocity
  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:,:) :: res_pv             ! Data residuals for phase velocity
  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:,:) :: res_gv             ! Data residuals for group velocity
  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:,:) :: Dar_pv             ! Autoregressive model replica for phase velocity
  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:,:) :: Dar_gv             ! Autoregressive model replica for group velocity
  REAL(KIND=SP)                            :: tcmp               ! 
  REAL(KIND=SP)                            :: acceptance = 0._SP ! Accepance rate for proposal  
END TYPE objstruc

!==============================================================================
!!
!!  Global variables
!!
!==============================================================================

REAL(KIND=SP),DIMENSION(:,:),ALLOCATABLE   :: sample             ! Posterior sample
REAL(KIND=SP),DIMENSION(:),ALLOCATABLE     :: tmpmap             ! Temporary for reading map
!==============================================================================
!!
!! Temporary data arrays
!!
!==============================================================================
REAL(KIND=SP),DIMENSION(NWF,NPTS)          :: dobstmp_pv = 0._SP ! Observed temporary data for phase vel
REAL(KIND=SP),DIMENSION(NWF,NPTS)          :: dobstmp_gv = 0._SP ! Observed temporary data for group vel
REAL(KIND=SP),DIMENSION(NWF,NPTS)          :: ttmp_pv    = 0._SP ! Observed time points for rf
REAL(KIND=SP),DIMENSION(NWF,NPTS)          :: ttmp_gv    = 0._SP ! Observed time points for rf
INTEGER(KIND=IB),DIMENSION(:),ALLOCATABLE  :: icount
!==============================================================================
!!  MPI global variables
!!
!==============================================================================
INTEGER(KIND=IB)                           :: rank,NTHREAD,ncount,ncount2,ierr
INTEGER(KIND=IB), PARAMETER                :: src = 0_IB
INTEGER                                    :: to,from,tag,COMM
INTEGER                                    :: stat(MPI_STATUS_SIZE)
REAL(KIND=RP)                              :: tsave1, tsave2     ! Overall time 
INTERFACE
   FUNCTION RANDPERM(num)
   USE data_type, ONLY : IB
   IMPLICIT NONE
   INTEGER(KIND=IB), INTENT(IN)            :: num
   INTEGER(KIND=IB), DIMENSION(num)        :: RANDPERM
   END FUNCTION RANDPERM
END INTERFACE

END MODULE RJMCMC_COM

!=======================================================================
!main program for rjmcmc sampling using metropolis hasting


PROGRAM  RJMCMC_PLANE

!=======================================================================
USE MPI
USE RJMCMC_COM
USE NR
IMPLICIT NONE
INTRINSIC RANDOM_NUMBER, RANDOM_SEED

INTEGER(KIND=IB)                      :: ipar,imcmc,ithin,isource1,isource2,ikeep,it2,it,ic
INTEGER(KIND=IB)                      :: i1,i2,idat,idat1,idat2, NN, I

TYPE (objstruc),DIMENSION(2)          :: objm     ! Objects on master node
TYPE (objstruc)                       :: obj      ! Object on slave node
TYPE (objstruc)                       :: objnew   ! Object on slave node
REAL(KIND=RP)                         :: ran_uni  ! Uniform random number
REAL(KIND=RP)                         :: ran_nor  ! Normal random number
REAL(KIND=RP),DIMENSION(:),ALLOCATABLE:: tmpvoro
!==============================================================================
!!
!!     MPI stuff:
!!
!==============================================================================
!INTEGER(KIND=IB)                              :: iseedsize
!INTEGER(KIND=IB), DIMENSION(:),   ALLOCATABLE :: iseed
!INTEGER(KIND=IB), DIMENSION(:,:), ALLOCATABLE :: iseeds
!REAL(KIND=SP),    DIMENSION(:,:), ALLOCATABLE :: rseeds
REAL(KIND=RP)                                 :: tstart, tend              ! Overall time 
REAL(KIND=RP)                                 :: tstart2, tend2            ! Time for one forward model computation
!REAL(KIND=RP)                                 :: tstartsnd, tendsnd        ! Communication time
!REAL(KIND=RP)                                 :: tstartcmp, tendcmp, tcmp  ! Forward computation time

CALL MPI_INIT( ierr )                               !INITIALIZATION  OF MPI
CALL MPI_COMM_RANK( MPI_COMM_WORLD, rank, ierr )    !DETERMINES THE RANK OF CALLING PROCESS IN THE COMMUNICATOR
CALL MPI_COMM_SIZE( MPI_COMM_WORLD, NTHREAD, ierr ) !DETERMINES THE SIZE OF THE GROUP ASSOCIATED WITH A COMMUNICATOR
                                                    !I.E. TOTAL NO. OF CORES USED 
OPEN(UNIT=20,FILE=filebasefile,FORM='formatted',STATUS='OLD',ACTION='READ')
READ(20,*) filebaselen !LENGTH OF noisy_data
READ(20,*) filebase    !name of the file, noisy_data
CLOSE(20)
!REPRESENTATION OF ALL INPUT AND OUTPUT FILES
infile_pv      = filebase(1:filebaselen) // '_pv.txt'
infile_gv      = filebase(1:filebaselen) // '_gv.txt'
logfile        = filebase(1:filebaselen) // '_RJMH.log'
seedfile       = filebase(1:filebaselen) // '_seeds.log'
mapfile        = filebase(1:filebaselen) // '_map.dat'
repfile_pv     = filebase(1:filebaselen) // '_rep_pv.dat'
repfile_gv     = filebase(1:filebaselen) // '_rep_gv.dat'
obsfile_pv     = filebase(1:filebaselen) // '_obs_pv.dat'
obsfile_gv     = filebase(1:filebaselen) // '_obs_gv.dat'
covfile        = filebase(1:filebaselen) // '_cov.txt'
sdfile         = filebase(1:filebaselen) // '_sigma.txt'
samplefile     = filebase(1:filebaselen) // '_sample.txt'
stepsizefile   = filebase(1:filebaselen) // '_stepsize.dat'
refmodel       = filebase(1:filebaselen) // '_ref_ak135.model'
!DISPLAYED IN THE SCREEN WHEN RANK == 0 (SRC)
305 FORMAT(I8)
304 FORMAT(A100)

OPEN(ulog,FILE=logfile,STATUS='REPLACE')
IF(rank == src)THEN
    WRITE(ulog,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  '
    WRITE(ulog,*) '  ~~~~~~  '
    WRITE(ulog,*) '  ~~~           Reversible Jump MCMC Sampling ~~~  '
    WRITE(ulog,*) '  ~~~~~~  '
    WRITE(ulog,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  '
    WRITE(ulog,*) '...running on ',NTHREAD,' cores'
ENDIF
CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
ncount = NKEEP*(NFPMX+NWF+(NARFP*NWF)+NMISC)
ncount2= 4+NFPMX+NMISC+2*NWF+2*NWF+6
ALLOCATE( sample(NKEEP,ncount2) )
ALLOCATE( tmpmap(1+NFPMX+NMISC+2*NWF+2*NWF) )
sample = 0._SP
tmpmap = 0._SP
!==============================================================================
!!
!! Population / parallel tempering:
!!
!==============================================================================
IF(NTHREAD > 1)THEN
  NPTCHAINS  = NTHREAD-1    ! Number of chains (equal to number of slaves)
  NT = NPTCHAINS-NPTCHAINS1+1
  ALLOCATE( NCHAINT(NT),beta_pt(NPTCHAINS) )
  NCHAINT = 1
  NCHAINT(1) = NPTCHAINS1

!==============================================================================
!!  Set up tempering schedule
!==============================================================================
  IF(rank == src)WRITE(ulog,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
  IF(rank == src)WRITE(ulog,*)'NT',NT
  IF(rank == src)WRITE(ulog,*)'NCHAINT',NCHAINT
  it2 = 1_IB
  DO it=1,NT
    DO ic=1,NCHAINT(it)
      beta_pt(it2) = 1._SP/dTlog**REAL(it-1_IB,SP)
      IF(rank == src)WRITE(ulog,219) 'Chain ',it2,':  beta = ',beta_pt(it2),'  T = ',1._SP/beta_pt(it2)
      it2 = it2 + 1
    ENDDO
  ENDDO
ELSE
   ALLOCATE( NCHAINT(1) )
ENDIF
IF(rank == src)WRITE(ulog,*) ''
IF(rank == src)WRITE(ulog,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
219 FORMAT(a6,i5,a10,F15.5,a6,F15.5)

IF(rank == src)WRITE(ulog,*) 'k start   = ',obj%k

IF(rank /= src)obj%beta = beta_pt(rank)
  OPEN(20,file=infile_pv,status='OLD')
  idat1      = 0
  ttmp_pv    = 0.0
  dobstmp_pv = 0.0
  DO idat = 1,NPTS
    READ(20,*, END=299) ttmp_pv(1,idat),dobstmp_pv(1,idat)
    idat1 = idat1 + 1
  ENDDO
  299   CLOSE(20)
  obj%NDAT_pv = idat1

  OPEN(20,file=infile_gv,status='OLD')
  idat1      = 0
  ttmp_gv    = 0.0
  dobstmp_gv = 0.0
  DO idat = 1,NPTS
     READ(20,*, END=399) ttmp_gv(1,idat),dobstmp_gv(1,idat)
     idat1 = idat1 + 1
  ENDDO
  399   CLOSE(20)
  obj%NDAT_gv = idat1
!======================================================================
  objm(1)%NDAT_pv  = obj%NDAT_pv
  objm(2)%NDAT_pv  = obj%NDAT_pv
  objnew%NDAT_pv   = obj%NDAT_pv

  objm(1)%NDAT_gv  = obj%NDAT_gv
  objm(2)%NDAT_gv  = obj%NDAT_gv
  objnew%NDAT_gv   = obj%NDAT_gv
!!
!!  Allocate Objects
!!
CALL ALLOC_OBJ(objm(1))
CALL ALLOC_OBJ(objm(2))
CALL ALLOC_OBJ(obj)
CALL ALLOC_OBJ(objnew)
!!
!!  Make MPI structure to match objstruc
!!
CALL MAKE_MPI_STRUC_SP(objm(1),objtype1)
CALL MAKE_MPI_STRUC_SP(objm(2),objtype2)
CALL MAKE_MPI_STRUC_SP(obj,objtype3)

!=======================================================================
!!
!!Store all data and variable in the allocated memory
!!
!=======================================================================
obj%Dobs_pv = 0.0
obj%t_pv    = 0.0
obj%Dobs_gv = 0.0
obj%t_gv    = 0.0
DO I = 1, obj%NDAT_pv
   obj%Dobs_pv(1,I) = dobstmp_pv(1,I)
   obj%t_pv(1,I)    = ttmp_pv(1,I)
ENDDO
DO I = 1, obj%NDAT_gv
   obj%Dobs_gv(1,I) = dobstmp_gv(1,I)
   obj%t_gv(1,I)    = ttmp_gv(1,I)
ENDDO

IF(rank == 0)THEN  
   WRITE(ulog,*) '=================================='
   WRITE(ulog,*) 'DATA WINDOWING and RESAMPLING DONE'
   WRITE(ulog,*) '=================================='
    !!
    !!  File to save posterior samples
    !!
    OPEN(UNIT=usample,FILE=samplefile,FORM='formatted',STATUS='REPLACE', &
    ACTION='WRITE')
    OPEN(UNIT=ustep,FILE=stepsizefile,FORM='formatted',STATUS='REPLACE', &
    ACTION='WRITE')   
ENDIF  
!=======================================================================

!-----------------------------------------------------------------------
!Open and read the ak135 1-D Earth model
!-----------------------------------------------------------------------
nreflyr = 0
OPEN(38,FILE=refmodel,STATUS='OLD')
DO I = 1,nlaymax
  READ(38,*,END=400)dthrf(I),vsrf(I)
  nreflyr = nreflyr+1
ENDDO
400  CLOSE(38)
ALLOCATE(vel_ak135(nreflyr,2))
vel_ak135 = 0.0
DO I = 1, nreflyr
   vel_ak135(I,1) = dthrf(I)
   vel_ak135(I,2) = vsrf(I)
ENDDO
!=======================================================================
!!  Prior bounds
!=======================================================================
!!           h      dvs 
minlim = (/ hmin, -1.00 /)
maxlim = (/ hmx,   1.00 /)
!----------------------------------------------------------------------

!! MAximum perturbation size
maxpert = maxlim-minlim
pertsd = maxpert/pertsdsc

IF(ICOV == 1)THEN
  !! Set prior and proposal scaling for data error standard deviations:
  maxlimsd   = 1.0000_SP  !! MAXLIM IS NOW READ FROM DATA FILE!!
  pertsdsdsc = 50._SP
  maxpertsd  = maxlimsd-minlimsd
  pertsdsd   = maxpertsd/pertsdsdsc
ENDIF
IF(IAR == 1)THEN
  !! Set prior and proposal scaling for data error standard deviations:
  pertarsdsc = 100._SP
  maxpertar  = maxlimar-minlimar
  pertarsd   = maxpertar/pertarsdsc
ENDIF

IF(IMISC == 1)THEN
   !! Set prior and proposal scaling for data error standard deviations:
   minlimmisc   =  0.1_SP
   maxlimmisc   = 10.0_SP     !! MAXLIM IS NOW READ FROM DATA FILE!!
   pertmiscsdsc = 40._SP
   maxpertmisc  = maxlimmisc-minlimmisc
   pertmiscsd   = maxpertmisc/pertmiscsdsc
ENDIF

!=======================================================================
!!  Print sampling parameters to screen for logging
!=======================================================================
IF(rank == src)THEN
         WRITE(ulog,*) 'ICOV                   :   ',ICOV  !210
         IF(ICOV == 1)THEN
            WRITE(ulog,*) '...sampling over standard deviations.' !204
         ENDIF
         WRITE(ulog,*) 'NKEEP                  :   ',NKEEP
         WRITE(ulog,*) 'ICHAINTHIN             :   ',ICHAINTHIN
         WRITE(ulog,*) 'IEXCHANGE              :   ',IEXCHANGE
         WRITE(ulog,*) 'NT                     :   ',NT
         WRITE(ulog,*) 'NPTCHAINS              :   ',NPTCHAINS
         WRITE(ulog,*) 'NCHAINT                :   ',NCHAINT
         WRITE(ulog,*) 'IAR                    :   ',IAR
         WRITE(ulog,*) 'IMISC                  :   ',IMISC
         WRITE(ulog,*) '# procs                :   ',NWF
         WRITE(ulog,*) 'Order                  :   ',NARFP
         WRITE(ulog,*) 'Number of data in PV   :   ',obj%NDAT_pv
         WRITE(ulog,*) 'Number of data in GV   :   ',obj%NDAT_gv
         WRITE(ulog,*) 'Number of Waveforms    :   ',NWF
         CALL FLUSH(ulog)
         WRITE(ulog,209) 'Sample file:             ',samplefile
         WRITE(ulog,*) ''
         WRITE(ulog,*) 'NFPMX            = ',NFPMX
         WRITE(ulog,*) 'IMAP             = ',IMAP
         WRITE(ulog,*) 'kmin             = ',kmin
         WRITE(ulog,*) 'kmax             = ',kmax
         WRITE(ulog,*) 'hmin             = ',hmin
         WRITE(ulog,*) 'hmax             = ',hmx
         WRITE(ulog,*) 'minlim:  '
         WRITE(ulog,*) minlim !201
         WRITE(ulog,*) 'maxlim:  '
         WRITE(ulog,*) maxlim !201
         WRITE(ulog,*) 'pertsdsc:'
         WRITE(ulog,*) pertsdsc !201
         WRITE(ulog,*) 'maxlim sigma:  '
         WRITE(ulog,*) maxlimsd !201
         WRITE(ulog,*) ''
         WRITE(ulog,*) 'Done reading data.'
         WRITE(ulog,*) ''
         CALL FLUSH(ulog)
         WRITE(ulog,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  '
         WRITE(ulog,*) ''
ENDIF
!!
!! Initialize random seeds on each core (Call RANDOM_SEED only once in the whole code. PARALLEL_SEED calls it)
!!
!CALL RANDOM_SEED
CALL PARALLEL_SEED()

!!--------------------
!Initialize the MCMC
!---------------------
!!
!! Read mapfile to start:
!!
OPEN(UNIT=20,FILE=mapfile,FORM='formatted',STATUS='OLD',ACTION='READ')
!READ(20,*)  !skip header information
READ(20,*) tmpmap
CLOSE(20)
obj%k = tmpmap(1)
!obj%k   = 1 !INT(tmpmap(1),IB)
obj%NFP = (obj%k*NPL)+NPL-1
obj%par = 0._SP
obj%sdpar_pv = 0._SP
obj%sdpar_gv = 0._SP
obj%misc= 0._SP
obj%h   = 0._SP
obj%sdpar_pv = 0._SP
obj%sdpar_gv = 0._SP
obj%arpar_pv = 0._SP
obj%arpar_gv = 0._SP
obj%par(1:obj%NFP)= tmpmap(2:obj%NFP+1)
IF(IMISC == 1)obj%misc(1:NMISC) = tmpmap(1+NFPMX+NMISC:1+NFPMX+NMISC)
IF(ICOV == 1) THEN
 IF(IPV==1) obj%sdpar_pv = tmpmap(1+NFPMX+NMISC+1:1+NFPMX+NMISC+NWF)
 IF(IGV==1) obj%sdpar_gv = tmpmap(1+NFPMX+NMISC+NWF+1:1+NFPMX+NMISC+NWF+NWF)
ENDIF
IF(IAR == 1)THEN
 IF(IPV==1) obj%arpar_pv = tmpmap(1+NFPMX+NMISC+2*NWF+1:1+NFPMX+NMISC+2*NWF+NWF)
 IF(IGV==1) obj%arpar_gv = tmpmap(1+NFPMX+NMISC+2*NWF+NWF+1:1+NFPMX+NMISC+2*NWF+NWF+NWF)
ENDIF
CALL CALCH(obj)
ALLOCATE( icount(NTHREAD) )
icount = 0
tstart = MPI_WTIME()
!=======================================================================
!Inner Do loop is for sweeping
!=======================================================================
IF(IMAP == 1)THEN
   IF(rank == src)WRITE(ulog,*) 'IMAP activated, exiting after computing replica for MAP.'
   CALL CHECKBOUNDS(obj)
   IF(ioutside == 1) WRITE(ulog,*)'outside the parameter bounds',ioutside
   tstart2 = MPI_WTIME()
   IF(ISMPPRIOR == 0) CALL LOGLHOOD(obj)
   IF(ISMPPRIOR == 1) CALL LOGLHOOD2(obj)
!=======================================================================
!Here this do loop is for sweeping
!=======================================================================
   !DO it = 1,200
   !    obj%hb = minlim(1)+REAL(it)*(hmx-minlim(1))/500    
   !   IF(ISMPPRIOR == 0) CALL LOGLHOOD(obj)
   !   IF(ISMPPRIOR == 1) CALL LOGLHOOD2(obj)
   !   tend2 = MPI_WTIME()
   !   WRITE(77,*)obj%hb,obj%logL
   !   WRITE(6,*)it,obj%hb,obj%logL
   !ENDDO
   !STOP
   tend2 = MPI_WTIME()
   IF(rank == src)CALL SAVEREPLICA(obj,repfile_pv,repfile_gv,obsfile_pv,obsfile_gv)
   IF(rank == src)WRITE(ulog,*) 'time for one logL evaluation:',real(tend2-tstart2)

   STOP
ELSE !IMAP ELSE

!=======================================================================
!IF IMAP .NE. 1 THEN DISPLAY THE STARTING MODEL PARAMETERS AND COMPUTE THE 
!LOG LIKELIHOOD AND CHECK THE BOUNDARY OF STARTING MODEL AS IT IS STARTED
!RANDOMLY !PRINTS ONLY BY ZERO RANKED CPU
!=======================================================================
     tstart2 = MPI_WTIME()
     IF(ISMPPRIOR == 0) CALL LOGLHOOD(obj)
     IF(ISMPPRIOR == 1) CALL LOGLHOOD2(obj)
     IF(rank == src)WRITE(ulog,*) 'Starting model:'
     IF(rank == src)CALL PRINTPAR(obj)
     IF(rank == src)WRITE(ulog,*) 'logL = ',obj%logL
     tend2 = MPI_WTIME()
     IF(rank == src)WRITE(ulog,*) 'time for one logL evaluation: ',tend2-tstart2
     CALL CHECKBOUNDS(obj)
     IF(ioutside == 1)WRITE(*,*)'FAILED on bounds'
     IF(ioutside == 1)STOP
ENDIF !END IMAP
CALL FLUSH(ulog)
icount(rank+1) = icount(rank+1) + 1   !COUNTING 
!=======================================================================
!
!          ************ RJMCMC Sampling ************
!=======================================================================

!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!
!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!
!!
!!     MASTER PART
!!
!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!
!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!
IF(rank==src)THEN
  ncswap = 0
  ncswapprop = 0
  ikeep = 1
  tsave1 = MPI_WTIME()
  DO imcmc = 1,NCHAIN
    !!
    !!  Receiving samples from any two slave (no particular order -> auto load
    !!  balancing) and propose swap:
    !!
    DO ithin = 1,ICHAINTHIN
      CALL MPI_RECV(objm(1), 1, objtype1, MPI_ANY_SOURCE,MPI_ANY_TAG, MPI_COMM_WORLD, stat, ierr )
      isource1 = stat(MPI_SOURCE)  !! This saves the slave id for the following communication
      tstart = MPI_WTIME()
      CALL MPI_RECV(objm(2), 1, objtype2, MPI_ANY_SOURCE,MPI_ANY_TAG, MPI_COMM_WORLD, stat, ierr )
      isource2 = stat(MPI_SOURCE)  !! This saves the slave id for the following communication
      IF(IEXCHANGE == 1)CALL TEMPSWP_MH(objm(1),objm(2))
      CALL MPI_SEND(objm(1), 1,objtype1, isource1, rank, MPI_COMM_WORLD, ierr)
      CALL MPI_SEND(objm(2), 1,objtype2, isource2, rank, MPI_COMM_WORLD, ierr)
      tend = MPI_WTIME()
    ENDDO
    !!Checkpoint every .... seconds:
    CALL SAVESAMPLE(objm,ikeep,isource1,isource2,REAL(tend-tstart,SP))
  ENDDO
!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
!!
!!    SLAVE PART
!!
!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!
!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!
ELSE
  DO imcmc = 1,NCHAIN
    !! Apply chain thinning
    !DO ithin = 1,5
       CALL EXPLORE_MH(obj,objnew,obj%beta)
    !ENDDO
    tend = MPI_WTIME()
    obj%tcmp = REAL(tend-tstart,RP)
    CALL MPI_SEND(obj, 1,objtype3, src, rank, MPI_COMM_WORLD, ierr)
    CALL MPI_RECV(obj, 1,objtype3, src, MPI_ANY_TAG, MPI_COMM_WORLD, stat, ierr )
    tstart = MPI_WTIME()
    obj%acceptance =REAL(obj%iaccept,SP)/REAL(obj%iaccept+obj%ireject,SP)
  ENDDO  !END imcmc
ENDIF !! MPI ENDIF

201   FORMAT(100F12.4)
202   FORMAT(I8,3F16.6,I8,F16.6,I8,1F13.3,I8,1F13.3)
203   FORMAT(A124)
204   FORMAT(A56)
205   FORMAT(10F10.4)
209   FORMAT(A26,A40)
210   FORMAT(A26,20I10)
211   FORMAT(20000ES12.4)
212   FORMAT(I4,A8,F12.4,A5,I3)
213   FORMAT(I4,I6,A8,F12.4,A5,I3,A14,F12.4,A11,I3)
214   FORMAT(a21,18I6)
215   FORMAT(a21,F8.4)
216   FORMAT(a21,F8.4)
217   FORMAT(a21,F8.4)
218   FORMAT(4F15.4)
619   FORMAT(6F15.6)
      CALL MPI_FINALIZE( ierr ) 

END PROGRAM RJMCMC_PLANE

!=======================================================================

SUBROUTINE LOGLHOOD(obj)
!=======================================================================
USE RJMCMC_COM
USE MPI
IMPLICIT NONE
INTEGER(KIND=IB)                  :: i,i1,irf
TYPE (objstruc)                   :: obj,objtmp
REAL(KIND=SP), DIMENSION(NPTS)    :: yrf
REAL(KIND=SP), DIMENSION(NWF)     :: Etmp_pv, Etmp_gv
!----------------------------------------------------------------------
!For surface wave dispersion
!---------------------------------------------------------------------
REAL(KIND=SP),DIMENSION(obj%NDAT_pv):: yrf_pv
REAL(KIND=SP),DIMENSION(obj%NDAT_gv):: yrf_gv
REAL(KIND=SP),DIMENSION(obj%k+1)    :: thick,vp,vs,den
REAL(KIND=SP),DIMENSION(obj%k+1)    :: thick1,vs1,vp1,den1
REAL(KIND=SP)                       :: tmppar(obj%NFP)

CALL CALCH(obj)
tmppar = 0.0

!---------------------------------------------------------------------
!Prepare model for forward calculation applying different algorithms
!---------------------------------------------------------------------
!---------------------------------------------------------------------
!For half space
!---------------------------------------------------------------------
CALL ADD_REF_MDL(obj,tmppar)
vs    = 0.0
thick = 0.0
vp    = 0.0
den   = 0.0
IF(obj%k == 0)THEN
  thick(obj%k+1) = 0.0
  vs(obj%k+1)    = tmppar(1)
  vp(obj%k+1)    = vpvs*vs(obj%k+1)
  den(obj%k+1)   = (vp(obj%k+1)+1.87)/3.05
ELSE
  thick(obj%k+1)= 0.0
  vs(obj%k+1)   = tmppar(obj%k*NPL+1)
  vp(obj%k+1)   = vpvs*vs(obj%k+1)
  den(obj%k+1)  = (vp(obj%k+1)+1.87)/3.05
  DO i1 = 1, obj%k
    thick(i1)   = obj%h(i1)
    vs(i1)      = tmppar(NPL*i1)
    vp(i1)      = vpvs*vs(i1)
    den(i1)     = (vp(i1)+1.87)/3.05
  ENDDO
ENDIF
thick1 = thick
vp1 = vp
vs1 = vs
den1= den
!do i1 = 1, obj%k+1
!  print*,thick(i1),vp(i1),vs(i1),den(i1)
!enddo
!-----------------------------------------------------------------------
!For phase velocity calculation using raydsp code from Doornbos
!-----------------------------------------------------------------------
IF(IPV==1)THEN
 DO irf=1,NWF
   yrf_pv             = 0._SP
   obj%Drep_pv(irf,:) = 0._SP
   obj%res_pv(irf,:)  = 0._SP
  IF(IHRMN == 1)THEN
   CALL compute_phvel_hrmn(obj%k+1,thick,vp,vs,den,obj%t_pv,obj%NDAT_pv,yrf_pv)
  ELSE
   CALL dispersion_pv(obj%k+1,den,vp,vs,thick,yrf_pv,obj%t_pv,obj%NDAT_pv)
  ENDIF
!   CALL COMPUTE_SWD(thick,vp,vs,den,obj%t_pv,obj%k+1,obj%NDAT_pv,yrf_pv,yrf_gv)
   obj%Drep_pv(irf,1:obj%NDAT_pv) = yrf_pv(1:obj%NDAT_pv)
   obj%res_pv(irf,1:obj%NDAT_pv) = obj%Dobs_pv(irf,1:obj%NDAT_pv)-obj%Drep_pv(irf,1:obj%NDAT_pv)
!======================================================================
   iarfail_pv = 0
   IF(IAR == 1)THEN
     obj%Dar_pv  = 0._SP
!=======================================================================
   !!  Compute autoregressive model
!=======================================================================
      CALL ARFB_PV(obj,1,obj%NDAT_pv)
!=======================================================================
! Recompute replica data as ith autoregressive model
!=======================================================================
     obj%res_pv(irf,1:obj%NDAT_pv)=obj%res_pv(irf,1:obj%NDAT_pv)-obj%Dar_pv(irf,1:obj%NDAT_pv)

!=======================================================================
!Check if predicted AR model data are outside max allowed bounds
!============================================================  
     CALL CHECKBOUNDS_ARMX(obj)
   ENDIF
 ENDDO
ENDIF
!-----------------------------------------------------------------------
!For group velocity surface wave dispersion using raydisp
!-----------------------------------------------------------------------
IF(IGV == 1)THEN
 DO irf=1,NWF
  obj%Drep_gv(irf,:) = 0._SP
  obj%res_gv(irf,:)  = 0._SP
  CALL dispersion_gv(obj%k+1,den1,vp1,vs1,thick1,yrf_gv,obj%t_gv,obj%NDAT_gv)
!  CALL COMPUTE_SWD(thick1,vp1,vs1,den1,obj%t_pv,obj%k+1,obj%NDAT_pv,yrf_pv,yrf_gv)
  obj%Drep_gv(irf,1:obj%NDAT_gv) = yrf_gv(1:obj%NDAT_gv)
  obj%res_gv(irf,1:obj%NDAT_gv) = obj%Dobs_gv(irf,1:obj%NDAT_gv)-obj%Drep_gv(irf,1:obj%NDAT_gv)
!======================================================================
  iarfail_gv = 0
  IF(IAR == 1)THEN
    obj%Dar_gv  = 0._SP
!=======================================================================
   !!  Compute autoregressive model
!=======================================================================
    CALL ARFB_GV(obj,1,obj%NDAT_gv)
!=======================================================================
! Recompute replica data as ith autoregressive model
!=======================================================================
    obj%res_gv(irf,1:obj%NDAT_gv)=obj%res_gv(irf,1:obj%NDAT_gv)-obj%Dar_gv(irf,1:obj%NDAT_gv)

!=======================================================================
!Check if predicted AR model data are outside max allowed bounds
!============================================================  
    CALL CHECKBOUNDS_ARMX(obj)
  ENDIF
 ENDDO
ENDIF
!=======================================================================
!!  Compute log likelihood for Phase velocity(PV) and, Group Velocity (GV)
!=======================================================================
IF(IPV == 1)THEN
 IF(iarfail_pv == 0)THEN
  IF (ICOV == 0) THEN
      !!
      !! Implicitly sample over sigma
      !!
    Etmp_pv = 0._SP
    DO irf = 1,NWF
        Etmp_pv(irf) = REAL(obj%NDAT_pv,SP)/2._SP*LOG(SUM(obj%res_pv(irf,1:obj%NDAT_pv)*obj%res_pv(irf,1:obj%NDAT_pv)))
        obj%sdpar_pv(irf) = sqrt((1._SP/REAL(obj%NDAT_pv))*SUM(obj%res_pv(irf,1:obj%NDAT_pv)**2._SP))
    ENDDO
    Etmp_pv = -Etmp_pv
  ELSEIF(ICOV == 1)THEN
      !!
      !! Sample over sigma (one for all freqs)
      !!
    Etmp_pv = 0._SP
    DO irf = 1,NWF
      Etmp_pv(irf)=LOG(1._SP/(2._SP*PI2))*(REAL(obj%NDAT_pv,SP)/2._SP)-REAL(obj%NDAT_pv)*LOG(obj%sdpar_pv(irf)) &
                  -SUM(obj%res_pv(irf,1:obj%NDAT_pv)*obj%res_pv(irf,1:obj%NDAT_pv))/(2._SP*obj%sdpar_pv(irf)**2._SP)
    ENDDO
  ENDIF
 ELSE
  Etmp_pv   = -HUGE(1._SP) 
  iarfail_pv = 0
 ENDIF
ELSE
 Etmp_pv = 0._SP
ENDIF

IF(IGV == 1)THEN
 IF(iarfail_gv == 0)THEN
  IF (ICOV == 0) THEN
      !!
      !! Implicitly sample over sigma
      !!
    Etmp_gv = 0._SP
    DO irf = 1,NWF
        Etmp_gv(irf) = REAL(obj%NDAT_gv,SP)/2._SP*LOG(SUM(obj%res_gv(irf,1:obj%NDAT_gv)*obj%res_gv(irf,1:obj%NDAT_gv)))
        obj%sdpar_gv(irf) = sqrt((1._SP/REAL(obj%NDAT_gv))*SUM(obj%res_gv(irf,1:obj%NDAT_gv)**2._SP))
    ENDDO
    Etmp_gv = -Etmp_gv
  ELSEIF(ICOV == 1)THEN
      !!
      !! Sample over sigma (one for all freqs)
      !!
    Etmp_gv = 0._SP
    DO irf = 1,NWF
      Etmp_gv(irf)=LOG(1._SP/(2._SP*PI2))*(REAL(obj%NDAT_gv,SP)/2._SP)-REAL(obj%NDAT_gv)*LOG(obj%sdpar_gv(irf)) &
                  -SUM(obj%res_gv(irf,1:obj%NDAT_gv)*obj%res_gv(irf,1:obj%NDAT_gv))/(2._SP*obj%sdpar_gv(irf)**2._SP)
    ENDDO
  ENDIF
 ELSE
  Etmp_gv   = -HUGE(1._SP) 
  iarfail_gv= 0
 ENDIF
ELSE
 Etmp_gv = 0._SP
ENDIF
obj%logL = SUM(Etmp_pv)+SUM(Etmp_gv)
RETURN
207   FORMAT(5000ES18.8)
END SUBROUTINE LOGLHOOD

SUBROUTINE LOGLHOOD2(obj)
!=======================================================================
USE RJMCMC_COM
USE MPI
IMPLICIT NONE
TYPE (objstruc)  :: obj
obj%logL = 1._SP
RETURN
207   FORMAT(5000ES18.8)
!==============================================================================
END SUBROUTINE LOGLHOOD2

!==============================================================================
SUBROUTINE TEMPSWP_MH(obj1,obj2)
!==============================================================================
USE DATA_TYPE
USE RJMCMC_COM
IMPLICIT NONE
TYPE(objstruc)              :: obj1
TYPE(objstruc)              :: obj2
TYPE(objstruc)              :: objtmp1,objtmp2
REAL(KIND=SP)               :: ran_uni1
REAL(KIND=SP)               :: logratio
REAL(KIND=SP)               :: betaratio
betaratio = obj2%beta-obj1%beta
logratio  = betaratio*(obj1%logL-obj2%logL)
CALL RANDOM_NUMBER(ran_uni1)
IF(ran_uni1 <= EXP(logratio))THEN
  objtmp1       = obj1
  objtmp2       = obj2
  obj1   = objtmp2
  obj2   = objtmp1
  !! Temperature does not swap
  obj1%beta = objtmp1%beta
  obj2%beta = objtmp2%beta
  !! Acceptance death, birth, proposal does not swap
  obj1%iaccept_bd = objtmp1%iaccept_bd
  obj2%iaccept_bd = objtmp2%iaccept_bd
  obj1%ireject_bd = objtmp1%ireject_bd
  obj2%ireject_bd = objtmp2%ireject_bd
  obj1%ireject    = objtmp1%ireject
  obj2%ireject    = objtmp2%ireject
  obj1%iaccept    = objtmp1%iaccept
  obj2%iaccept    = objtmp2%iaccept
  obj1%acceptance = objtmp1%acceptance
  obj2%acceptance = objtmp2%acceptance
  ncswap    = ncswap+1_IB
ENDIF
ncswapprop    = ncswapprop+1_IB

RETURN
END SUBROUTINE TEMPSWP_MH
!=======================================================================


SUBROUTINE EXPLORE_MH(obj,objnew1,beta_mh)
!==============================================================================
USE DATA_TYPE
USE RJMCMC_COM
IMPLICIT NONE
INTEGER(KIND=IB)                            :: iwhich,ipar
TYPE(objstruc)                              :: obj,objnew1
REAL(KIND=SP)                               :: logPLratio,logPratio
REAL(KIND=SP)                               :: ran_uni,ran_uni_BD, ran_unik,ran_uni_ar
REAL(KIND=SP)                               :: beta_mh
INTEGER(KIND=IB),DIMENSION(NFPMX)           :: idxrand
INTEGER(KIND=IB),DIMENSION((kmax+1)*(NPL-1)):: idxvpvsr
REAL(KIND=SP),DIMENSION(obj%k)              :: ztmp
REAL(KIND=SP)                               :: logarp
INTEGER(KIND=IB)                            :: arptype
objnew1 = obj
!! Draw uniform Birth-Death probability
CALL RANDOM_NUMBER(ran_uni_BD)
i_bd = 0
!! Do BIRTH-DEATH MCMC with 0.5 probability
IF(kmin /= kmax)THEN
  !! Perturbing k:
  CALL RANDOM_NUMBER(ran_unik)
  IF(obj%k == kmax)THEN
    !! If k == kmax, no birth allowed, so 1/3 death, 2/3 stay local
    IF(ran_unik<=0.3333_SP)i_bd = 2
  ELSEIF(obj%k == kmin)THEN
    !! If k == kmin, no death allowed, so 1/3 birth, 2/3 stay local
    IF(ran_unik<=0.3333_SP)i_bd = 1
  ELSE
    !! 1/3 birth, 1/3 death, 1/3 stay local
    IF((ran_unik <= 0.3333_SP))i_bd = 1
    IF((ran_unik > 0.6666_SP))i_bd = 2
  ENDIF
  IF(i_bd == 1)THEN
       CALL BIRTH(obj,objnew1)
  ELSEIF(i_bd == 2)THEN 
       CALL DEATH(obj,objnew1)
  ENDIF
ENDIF
!==================================================================
!!
!! Do Metropolis-Hastings
!!
!==================================================================
IF(obj%k /= objnew1%k)THEN !if k changes
    !!
    !! If k changed, check BD acceptance
    !!
    CALL CHECKBOUNDS(objnew1)
    IF(ioutside == 0)THEN
      IF(ISMPPRIOR == 0) CALL LOGLHOOD(objnew1)
      IF(ISMPPRIOR == 1) CALL LOGLHOOD2(objnew1)
      IF(i_bd == 1)THEN 
        !! BIRTH  ACCEPTANCE RATIO
         logPLratio =  (objnew1%logL - obj%logL)*beta_mh
      ELSEIF(i_bd == 2)THEN
        !! DEATH ACCEPTANCE RATIO
         logPLratio =  (objnew1%logL - obj%logL)*beta_mh
      ENDIF
      CALL RANDOM_NUMBER(ran_uni)
      IF(ran_uni >= EXP(logPLratio))THEN
        objnew1 = obj
        obj%ireject_bd = obj%ireject_bd + 1
      ELSE
        obj = objnew1
        obj%iaccept_bd = obj%iaccept_bd + 1
      ENDIF
    ELSE
      objnew1 = obj
      obj%ireject_bd = obj%ireject_bd + 1
      ioutside = 0
    ENDIF
ELSE  !if k does not change
!!------------------------------------------------------
!! If k did not change, carry out MH sweep
!!------------------------------------------------------
   idxrand = 0
   idxrand(1:obj%NFP) = RANDPERM(obj%NFP)
!!-------------------------------------------------------
!! Do Metropolis-Hastings update on vp, vs and rho on the 
!!-------------------------------------------------------
   DO ipar = 1,obj%NFP
        iwhich = idxrand(ipar)
        CALL PROPOSAL(obj,objnew1,iwhich,1._SP)
        CALL CHECKBOUNDS(objnew1)
        IF(ioutside == 0)THEN
          !CALL LOGLHOOD(objnew1)
          IF(ISMPPRIOR == 0) CALL LOGLHOOD(objnew1)
          IF(ISMPPRIOR == 1) CALL LOGLHOOD2(objnew1)
          logPLratio = (objnew1%logL - obj%logL)*beta_mh       !likelihood ratio with parallel tempering
          CALL RANDOM_NUMBER(ran_uni)
          IF(ran_uni >= EXP(logPLratio))THEN
            objnew1 = obj
            obj%ireject = obj%ireject + 1
          ELSE
            obj = objnew1            
            obj%iaccept = obj%iaccept + 1
          ENDIF
        ELSE !! outside before delayed rejection
          objnew1 = obj
          obj%ireject = obj%ireject + 1
          ioutside = 0
        ENDIF
   ENDDO    !ipar enddo

!!----------------------------------------------------
!! Do Metropolis-Hastings on misc parameters
!!----------------------------------------------------
   IF(IMISC == 1)THEN
   !! Perturb std devs with .25 probability
       CALL RANDOM_NUMBER(ran_uni_ar)
       IF(ran_uni_ar<=0.5_SP)THEN
         DO ipar = 1,NWF
            CALL PROPOSAL_MISC(obj,objnew1,ipar)
!           i_sdpert = 1
            IF(ioutside == 0)THEN
               !CALL LOGLHOOD(objnew1)
               IF(ISMPPRIOR == 0) CALL LOGLHOOD(objnew1)
               IF(ISMPPRIOR == 1) CALL LOGLHOOD2(objnew1)
               logPLratio = (objnew1%logL - obj%logL)*beta_mh
               CALL RANDOM_NUMBER(ran_uni)
               IF(ran_uni >= EXP(logPLratio))THEN
                  objnew1 = obj
                  obj%ireject = obj%ireject + 1
               ELSE
                   obj = objnew1
                   obj%iaccept = obj%iaccept + 1
               ENDIF
            ELSE    !ioutside else
               objnew1 = obj
               obj%ireject = obj%ireject + 1
               ioutside = 0
            ENDIF   !ioutside endif
               i_sdpert = 0
         ENDDO      !ipar enddo
       ENDIF        !ran_un_ar endif
   ENDIF ! IMISC if

!!
!! Do Metropolis-Hastings on data-error standard deviations
!!
   IF(ICOV == 1)THEN
       !! Perturb std devs with .25 probability
     IF(IPV == 1)THEN 
        CALL RANDOM_NUMBER(ran_uni_ar)
        IF(ran_uni_ar<=0.25_SP)THEN
           DO ipar = 1,NWF
              CALL PROPOSAL_SDPV(obj,objnew1,ipar)
!              i_sdpert = 1
              IF(ioutside == 0)THEN
               !CALL LOGLHOOD(objnew1)
                 IF(ISMPPRIOR == 0) CALL LOGLHOOD(objnew1)
                 IF(ISMPPRIOR == 1) CALL LOGLHOOD2(objnew1)
                 logPLratio = (objnew1%logL - obj%logL)*beta_mh
                 CALL RANDOM_NUMBER(ran_uni)
                 IF(ran_uni >= EXP(logPLratio))THEN
                    objnew1 = obj
                    obj%ireject = obj%ireject + 1
                 ELSE
                    obj = objnew1
                    obj%iaccept = obj%iaccept + 1
                 ENDIF
              ELSE
                 objnew1 = obj
                 obj%ireject = obj%ireject + 1
                 ioutside = 0
              ENDIF
                 i_sdpert = 0
           ENDDO
        ENDIF
      ENDIF
      IF(IGV == 1)THEN
       !! Perturb std devs with .25 probability
        CALL RANDOM_NUMBER(ran_uni_ar)
        IF(ran_uni_ar<=0.25_SP)THEN
           DO ipar = 1,NWF
              CALL PROPOSAL_SDGV(obj,objnew1,ipar)
!              i_sdpert = 1
              IF(ioutside == 0)THEN
               !CALL LOGLHOOD(objnew1)
                 IF(ISMPPRIOR == 0) CALL LOGLHOOD(objnew1)
                 IF(ISMPPRIOR == 1) CALL LOGLHOOD2(objnew1)
                 logPLratio = (objnew1%logL - obj%logL)*beta_mh
                 CALL RANDOM_NUMBER(ran_uni)
                 IF(ran_uni >= EXP(logPLratio))THEN
                    objnew1 = obj
                    obj%ireject = obj%ireject + 1
                 ELSE
                    obj = objnew1
                    obj%iaccept = obj%iaccept + 1
                 ENDIF
              ELSE
                 objnew1 = obj
                 obj%ireject = obj%ireject + 1
                 ioutside = 0
              ENDIF
                 i_sdpert = 0
           ENDDO
        ENDIF
      ENDIF
   ENDIF ! ICOV if
!!
!! Do Metropolis-Hastings on autoregressive model
!!
   IF(IAR == 1)THEN
     IF(IPV == 1)THEN
      !! Perturb AR model with .25 probability
         CALL RANDOM_NUMBER(ran_uni_ar)
         IF(ran_uni_ar>=0.25_SP)THEN
            DO ipar = 1,NARFP*NWF
              CALL PROPOSAL_ARPV(obj,objnew1,ipar)
!              i_sdpert = 1
               IF(ioutside == 0)THEN
                  IF(ISMPPRIOR == 0) CALL LOGLHOOD(objnew1)
                  IF(ISMPPRIOR == 1) CALL LOGLHOOD2(objnew1)
                  logPLratio = (objnew1%logL - obj%logL)*beta_mh
                  CALL RANDOM_NUMBER(ran_uni)
                  IF(ran_uni >= EXP(logPLratio))THEN
                     objnew1 = obj
                     obj%ireject = obj%ireject + 1
                  ELSE
                     obj = objnew1
                     obj%iaccept = obj%iaccept + 1
                  ENDIF
               ELSE
                  objnew1 = obj
                  obj%ireject = obj%ireject + 1
                  ioutside = 0
               ENDIF     ! ioutside endif
            ENDDO        ! ipar enddo
         ENDIF           ! ran_uni_ar endif
       ENDIF
       IF(IGV == 1)THEN
      !! Perturb AR model with .25 probability
         CALL RANDOM_NUMBER(ran_uni_ar)
         IF(ran_uni_ar>=0.25_SP)THEN
            DO ipar = 1,NARFP*NWF
              CALL PROPOSAL_ARGV(obj,objnew1,ipar)
!              i_sdpert = 1
               IF(ioutside == 0)THEN
                  IF(ISMPPRIOR == 0) CALL LOGLHOOD(objnew1)
                  IF(ISMPPRIOR == 1) CALL LOGLHOOD2(objnew1)
                  logPLratio = (objnew1%logL - obj%logL)*beta_mh
                  CALL RANDOM_NUMBER(ran_uni)
                  IF(ran_uni >= EXP(logPLratio))THEN
                     objnew1 = obj
                     obj%ireject = obj%ireject + 1
                  ELSE
                     obj = objnew1
                     obj%iaccept = obj%iaccept + 1
                  ENDIF
               ELSE
                  objnew1 = obj
                  obj%ireject = obj%ireject + 1
                  ioutside = 0
               ENDIF     ! ioutside endif
            ENDDO        ! ipar enddo
         ENDIF           ! ran_uni_ar endif
      ENDIF
   ENDIF                 ! IAR if
ENDIF                    ! k-change endif
END SUBROUTINE EXPLORE_MH
!=======================================================================

SUBROUTINE CHECKBOUNDS(obj)
USE DATA_TYPE
USE RJMCMC_COM
USE MPI
IMPLICIT NONE
INTEGER(KIND=IB)                            :: ip,iwhich,ipar,ilay,nvpvsr
TYPE(objstruc)                              :: obj
INTEGER(KIND=IB),DIMENSION((kmax+1)*(NPL-1)):: idxvpvsr
DO ilay = 1, obj%k*NPL, NPL
  IF(obj%par(ilay) < minlim(1))ioutside = 1
  IF(obj%par(ilay)  >  maxlim(1))ioutside = 1
ENDDO
DO ilay = 1,obj%k
   IF(hmin > obj%h(ilay))ioutside = 1
   IF(maxlim(1) < obj%h(ilay))ioutside = 1
ENDDO
CALL GETIDXVPVSR(obj,idxvpvsr,nvpvsr)

DO ip = 1,nvpvsr
   iwhich = idxvpvsr(ip)
   ilay = CEILING(REAL(iwhich,SP)/REAL(NPL,SP))
   ipar = iwhich-(ilay-1)*NPL
   IF(ilay > obj%k) ipar = ipar + 1
   IF(((obj%par(iwhich) - minlim(ipar)) < 0._SP).OR.((maxlim(ipar) - obj%par(iwhich)) < 0._SP))THEN
      ioutside = 1
      IF(rank==src)PRINT*,iwhich,obj%par(iwhich), minlim(ipar),maxlim(ipar)
   ENDIF
ENDDO 
RETURN
END SUBROUTINE CHECKBOUNDS
!=======================================================================

SUBROUTINE PROPOSAL(obj,objnew,iwhich,factor)
!=======================================================================
USE DATA_TYPE
USE RJMCMC_COM
IMPLICIT NONE
INTEGER(KIND=IB) :: iwhich,ipar,ilay
TYPE(objstruc) :: obj,objnew
REAL(KIND=SP)  :: ran_uni,factor

objnew = obj
ilay  = CEILING(REAL(iwhich,SP)/REAL(NPL,SP))
ipar  = iwhich-(ilay-1)*NPL

CALL RANDOM_NUMBER(ran_uni)
!! Gaussian:
!objnew%par(iwhich) = obj%par(iwhich) + fact/factor*pertsd(ipar)*ran_nor
!! Cauchy:
objnew%par(iwhich) = obj%par(iwhich) + fact/factor*pertsd(ipar)*TAN(PI2*(ran_uni-0.5_SP))
CALL CALCH(objnew)
201 FORMAT(4F14.6)
RETURN
END SUBROUTINE PROPOSAL
!=======================================================================

SUBROUTINE PROPOSAL_ARPV(obj,objnew,iwhich)
!=======================================================================
USE DATA_TYPE
USE RJMCMC_COM
IMPLICIT NONE
INTEGER(KIND=IB) :: iwhich
TYPE(objstruc) :: obj,objnew
REAL(KIND=SP)  :: ran_nor
!!
!! Gaussian proposal
!!
CALL GASDEVJ(ran_nor)
objnew%arpar_pv(iwhich) = obj%arpar_pv(iwhich) + pertarsd(iwhich)*ran_nor
IF(((objnew%arpar_pv(iwhich) - minlimar(iwhich)) < 0._SP).OR.((maxlimar(iwhich) - objnew%arpar_pv(iwhich)) < 0._SP))ioutside = 1
RETURN
END SUBROUTINE PROPOSAL_ARPV


SUBROUTINE PROPOSAL_ARGV(obj,objnew,iwhich)
!=======================================================================
USE DATA_TYPE
USE RJMCMC_COM
IMPLICIT NONE
INTEGER(KIND=IB) :: iwhich
TYPE(objstruc) :: obj,objnew
REAL(KIND=SP)  :: ran_nor
!!
!! Gaussian proposal
!!
CALL GASDEVJ(ran_nor)
objnew%arpar_gv(iwhich) = obj%arpar_gv(iwhich) + pertarsd(iwhich)*ran_nor
IF(((objnew%arpar_gv(iwhich) - minlimar(iwhich)) < 0._SP).OR.((maxlimar(iwhich)-objnew%arpar_gv(iwhich)) < 0._SP))ioutside = 1
RETURN
END SUBROUTINE PROPOSAL_ARGV

!=======================================================================

SUBROUTINE CHECKBOUNDS_ARMX(obj)
!=======================================================================
USE DATA_TYPE
USE RJMCMC_COM
IMPLICIT NONE
INTEGER(KIND=IB) :: ifr
TYPE(objstruc):: obj

DO ifr = 1,NWF
  IF(IPV == 1)THEN
    IF(MAXVAL(obj%Dar_pv(ifr,:)) > armx)THEN
      iarfail_pv = 1
    ENDIF
    IF(MINVAL(obj%Dar_pv(ifr,:)) < -armx)THEN
      iarfail_pv = 1
    ENDIF
  ENDIF
  IF(IGV == 1)THEN
    IF(MAXVAL(obj%Dar_gv(ifr,:)) > armx)THEN
      iarfail_gv = 1
    ENDIF
    IF(MINVAL(obj%Dar_gv(ifr,:)) < -armx)THEN
      iarfail_gv = 1
    ENDIF
  ENDIF
ENDDO

RETURN
END SUBROUTINE CHECKBOUNDS_ARMX
!=======================================================================


SUBROUTINE PROPOSAL_MISC(obj,objnew,iwhich)
!=======================================================================
USE DATA_TYPE
USE RJMCMC_COM
IMPLICIT NONE
INTEGER(KIND=IB) :: iwhich
TYPE(objstruc) :: obj,objnew
REAL(KIND=SP)  :: ran_nor

!!
!! Gaussian proposal
!!
CALL GASDEVJ(ran_nor)
!PRINT*,'obj%imsc',obj%misc(iwhich),'perturbation',pertmiscsd(iwhich)
objnew%misc(iwhich) = obj%misc(iwhich) + pertmiscsd(iwhich)*ran_nor
IF(((objnew%misc(iwhich) - minlimmisc(iwhich)) < 0._SP).OR.((maxlimmisc(iwhich) - objnew%misc(iwhich)) < 0._SP))ioutside = 1
!PRINT*,'objnew%misc',obj%misc(iwhich),'iwhich',iwhich
RETURN
END SUBROUTINE PROPOSAL_MISC
!=======================================================================

SUBROUTINE PROPOSAL_SDPV(obj,objnew,iwhich)
!=======================================================================
USE DATA_TYPE
USE RJMCMC_COM
IMPLICIT NONE
INTEGER(KIND=IB) :: iwhich
TYPE(objstruc) :: obj,objnew
REAL(KIND=SP)  :: ran_nor
!!
!! Gaussian proposal
!!
CALL GASDEVJ(ran_nor)
objnew%sdpar_pv(iwhich) = obj%sdpar_pv(iwhich) + pertsdsd(iwhich)*ran_nor
IF(((objnew%sdpar_pv(iwhich) - minlimsd(iwhich)) < 0._SP).OR.((maxlimsd(iwhich) - objnew%sdpar_pv(iwhich)) < 0._SP))ioutside = 1
RETURN
END SUBROUTINE PROPOSAL_SDPV

SUBROUTINE PROPOSAL_SDGV(obj,objnew,iwhich)
!=======================================================================^M
USE DATA_TYPE
USE RJMCMC_COM
IMPLICIT NONE
INTEGER(KIND=IB) :: iwhich
TYPE(objstruc) :: obj,objnew
REAL(KIND=SP)  :: ran_nor
!!
!! Gaussian proposal
!!
CALL GASDEVJ(ran_nor)
objnew%sdpar_gv(iwhich) = obj%sdpar_gv(iwhich) + pertsdsd(iwhich)*ran_nor
IF(((objnew%sdpar_gv(iwhich) - minlimsd(iwhich)) < 0._SP).OR.((maxlimsd(iwhich)-objnew%sdpar_gv(iwhich)) < 0._SP))ioutside = 1
RETURN
END SUBROUTINE PROPOSAL_SDGV
!=======================================================================

SUBROUTINE GETIDXVPVSR(obj,idxvpvsr,nvpvsr)
!=======================================================================
USE DATA_TYPE
USE RJMCMC_COM
IMPLICIT NONE
INTEGER(KIND=IB) :: ipar,ilay,ipl,nvpvsr,i
TYPE(objstruc) :: obj
INTEGER(KIND=IB),DIMENSION((kmax+1)*(NPL-1)) :: idxvpvsr

IF(obj%k == 0)THEN
!  idxvpvsr = (/1,2,3/)
  idxvpvsr = (/(i, i=1,NPL-1,1)/)
  nvpvsr = NPL-1
ELSE
  !! Layers:
  idxvpvsr = 0
  ipar = 1
  DO ilay = 1,obj%k
    DO ipl = 2,NPL
      idxvpvsr(ipar) = (ilay-1)*NPL+ipl
      ipar = ipar + 1
    ENDDO
  ENDDO
  !! Half-space: (it was commented before)
  DO ipl = 1,NPL-1
    idxvpvsr(ipar) = obj%k*NPL+ipl
    ipar = ipar + 1
  ENDDO
  nvpvsr = (obj%k+1)*(NPL-1)
! nvpvsr = (obj%k)*(NPL-1)
ENDIF
END SUBROUTINE GETIDXVPVSR
!=======================================================================

SUBROUTINE CALCH(obj)
!THIS SUBROUTINE COMPUTES THE THICKNESS OF EACH LAYER AFTER DIVIDING THE
!A LAYER BY A NEW INTERFACE
!=======================================================================
USE DATA_TYPE
USE RJMCMC_COM
IMPLICIT NONE
INTEGER(KIND=IB) :: i
TYPE(objstruc) :: obj
INTEGER(KIND=IB),DIMENSION(obj%k) :: idxz

IF(obj%k == 0)THEN
  obj%h = 0._SP
  obj%z = 0._SP
ELSE
  idxz = 0
  idxz = (/(i, i=1,obj%k*NPL,NPL)/)
  obj%z = 0._SP
  obj%z(1:obj%k) = obj%par(idxz)
  obj%h = 0._SP
  obj%h(1) = obj%par(1)
  DO i = 2,obj%k
    obj%h(i) = obj%par(idxz(i))-obj%par(idxz(i-1))
  ENDDO
ENDIF
END SUBROUTINE CALCH
!=======================================================================

SUBROUTINE DEATH(obj,objnew)
!=======================================================================
USE DATA_TYPE
USE RJMCMC_COM
IMPLICIT NONE
INTEGER(KIND=IB)                 :: idel,idel2
TYPE(objstruc)                   :: obj,objnew
INTEGER(KIND=IB),DIMENSION(NFPMX):: idxdeath
REAL(KIND=SP)                    :: ran_uni
REAL(KIND=SP),DIMENSION(NPL-1)   :: vpvsr

objnew%k   = obj%k - 1
objnew%NFP = (objnew%k * NPL) + (NPL-1)
IF(objnew%k >  0)THEN
!! Pick random layer:
  idxdeath = 0
  idxdeath(1:obj%k) = RANDPERM(obj%k)
  idel = idxdeath(1)

!! Delete interface (setting new layer properties to one of old layers
!! properties at random ):
  CALL RANDOM_NUMBER(ran_uni)
  idel2 = idel
  IF(ran_uni>=0.5_SP)idel2 = idel+1

  objnew%par = 0._SP  ! Take care here that there are no old "left overs" from other dimensions
  vpvsr = 0._SP
  IF(idel == 1)THEN
!! Killed interface is first interface:
     vpvsr = obj%par((idel2-1)*NPL+2:(idel2-1)*NPL+NPL)
     objnew%par(1:objnew%NFP) = obj%par(NPL+1:obj%NFP)
     objnew%par(2:NPL) = vpvsr
  ELSEIF(idel == obj%k)THEN
!! Killed interface is last interface:
     IF(idel2 == obj%k)vpvsr = obj%par((idel2-1)*NPL+2:(idel2-1)*NPL+NPL)
     IF(idel2 > obj%k) vpvsr = obj%par(obj%NFP-(NPL-2):obj%NFP)
     objnew%par(1:objnew%NFP-(NPL-1)) = (/ obj%par(1:obj%NFP-(NPL*2-1)) /)
     objnew%par(objnew%NFP-(NPL-2):objnew%NFP) = vpvsr
  ELSE
!! Killed interface is in normal layer stack:
     vpvsr = obj%par((idel2-1)*NPL+2:(idel2-1)*NPL+NPL)
     objnew%par(1:objnew%NFP) = (/ obj%par(1:(idel-1)*NPL),obj%par(idel*NPL+1:obj%NFP) /)
     objnew%par((idel-1)*NPL+2:(idel-1)*NPL+NPL) = vpvsr
  ENDIF
ELSE !else objnew%k
  objnew%par = 0._SP
  idel = 1
  CALL RANDOM_NUMBER(ran_uni)
  IF(ran_uni >= 0.5) idel = idel+1
  objnew%par(1:objnew%NFP) = obj%par(idel*(NPL-1)+1-(NPL-2):idel*(NPL-1)+1)
ENDIF
CALL CALCH(objnew)
!PRINT*,'old',idel,idel2
!CALL PRINTPAR(obj)
!PRINT*,'new'
!CALL PRINTPAR(objnew)
!IF(idel2 == 1)STOP

END SUBROUTINE DEATH
!=======================================================================

SUBROUTINE BIRTH(obj,objnew)
!=======================================================================
USE DATA_TYPE
USE RJMCMC_COM
IMPLICIT NONE
INTEGER(KIND=IB)              :: i, iznew, ipert, iloop, iwhich, ipar
TYPE(objstruc)                :: obj,objnew,objnew2
REAL(KIND=SP)                 :: znew,hnew1,ran_uni
REAL(KIND=SP),DIMENSION(obj%k):: ztmp

objnew%k   = obj%k + 1
objnew%NFP = (objnew%k * NPL) + (NPL-1)
!!
!! Draw new z until new layer > hmin
!!
CALL RANDOM_NUMBER(ran_uni)
znew = maxpert(1)*ran_uni
IF(objnew%k > 1)THEN
  ztmp = obj%z(1:obj%k) - znew
  IF(hmin > MINVAL(ABS(ztmp))) ioutside = 1
  IF(hmin > znew) ioutside = 1
!!
!! Insert new interface
!!
  hnew1 = 0
  iznew = 0
  DO i = 1,obj%k
     IF(ztmp(i) > 0._SP) EXIT
     iznew = i
  ENDDO
  iznew = iznew + 1
  objnew%z   = 0._SP
  objnew%par = 0._SP
  IF(iznew == 1)THEN
!! New layer is first layer:
     objnew%par(1:objnew%NFP) = (/ obj%par(1:NPL), obj%par /)
     objnew%par(1)            = znew
  ELSEIF(iznew > obj%k)THEN
   !! New layer is created in half-space:
     objnew%par(1:objnew%NFP) = (/ REAL(obj%par(1:obj%NFP-(NPL-1)),SP),REAL(hnew1,SP), &
                                 REAL(obj%par(obj%NFP-(NPL-2):obj%NFP),SP), &
                                 REAL(obj%par(obj%NFP-(NPL-2):obj%NFP),SP) /)
     objnew%par(objnew%NFP-((NPL*2)-2)) = znew
  ELSE
   !! New layer is within layer stack:
     objnew%par(1:objnew%NFP) = (/ obj%par(1:(iznew)*NPL), &
                                 obj%par(((iznew-1)*NPL)+1:iznew*NPL), &
                                 obj%par(iznew*NPL+1:obj%NFP) /)
     objnew%par(((iznew-1)*NPL)+1) = znew
  ENDIF
ELSE ! else objnew%k
  iznew = 1
  objnew%z   = 0._SP
  objnew%par = 0._SP
  objnew%par(1) = znew
  objnew%par(2:objnew%NFP) = (/ obj%par(1:NPL-1), obj%par /)
ENDIF ! end objnew%k loop
 
CALL CALCH(objnew)
!!
!! Pick one of the two new layers at random and perturb:
!!
CALL RANDOM_NUMBER(ran_uni)
IF(ran_uni <= 0.5_SP)THEN
   ipert = iznew
ELSE
   ipert = iznew+1
ENDIF

!! Perturb only one parameter (determined in EXPLORE_MH call)
objnew2 = objnew
DO ipar = 2,NPL
   iloop = 0
   iwhich = (ipert-1)*NPL+ipar
   IF(ipert > objnew%k) iwhich = iwhich - 1
   !!
   !! Draw from unifomr prior (so that proposal and prior ratios cancel in BD acceptance):
   !!
   CALL RANDOM_NUMBER(ran_uni)
   objnew2%par(iwhich) = minlim(ipar) + maxpert(ipar)*ran_uni

   IF(((objnew2%par(iwhich) - minlim(ipar)) < 0._SP).OR. &
      ((maxlim(ipar) - objnew%par(iwhich)) < 0._SP)) ioutside = 1
ENDDO
objnew = objnew2
!PRINT*,'debug -- before birth'
!CALL PRINTPAR(obj)
!PRINT*,'debug---sp5 after birth'
!CALL PRINTPAR(objnew)

END SUBROUTINE BIRTH

!=======================================================================
!!
!=======================================================================
SUBROUTINE ARFB_PV(obj,idata,idatb)
!=======================================================================
!! Autoregressive model to model data error correlations.
!! AR process is computed forward and backward and average is used.
!!
USE MPI
USE RJMCMC_COM
IMPLICIT NONE
TYPE (objstruc)  :: obj
INTEGER          :: i,j,k,ifr,idata,idatb
REAL(KIND=SP),DIMENSION(idatb-idata+1)::dres1,dar1
DO ifr=1,NWF
   k = NARFP
   obj%Dar_pv(ifr,idata)=0._SP          ! Matlab sets first point to zero...
   !!
   !! Real part:
   !!
   dres1 = 0._SP
   dres1 = obj%res_pv(ifr,idata:idatb)
   dar1(1)=0._SP          ! Matlab sets first point to zero...
   DO i=2,idatb-idata+1
      dar1(i) = 0.0
      IF(k >= i)THEN
         DO j=1,i-1
            dar1(i) = dar1(i) + obj%arpar_pv(NARFP*(ifr-1)+j) * dres1(i-j)
         ENDDO
      ELSE
         DO j=1,k
            dar1(i) = dar1(i) + obj%arpar_pv(NARFP*(ifr-1)+j) * dres1(i-j)
         ENDDO
      ENDIF
   ENDDO
   obj%Dar_pv(ifr,idata:idatb) = dar1
   obj%Dar_pv(ifr,idata) = 0._SP
   obj%Dar_pv(ifr,idatb) = 0._SP
ENDDO
END SUBROUTINE ARFB_PV
!=======================================================================

SUBROUTINE ARFB_GV(obj,idata,idatb)
!=======================================================================
!! Autoregressive model to model data error correlations.
!! AR process is computed forward and backward and average is used.
!!
USE MPI
USE RJMCMC_COM
IMPLICIT NONE
TYPE (objstruc)  :: obj
INTEGER          :: i,j,k,ifr,idata,idatb
REAL(KIND=SP),DIMENSION(idatb-idata+1)::dres1,dar1
DO ifr=1,NWF
   k = NARFP
   obj%Dar_gv(ifr,idata)=0._SP          ! Matlab sets first point to zero...
   !!
   !! Real part:
   !!
   dres1 = 0._SP
   dres1 = obj%res_gv(ifr,idata:idatb)
   dar1(1)=0._SP          ! Matlab sets first point to zero...
   DO i=2,idatb-idata+1
      dar1(i) = 0.0
      IF(k >= i)THEN
         DO j=1,i-1
            dar1(i) = dar1(i) + obj%arpar_gv(NARFP*(ifr-1)+j) * dres1(i-j)
         ENDDO
      ELSE
         DO j=1,k
            dar1(i) = dar1(i) + obj%arpar_gv(NARFP*(ifr-1)+j) * dres1(i-j)
         ENDDO
      ENDIF
   ENDDO
   obj%Dar_gv(ifr,idata:idatb) = dar1
   obj%Dar_gv(ifr,idata) = 0._SP
   obj%Dar_gv(ifr,idatb) = 0._SP
ENDDO
END SUBROUTINE ARFB_GV
!=======================================================================


SUBROUTINE PARALLEL_SEED()
!!
!!  Ensure unique random seed for each CPU
!!
!=======================================================================
USE MPI
USE RJMCMC_COM
IMPLICIT NONE

INTEGER(KIND=IB) :: i
INTEGER(KIND=IB), DIMENSION(:),   ALLOCATABLE :: iseed
INTEGER(KIND=IB)                              :: iseedsize
INTEGER(KIND=IB), DIMENSION(:,:), ALLOCATABLE :: iseeds
REAL(KIND=SP),    DIMENSION(:,:), ALLOCATABLE :: rseeds
INTEGER(KIND=IB), DIMENSION(:),   ALLOCATABLE :: iseed1
REAL(KIND=SP) :: ran_uni

CALL RANDOM_SEED   !INITIALIZED A PSEUDO-RANDOM NUMBER SEQUENCE TO A DEFAULT STATE
CALL RANDOM_SEED(SIZE=iseedsize)  !MINIMUM SIZE OF THE ARRAY USED 
IF(iseedsize <= 0)THEN
  STOP
ENDIF
ALLOCATE( iseed1(iseedsize) )
IF(ISETSEED == 1)THEN
   iseed1 = (/2303055,     2435432,     5604058,     4289794,     3472290, &
      7717070,      141180,     3783525,     3087889,     4812786,     3028075, &
      3712062,     6316731,      436800,     7957708,     2055697,     1944360, &
      1222992,     7537775,     7769874,     5588112,     7590383,     1426393, &
      1753301,     7681841,     2842400,     4411488,     7304010,      497639, &
      4978920,     5345495,      754842,     7360599,     5776102/)
   CALL RANDOM_SEED(PUT=iseed1)
ELSE
   CALL RANDOM_SEED(GET=iseed1)
!   IF(rank==src)WRITE(6,*) 'Master seed:',iseed1
ENDIF
ALLOCATE( iseed(iseedsize), rseeds(iseedsize,NTHREAD), iseeds(iseedsize,NTHREAD) )
iseed = 0
rseeds = 1._SP
iseeds = 0
IF(rank == src)THEN
   CALL RANDOM_NUMBER(rseeds)
   iseeds = -NINT(rseeds*1000000._SP)
ENDIF
DO i = 1,iseedsize
   CALL FLUSH(6)
   CALL MPI_BCAST( iseeds(i,:), NTHREAD, MPI_INTEGER, src, MPI_COMM_WORLD, ierr )
ENDDO
iseed = iseeds(:,rank+1)
!!
!! Write seeds to seed logfile:
!!
IF(rank == src)THEN
   OPEN(UNIT=50,FILE=seedfile,FORM='formatted',STATUS='UNKNOWN', &
   ACTION='WRITE',POSITION='REWIND',RECL=1024)
   WRITE(50,*) 'Rank: ',rank
   WRITE(50,201) iseed
   WRITE(50,*) ''
   CLOSE(50)
ENDIF
CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
DO i = 1,NTHREAD-1
   IF(rank == i)THEN
      OPEN(UNIT=50,FILE=seedfile,FORM='formatted',STATUS='UNKNOWN', &
      ACTION='WRITE',POSITION='APPEND',RECL=1024)
      WRITE(50,*) 'Rank: ',rank
      WRITE(50,201) iseed
      WRITE(50,*) ''
      CLOSE(50)
   ENDIF
   CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
ENDDO
CALL RANDOM_SEED(PUT=iseed)

DO i = 1,100
   CALL RANDOM_NUMBER(ran_uni)
ENDDO
DO i = 1,CEILING(ran_uni*10000._SP)
   CALL RANDOM_NUMBER(ran_uni)
ENDDO
201   FORMAT(50I10)
END SUBROUTINE PARALLEL_SEED

!=======================================================================


SUBROUTINE SAVESAMPLE(objm,ikeep,isource1,isource2,tmaster)
!!=======================================================================
!!
!! Exchanging and saving posterior samples
!!
USE MPI
USE RJMCMC_COM
IMPLICIT NONE
TYPE (objstruc),DIMENSION(2)     :: objm     ! Objects on master node
INTEGER(KIND=IB)                 :: ic,ikeep,jidx,ik
INTEGER(KIND=IB)                 :: isource,isource1,isource2
!INTEGER(KIND=SP),DIMENSION(NFPMX):: tmpprop,tmpacc
REAL(KIND=SP) :: tmaster

!!
!! Save object into global sample array
!!
DO ic = 1,2
  IF(ic == 1)isource = isource1
  IF(ic == 2)isource = isource2
  IF(objm(ic)%beta > 1._SP/dTlog)THEN
    !!
    !! Write to stdout
    !!
    IF(MOD(imcmc1,20*NKEEP) == 0)THEN
      tsave2 = MPI_WTIME()
                        WRITE(ulog,213)'          iaccept_bd = ',objm(ic)%iaccept_bd,objm(ic)%beta
      IF(IEXCHANGE == 1)WRITE(ulog,215)'       T swap accept = ',REAL(ncswap,SP)/REAL(ncswapprop,SP),isource1,isource2
                        WRITE(ulog,216)'Total time for block = ',tsave2-tsave1,' s'
      WRITE(ulog,*) ''
      WRITE(ulog,203) '     imcmc1,         logL,  k,  iacc_bd, acceptance  source, time(slave)'
      WRITE(ulog,203) '------------------------------------------------------------------------'
      tsave1 = MPI_WTIME()
    ENDIF
!      tmpprop(ipar:ipar+NPL-1) = objm(ic)%iproposevoro(ivo,:)
!      tmpacc(ipar:ipar+NPL-1)  = objm(ic)%iacceptvoro(ivo,:)
    
    sample(ikeep,:) =  (/ REAL(objm(ic)%logL,SP), REAL(1._SP,SP), REAL(objm(ic)%beta,SP), REAL(objm(ic)%k,SP), & ! 4 parameters
                          objm(ic)%par,objm(ic)%misc,objm(ic)%sdpar_pv, objm(ic)%sdpar_gv,  & 
                          objm(ic)%arpar_pv,objm(ic)%arpar_gv,REAL(objm(ic)%acceptance,SP), &
                       REAL(objm(ic)%iaccept_bd,SP),REAL(objm(ic)%ireject_bd,SP),REAL(i_bd,SP),REAL(ic,SP),REAL(rank,SP)/)

    IF(MOD(imcmc1,NKEEP)==0)THEN
      WRITE(ulog,202) '       ',imcmc1,objm(ic)%logL,objm(ic)%k,objm(ic)%iaccept_bd,objm(ic)%acceptance,isource,objm(ic)%tcmp
    ENDIF
    ikeep = ikeep + 1_IB
    imcmc1 = imcmc1 + 1_IB

    !!
    !! Write to file
    !!
    IF(ikeep > NKEEP)THEN
      DO jidx=1,NKEEP
        WRITE(usample,207) sample(jidx,:)
      ENDDO
      CALL FLUSH(usample)
      !WRITE(ustep,207) pertsd(1:NPL)
      !CALL FLUSH(ustep)
      sample  = 0._SP
      ikeep   = 1
    ENDIF
  ENDIF
ENDDO
CALL FLUSH(6)
202 FORMAT(A3,I8,1F14.4,I4,I10,f12.5,I8,f12.4)
203 FORMAT(A69)
207 FORMAT(500ES18.8)
213 FORMAT(a23,I3,F8.4)
214 FORMAT(a23,I6,a9,I4)
215 FORMAT(a23,F8.4,I4,I4)
216 FORMAT(a23,F8.2,a)
RETURN
END SUBROUTINE SAVESAMPLE
!=======================================================================

SUBROUTINE SAVEREPLICA(obj,filename2,filename3,filename5,filename6)
!=======================================================================
USE MPI
USE RJMCMC_COM
IMPLICIT NONE

INTEGER(KIND=IB) :: i,j
TYPE (objstruc)  :: obj      ! Best object
CHARACTER(len=64):: filename2,filename3,filename5,filename6 ! replica file name
WRITE(ulog,*) 'MAP model:'
CALL PRINTPAR(obj)
WRITE(ulog,*) 'logL = ',obj%logL
OPEN(UNIT=50,FILE=filename2,FORM='formatted',STATUS='REPLACE',ACTION='WRITE',RECL=1024)
OPEN(UNIT=51,FILE=filename5,FORM='formatted',STATUS='REPLACE',ACTION='WRITE',RECL=1024)
DO i = 1,NWF
  DO j = 1,obj%NDAT_pv
    WRITE(50,208) obj%t_pv(i,j),obj%Drep_pv(i,j)
    WRITE(51,208) obj%t_pv(i,j),obj%Dobs_pv(i,j)
  ENDDO
ENDDO
CLOSE(50)
CLOSE(51)
OPEN(UNIT=50,FILE=filename3,FORM='formatted',STATUS='REPLACE',ACTION='WRITE',RECL=1024)
OPEN(UNIT=51,FILE=filename6,FORM='formatted',STATUS='REPLACE',ACTION='WRITE',RECL=1024)
DO i = 1,NWF
 DO j = 1,obj%NDAT_gv
   WRITE(50,208) obj%t_gv(i,j),obj%Drep_gv(i,j)
   WRITE(51,208) obj%t_gv(i,j),obj%Dobs_gv(i,j)
 ENDDO
ENDDO
CLOSE(50)
CLOSE(51)

IF(IAR == 1)THEN
  OPEN(29,file='ar_pv_vector.dat',status='replace')
  OPEN(30,file='ar_gv_vector.dat',status='replace')    
  DO I = 1, NWF
    DO j = 1,obj%NDAT_pv
       WRITE(29,*)obj%Dar_pv(i,j)
    ENDDO
    DO j = 1,obj%NDAT_gv
      WRITE(30,*)obj%Dar_gv(i,j)
    ENDDO
  ENDDO
ENDIF
CLOSE(29)
CLOSE(30)
208 FORMAT(2ES20.10)
RETURN
END SUBROUTINE SAVEREPLICA
!=======================================================================

SUBROUTINE MAKE_BETA(beta1,beta4,NTEMP1)
!!
!! Make cooling schedule (three geometrically spaced legs)
!!
!=======================================================================
USE MPI
USE RJMCMC_COM
IMPLICIT NONE
INTEGER(KIND=IB) :: itemp,NTEMP1
REAL(KIND=SP)                  :: beta1,beta2,beta3,beta4
REAL(KIND=SP),DIMENSION(NTEMP1):: logbeta
beta2 = beta1+0.2_SP*ABS(beta1-beta4)
beta3 = beta1+0.75_SP*ABS(beta1-beta4)
logbeta = 0._SP
beta = 0._SP
logbeta(1) = LOG(beta1)
DO itemp=2,NTEMP1/20
   logbeta(itemp) = logbeta(itemp-1)+(LOG(beta2)-LOG(beta1))/(REAL(NTEMP1/20,SP)-1)
ENDDO
DO itemp=NTEMP1/20+1,NTEMP1/3
   logbeta(itemp)=logbeta(itemp-1)+(LOG(beta3)-LOG(beta2))/(REAL(NTEMP1/3,SP)-REAL(NTEMP1/20,SP))
ENDDO
DO itemp=NTEMP1/3+1,NTEMP1
   logbeta(itemp)=logbeta(itemp-1)+(LOG(beta4)-LOG(beta3))/(REAL(NTEMP1,SP)-REAL(NTEMP1/3,SP))
ENDDO
beta = EXP(logbeta)

209 FORMAT(50ES16.8)
RETURN
END SUBROUTINE MAKE_BETA
!=======================================================================

SUBROUTINE PRINTPAR(obj)
!=======================================================================
USE DATA_TYPE
USE RJMCMC_COM
IMPLICIT NONE
INTEGER :: i,npar
TYPE(objstruc) :: obj
REAL(KIND=SP)  :: tmppar(obj%NFP)
WRITE(ulog,*)'Parameters on the layers bolow the Suraface to the k-th layer'
!IF(IMAP == 1)THEN
!  DO i=1,obj%k
!     WRITE(ulog,201) obj%par((i-1)*NPL+1:i*NPL)
!  ENDDO
!ELSE
IF(obj%k>0)THEN
  WRITE(ulog,*)'Parameters on the layers, depth, vs, perturbation'
  DO i=1,obj%k
   WRITE(ulog,201) obj%par((i-1)*NPL+1:i*NPL)
  ENDDO
ENDIF
WRITE(ulog,*)'Parameters on the half space, vs perturbation'
WRITE(ulog,202) '            ',obj%par(obj%k*NPL+1:obj%k*NPL+(NPL-1))

tmppar = 0.0
CALL ADD_REF_MDL(obj,tmppar)
IF(obj%k>0)THEN
  WRITE(ulog,*)'Parameters on the layers above half space after adding reference model'
  DO i=1,obj%k
   WRITE(ulog,201) tmppar((i-1)*NPL+1:i*NPL)
  ENDDO
ENDIF
WRITE(ulog,*)'Parameters on the half space after adding reference model'
WRITE(ulog,202) '            ',tmppar(obj%k*NPL+1:obj%k*NPL+(NPL-1))

WRITE(ulog,*) 'Partition:'
WRITE(ulog,203) obj%z
WRITE(ulog,*) 'Layers:'
WRITE(ulog,203) obj%h

IF(ICOV == 1)THEN
   WRITE(ulog,*) 'SD parameters:'
   IF(IPV==1) WRITE(ulog,203) obj%sdpar_pv
   IF(IGV==1) WRITE(ulog,203) obj%sdpar_gv
ENDIF
IF(IAR == 1)THEN
   WRITE(ulog,*) 'AR parameters:'
   IF(IPV==1) WRITE(ulog,203) obj%arpar_pv
   IF(IGV==1)  WRITE(ulog,203) obj%arpar_gv
ENDIF
IF(IMISC == 1)THEN
  WRITE(ulog,*) 'MISC parameters: Gaussian filter width'
  WRITE(ulog,*) obj%misc
ENDIF

201 FORMAT(5F12.4)
202 FORMAT(A12,8F12.4)
203 FORMAT(11F12.4)
205 FORMAT(11I3)
204 FORMAT(4F16.4)
END SUBROUTINE PRINTPAR
!=======================================================================

FUNCTION RANDPERM(num)
!==============================================================================
USE data_type, ONLY : IB, RP, SP
IMPLICIT NONE
INTEGER(KIND=IB), INTENT(IN) :: num
INTEGER(KIND=IB) :: numb, i, j, k
INTEGER(KIND=IB), DIMENSION(num) :: randperm
REAL(KIND=SP), DIMENSION(num) :: rand2
INTRINSIC RANDOM_NUMBER
CALL RANDOM_NUMBER(rand2)
DO i=1,num
   numb=1
   DO j=1,num
      IF (rand2(i) > rand2(j)) THEN
           numb=numb+1
      END IF
   END DO
   DO k=1,i-1
      IF (rand2(i) <= rand2(k) .AND. rand2(i) >= rand2(k)) THEN
           numb=numb+1
      END IF
   END DO
   randperm(i)=numb
END DO
RETURN
END FUNCTION RANDPERM
!====================================================================

SUBROUTINE GASDEVJ(harvest)
!====================================================================
USE DATA_TYPE
!USE nrtype
USE nr
IMPLICIT NONE
REAL(SP), INTENT(OUT) :: harvest
REAL(SP) :: rsq,v1,v2
REAL(SP), SAVE :: g
LOGICAL, SAVE :: gaus_stored=.FALSE.
IF (gaus_stored) THEN
   harvest=g
   gaus_stored=.FALSE.
ELSE
   DO
      CALL RANDOM_NUMBER(v1)
      CALL RANDOM_NUMBER(v2)
      v1=2.0_SP*v1-1.0_SP
      v2=2.0_SP*v2-1.0_SP
      rsq=v1**2+v2**2
      IF (rsq > 0.0_SP .AND. rsq < 1.0_SP) EXIT
   END DO
   rsq=SQRT(-2.0_SP*LOG(rsq)/rsq)
   harvest=v1*rsq
   g=v2*rsq
   gaus_stored=.TRUE.
END IF
END SUBROUTINE GASDEVJ
!=======================================================================

!Function LOGPLUS(x,y)
!=======================================================================
!USE DATA_TYPE, ONLY : IB, RP
!REAL(KIND=SP)    :: logplus,x,y

!IF(x > y)THEN
!   LOGPLUS = x+LOG(1.0_SP+EXP(y-x))
!ELSE
!   LOGPLUS = y+log(1.0_SP+EXP(x-y))
!ENDIF

!RETURN
!END FUNCTION LOGPLUS
!=======================================================================

SUBROUTINE CHOLESKY(A,n,L)
!=======================================================================
!  Cholesky decomposition of symmetric, positive-definite matix A
!  into lower-triangular matrix L. Modified from Numerical Recipes.
!-----------------------------------------------------------------------

IMPLICIT NONE
INTEGER      :: n,i,j,k
REAL         :: A(n,n),A2(n,n),L(n,n),p(n),summ

A2 = A
DO i=1,n
   DO j=1,n
      IF (i > j) A2(i,j) = 0.
   ENDDO
ENDDO

DO i=1,n
   DO j=i,n
      summ = A2(i,j)
      DO k=i-1,1,-1
         summ = summ-A2(i,k)*A2(j,k)
      ENDDO
      IF (i == j) THEN
        IF (summ <= 0.) THEN
          WRITE(6,*) 'Cholesky Decomp Failed ',summ
          STOP
        ENDIF
        p(i) = SQRT(summ)
      ELSE
        A2(j,i) = summ/p(i)
      ENDIF
  ENDDO
ENDDO

DO i=1,n
   DO j=1,n
      IF (i > j) L(i,j) = A2(i,j)
      IF (i == j) L(i,i) = p(i)
   ENDDO
ENDDO

RETURN
END SUBROUTINE     
!=======================================================================

SUBROUTINE ALLOC_OBJ(obj)
!!==============================================================================
!!
!! Reads observed data.
!!
USE RJMCMC_COM
IMPLICIT NONE
TYPE (objstruc) :: obj

ALLOCATE( obj%par((NLMX*NPL)+NPL-1) )
ALLOCATE( obj%z(NLMX) )
ALLOCATE( obj%h(NLMX) )
ALLOCATE( obj%arpar_pv(NARFP*NWF) )
ALLOCATE( obj%arpar_gv(NARFP*NWF) )
ALLOCATE( obj%misc(NMISC) )
ALLOCATE( obj%sdpar_pv(NWF) )
ALLOCATE( obj%sdpar_gv(NWF) )
ALLOCATE( obj%sdave(NWF) )
ALLOCATE( obj%Dobs_pv(NWF,obj%NDAT_pv) )
ALLOCATE( obj%Dobs_gv(NWF,obj%NDAT_gv) )
ALLOCATE( obj%t_pv(NWF,obj%NDAT_pv) )
ALLOCATE( obj%t_gv(NWF,obj%NDAT_gv) )
ALLOCATE( obj%Drep_pv(NWF,obj%NDAT_pv) )
ALLOCATE( obj%Drep_gv(NWF,obj%NDAT_gv) )
ALLOCATE( obj%res_pv(NWF,obj%NDAT_pv) )
ALLOCATE( obj%res_gv(NWF,obj%NDAT_gv) )
ALLOCATE( obj%Dar_pv(NWF,obj%NDAT_pv) )
ALLOCATE( obj%Dar_gv(NWF,obj%NDAT_gv) )
obj%par = 0._SP
obj%z = 0._SP
obj%h = 0._SP
obj%arpar_pv = 0._SP
obj%arpar_gv = 0._SP
obj%misc = 0._SP
obj%sdpar_pv = 0._SP
obj%sdpar_gv = 0._SP
obj%sdave = 0._SP
obj%Dobs_pv = 0._SP
obj%Dobs_gv = 0._SP
obj%t_pv = 0._SP
obj%t_gv = 0._SP
obj%Drep_pv = 0._SP
obj%Drep_gv = 0._SP
obj%res_pv = 0._SP
obj%res_gv = 0._SP
obj%Dar_pv = 0._SP
obj%Dar_gv = 0._SP

END SUBROUTINE ALLOC_OBJ
!=======================================================================

SUBROUTINE MAKE_MPI_STRUC_SP(obj,objtype)
!!=======================================================================
!!
!! 
!!
USE MPI
USE RJMCMC_COM
IMPLICIT NONE
TYPE (objstruc) :: obj
INTEGER(KIND=IB) :: ifield,objtype
INTEGER(KIND=IB) :: oldtypes(NFIELD), blockcounts(NFIELD)
INTEGER(KIND=MPI_ADDRESS_KIND) :: offsets(NFIELD)

!!
!!
!!                  k,    NFP, ireject_bd iaccept_bd  iaccept ireject  NDAT_pv  NDAT_gv 
   blockcounts = (/ 1,    1,      1,         1,         1,      1,      1,       1,  &  
!                       par             z       h        arpar_pv     arpar_gv        
                   (NLMX*NPL)+NPL-1,   NLMX , NLMX,     NARFP*NWF,  NARFP*NWF,  &
!!                    misc  sdpar_pv sdpar_gv sdave  beta   logL     Dobs_pv 
                     NMISC,   NWF,    NWF,     NWF,   1,     1,  NWF*obj%NDAT_pv,&
!!                      Dobs_gv            t_pv              t_gv           Drep_pv
                    NWF*obj%NDAT_gv,  NWF*obj%NDAT_pv,  NWF*obj%NDAT_gv, NWF*obj%NDAT_pv, &
!                        Drep_gv          res_pv         res_gv            Dar_pv  
                    NWF*obj%NDAT_gv,  NWF*obj%NDAT_pv, NWF*obj%NDAT_gv, NWF*obj%NDAT_pv, &
!!                        Dar_gv        tcmp, acceptance
                     NWF*obj%NDAT_gv,     1,     1 /)

!! Oldtypes

!  INTEGER                                  :: k                  ! Layerdimension (number of layers)
!  INTEGER                                  :: NFP                ! Numberforward parameters
!  INTEGER                                  :: ireject_bd = 0     ! Rejectbirth/death
!  INTEGER                                  :: iaccept_bd = 0     ! Acceptbirth/death
!  INTEGER                                  :: iaccept    = 0     ! MHS accept
!  INTEGER                                  :: ireject    = 0     ! MHS reject
!  INTEGER                                  :: NDAT_pv = 0        ! No. of datapoints for rf
!  INTEGER                                  :: NDAT_gv = 0        ! No. of datapoints for rf
!  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:)   :: par                ! Forwardparameters for layers above half space
!  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:)   :: z                  ! Forwardparameters (interface positions)
!  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:)   :: h                  ! Depthparameters
!  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:)   :: arpar_pv           ! AR for PVforward parameters
!  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:)   :: arpar_gv           ! AR for GVforward parameters
!  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:)   :: misc               ! Miscellaneousparameter (for e.g., gauss filter width)
!  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:)   :: sdpar_pv           ! Std dev forPV forward parameters
!  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:)   :: sdpar_gv           ! Std dev forGV forward parameters
!  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:)   :: sdave              ! Std dev of1000/.. of sdpar
!  REAL(KIND=SP)                            :: beta               ! Temperedlevel
!  REAL(KIND=SP)                            :: logL               ! Loglikelihood
!  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:,:) :: Dobs_pv            ! Observedphase velocity
!  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:,:) :: Dobs_gv            ! Observedgroup velocity
!  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:,:) :: t_pv               !Observed/samples times phase velocity
!  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:,:) :: t_gv               !Observed/samples times group velocity
!  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:,:) :: Drep_pv            ! Replica datafor phase velocity
!  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:,:) :: Drep_gv            ! Replica datafor group velocity
!  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:,:) :: res_pv             ! Dataresiduals for phase velocity
!  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:,:) :: res_gv             ! Dataresiduals for group velocity
!  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:,:) :: Dar_pv             !Autoregressive model replica for phase velocity
!  REAL(KIND=SP),ALLOCATABLE,DIMENSION(:,:) :: Dar_gv             !Autoregressive model replica for group velocity
!  REAL(KIND=SP)                            :: tcmp               ! 
!  REAL(KIND=SP)                            :: acceptance = 0._SP ! Accepancerate for proposal

oldtypes(1:8)    = MPI_INTEGER
oldtypes(9:31)   = MPI_REAL
call mpi_get_address(obj%k,offsets(1),ierr)
call mpi_get_address(obj%NFP,offsets(2),ierr)
call mpi_get_address(obj%ireject_bd,offsets(3),ierr)
call mpi_get_address(obj%iaccept_bd,offsets(4),ierr)
call mpi_get_address(obj%iaccept,offsets(5),ierr)
call mpi_get_address(obj%ireject,offsets(6),ierr)
call mpi_get_address(obj%NDAT_pv,offsets(7),ierr)
call mpi_get_address(obj%NDAT_gv,offsets(8),ierr)
call mpi_get_address(obj%par,offsets(9),ierr)
call mpi_get_address(obj%z,offsets(10),ierr)
call mpi_get_address(obj%h,offsets(11),ierr)
call mpi_get_address(obj%arpar_pv,offsets(12),ierr)
call mpi_get_address(obj%arpar_gv,offsets(13),ierr)
call mpi_get_address(obj%misc,offsets(14),ierr)
call mpi_get_address(obj%sdpar_pv,offsets(15),ierr)
call mpi_get_address(obj%sdpar_gv,offsets(16),ierr)
call mpi_get_address(obj%sdave,offsets(17),ierr)
call mpi_get_address(obj%beta,offsets(18),ierr)
call mpi_get_address(obj%logL,offsets(19),ierr)
call mpi_get_address(obj%Dobs_pv,offsets(20),ierr)
call mpi_get_address(obj%Dobs_gv,offsets(21),ierr)
call mpi_get_address(obj%t_pv,offsets(22),ierr)
call mpi_get_address(obj%t_gv,offsets(23),ierr)
call mpi_get_address(obj%Drep_pv,offsets(24),ierr)
call mpi_get_address(obj%Drep_gv,offsets(25),ierr)
call mpi_get_address(obj%res_pv,offsets(26),ierr)
call mpi_get_address(obj%res_gv,offsets(27),ierr)
call mpi_get_address(obj%Dar_pv,offsets(28),ierr)
call mpi_get_address(obj%Dar_gv,offsets(29),ierr)
call mpi_get_address(obj%tcmp,offsets(30),ierr)
call mpi_get_address(obj%acceptance,offsets(31),ierr)
DO ifield=2,SIZE(offsets)
  offsets(ifield) = offsets(ifield) - offsets(1)
ENDDO
offsets(1) = 0
!IF(rank == src)PRINT*,'offsets:',offsets

!  Now define structured type and commit it
call MPI_TYPE_CREATE_STRUCT( NFIELD, blockcounts, offsets, oldtypes, objtype,ierr)
call MPI_TYPE_COMMIT(objtype, ierr)

RETURN
END SUBROUTINE MAKE_MPI_STRUC_SP


SUBROUTINE ADD_REF_MDL(obj,tmppar)
!----------------------------------------------------------------
!This subroutine adds the reference model values on the perturbed 
!Then new model values are used to compute the prediction
!----------------------------------------------------------------
USE RJMCMC_COM
TYPE(objstruc)              :: obj
REAL(KIND=SP)               :: vs(obj%k+1),depth(obj%k+1)
REAL(KIND=SP)               :: vsref
REAL(KIND=SP)               :: tmppar(obj%NFP)
!----------------------------------------------------------------
!For half space
!----------------------------------------------------------------
depth = 0._SP
tmppar = 0._SP
vsref = 0._SP
tmppar(1:obj%NFP) = obj%par(1:obj%NFP)

IF(obj%k == 0)THEN
 depth(obj%k+1) = hmx        !Close to the surface reference values
 vs(obj%k+1)    = obj%par(1)
ELSE
  depth(obj%k+1)= hmx
  vs(obj%k+1)   = obj%par(obj%k*NPL+1)
  DO i1 = 1, obj%k
    depth(i1)   = obj%z(i1)
    vs(i1)      = obj%par(NPL*i1)
  ENDDO
ENDIF

IF(obj%k == 0)THEN
  vsref   = 0._SP
  vpvsref = 0._SP
  CALL GETREFVAL(vsref,depth(obj%k+1))
  tmppar(1) = vs(obj%k+1) + vsref
ELSE
  vsref   = 0._SP
  CALL GETREFVAL(vsref,depth(obj%k+1))
  tmppar(obj%k*NPL+1) = vs(obj%k+1) + vsref
  DO i1 = 1, obj%k
    vsref = 0._SP
    vpvsref = 0._SP
    CALL GETREFVAL(vsref,depth(i1))
    tmppar(NPL*i1) = vs(i1) + vsref
  ENDDO
ENDIF 
END SUBROUTINE ADD_REF_MDL

SUBROUTINE GETREFVAL(vsref,z)
!-----------------------------------------------------------------------
!For a given thickness, this subroutine computes the S-velocity and Vp/Vs
!ratio by interpolating the given model
!-----------------------------------------------------------------------
USE RJMCMC_COM
IMPLICIT NONE
REAL(KIND=SP)   :: z,vsref,grad,dz
INTEGER(KIND=IB):: ipar
vsref = 0.0
grad = 0.0
dz = 0.0
IF(z >= vel_ak135(nreflyr,1))THEN
  vsref = vel_ak135(nreflyr,2)
ELSEIF(z == 0.)THEN
  vsref = vel_ak135(1,2)
ELSE
  DO ipar=1,nreflyr
    IF((z == vel_ak135(ipar,1)))THEN
      vsref = vel_ak135(ipar,2)
    ELSEIF(z > vel_ak135(ipar,1))THEN
!Linear interpolation of velocity and Vp/Vs ratio 
      grad = (vel_ak135(ipar+1,2)-vel_ak135(ipar,2))/(vel_ak135(ipar+1,1)-vel_ak135(ipar,1))
      dz = (z-vel_ak135(ipar,1))
      vsref = vel_ak135(ipar,2) + dz*grad
    ENDIF
 ENDDO
ENDIF
RETURN
END SUBROUTINE GETREFVAL

!-----------------------------------------------------------------------------
!SUBROUTINE TO GENERATE THE GAUSSIAN RANDOM NOISE OF STANDARD DEVIATION 1 AND 
!MEAN ZERO USING GASSDEV FUNCTION FROM FORTRAN
!-----------------------------------------------------------------------------
FUNCTION gasdev(idum)
INTEGER idum
REAL gasdev
!CU    USES ran1
INTEGER iset
REAL fac,gset,rsq,v1,v2,ran1
SAVE iset,gset
DATA iset/0/
if (iset.eq.0) then
1   v1=2.*ran1(idum)-1.
    v2=2.*ran1(idum)-1.
    rsq=v1**2+v2**2
    if(rsq.ge.1..or.rsq.eq.0.)goto 1
    fac=sqrt(-2.*log(rsq)/rsq)
    gset=v1*fac
    gasdev=v2*fac
    iset=1
else
    gasdev=gset
    iset=0
endif
return
END
!C  (C) Copr. 1986-92 Numerical Recipes Software '%1&9p#!.
FUNCTION ran1(idum)
INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV
REAL ran1,AM,EPS,RNMX
PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836, &
           NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
INTEGER j,k,iv(NTAB),iy
SAVE iv,iy
DATA iv /NTAB*0/, iy /0/
if (idum.le.0.or.iy.eq.0) then
  idum=max(-idum,1)
  do 11 j=NTAB+8,1,-1
      k=idum/IQ
      idum=IA*(idum-k*IQ)-IR*k
      if (idum.lt.0) idum=idum+IM
      if (j.le.NTAB) iv(j)=idum
11 continue
    iy=iv(1)
endif

k=idum/IQ
idum=IA*(idum-k*IQ)-IR*k
if (idum.lt.0) idum=idum+IM
j=1+iy/NDIV
iy=iv(j)
iv(j)=idum
ran1=min(AM*iy,RNMX)
return
END
!C  (C) Copr. 1986-92 Numerical Recipes Software '%1&9p#!.
!===========================================================
! This is the end my fiend...
! EOF
