diff --git a/CMakeLists.txt b/CMakeLists.txt
index 373559f..c2678c9 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.10)
 
 # set the project name
 set(CMAKE_PROJECT_NAME "SurfATT")
-project(${CMAKE_PROJECT_NAME} VERSION 1.5.1 LANGUAGES Fortran CXX)
+project(${CMAKE_PROJECT_NAME} VERSION 1.6.0 LANGUAGES Fortran CXX)
 
 # set the version number and git commit hash
 execute_process(
diff --git a/examples/00_checkerboard_iso/input_params.yml b/examples/00_checkerboard_iso/input_params.yml
index 867ddb8..fe8dd21 100644
--- a/examples/00_checkerboard_iso/input_params.yml
+++ b/examples/00_checkerboard_iso/input_params.yml
@@ -28,6 +28,9 @@ topo:
   wavelen_factor: 2.5
 
 inversion:
+  # ------- Inversion parameters --------
+  use_alpha_beta_rho: false
+  rho_scaling: false
   # ------- Initial model --------
   init_model_type: 0 # 0: increase from v0 to v1; 1: 1D inversion for average tt; 2: use a 3D vs model
   vel_range: [1.8, 4.2]
diff --git a/examples/04_hawaii_topo/input_params.yml b/examples/04_hawaii_topo/input_params.yml
index a3d998f..d33949c 100644
--- a/examples/04_hawaii_topo/input_params.yml
+++ b/examples/04_hawaii_topo/input_params.yml
@@ -22,6 +22,9 @@ topo:
   wavelen_factor: 2.5
 
 inversion:
+  # ------- Inversion parameters --------
+  use_alpha_beta_rho: true
+  rho_scaling: true
   # ------- Initial model --------
   init_model_type: 1 # 0: increase from v0 to v1; 1: 1D inversion for average tt; 2: use a 3D vs model
   vel_range: [1.8, 4.2]
diff --git a/examples/04_hawaii_topo/run_this_example.sh b/examples/04_hawaii_topo/run_this_example.sh
index 0861a2e..3e76b4f 100644
--- a/examples/04_hawaii_topo/run_this_example.sh
+++ b/examples/04_hawaii_topo/run_this_example.sh
@@ -6,7 +6,7 @@ pos_str=19.5/-155.5
 angle=-30
 
 # download topography
-gmt grdcut @earth_relief_01m -R-157/-152/18/21 -Ghawaii.nc
+# gmt grdcut @earth_relief_01m -R-157/-152/18/21 -Ghawaii.nc
 
 # rotate source receiver file
 ../../bin/surfatt_rotate_src_rec -i src_rec_file_raw.csv -a $angle -c $pos_str -o src_rec_file_rotated.csv
diff --git a/src/acqui.f90 b/src/acqui.f90
index 4030fc5..acaf690 100644
--- a/src/acqui.f90
+++ b/src/acqui.f90
@@ -32,7 +32,7 @@ module acqui
     real(kind=dp), dimension(:,:,:), pointer               :: adj_s, adj_density
     real(kind=dp), dimension(:,:,:), allocatable           :: adj_s_local, adj_density_local
     real(kind=dp), dimension(:,:,:,:), pointer             :: sen_vsRc, sen_vpRc, sen_rhoRc
-    real(kind=dp), dimension(:,:,:), pointer               :: ker_beta, ker_density
+    real(kind=dp), dimension(:,:,:,:), pointer             :: ker, ker_density
     real(kind=dp)                                          :: chi0
     contains
     procedure :: init => att_acqui_init, post_proc => post_processing_for_kernel
@@ -43,7 +43,7 @@ module acqui
   type(att_acqui), target, public                          :: att_acqui_global_ph, att_acqui_global_gr
   integer :: win_adj_s,win_adj_density,win_isrcs,&
             win_sen_vsRc, win_sen_vpRc, win_sen_rhoRc,&
-            win_ker_beta, win_ker_density
+            win_ker,win_ker_density
 contains
   subroutine att_acqui_init(this, itype, is_fwd)
     class(att_acqui), intent(inout) :: this
@@ -82,7 +82,7 @@ subroutine prepare_inv(this)
       this%sen_vsRc = 0._dp
       this%sen_vpRc = 0._dp
       this%sen_rhoRc = 0._dp
-      this%ker_beta = 0._dp
+      this%ker = 0._dp
       this%ker_density = 0._dp
     endif
     call synchronize_all()
@@ -206,23 +206,42 @@ subroutine combine_kernels(this)
     call sum_all(this%adj_density_local, this%adj_density, this%sr%nperiod, am%n_xyz(1), am%n_xyz(2))
     call write_log('Combining eikonal and surface wave kernels...',1,this%module)
     if (myrank == 0) then
-      do ip = 1, this%sr%nperiod
-        do i = 1, am%n_xyz(3)
-          this%ker_beta(:,:,i) = this%ker_beta(:,:,i) -this%adj_s(ip,:,:) * this%sen_vsRc(ip,:,:,i) &
-                               - this%adj_s(ip,:,:) * this%sen_vpRc(ip,:,:,i) * dalpha_dbeta(am%vs3d_opt(:,:,i)) &
-                               - this%adj_s(ip,:,:) * this%sen_rhoRc(ip,:,:,i) * &
-                                 drho_dalpha(am%vp3d_opt(:,:,i)) * dalpha_dbeta(am%vs3d_opt(:,:,i))
-          this%ker_density(:,:,i) = this%ker_density(:,:,i) - this%adj_density(ip,:,:) * this%sen_vsRc(ip,:,:,i) &
-                                  - this%adj_density(ip,:,:) * this%sen_vpRc(ip,:,:,i) * dalpha_dbeta(am%vs3d_opt(:,:,i)) &
-                                  - this%adj_density(ip,:,:) * this%sen_rhoRc(ip,:,:,i) * &
-                                    drho_dalpha(am%vp3d_opt(:,:,i)) * dalpha_dbeta(am%vs3d_opt(:,:,i))
+      if (ap%inversion%use_alpha_beta_rho) then
+        do ip = 1, this%sr%nperiod
+          do i = 1, am%n_xyz(3)
+            this%ker(1,:,:,i) = this%ker(1,:,:,i) - this%adj_s(ip,:,:) * this%sen_vsRc(ip,:,:,i)
+            this%ker_density(1,:,:,i) = this%ker_density(1,:,:,i) - this%adj_density(ip,:,:) * this%sen_vsRc(ip,:,:,i)
+            this%ker(2,:,:,i) = this%ker(2,:,:,i) - this%adj_s(ip,:,:) * this%sen_vpRc(ip,:,:,i)
+            this%ker_density(2,:,:,i) = this%ker_density(1,:,:,i) - this%adj_density(ip,:,:) * this%sen_vpRc(ip,:,:,i)
+            if (.not. ap%inversion%rho_scaling) then
+              this%ker(3,:,:,i) = this%ker(3,:,:,i) - this%adj_s(ip,:,:) * this%sen_rhoRc(ip,:,:,i)
+              this%ker_density(3,:,:,i) = this%ker_density(3,:,:,i) - this%adj_density(ip,:,:) * this%sen_rhoRc(ip,:,:,i)
+            endif
+          enddo
         enddo
-      enddo
-      this%ker_beta = this%ker_beta / this%chi0
+        if (ap%inversion%rho_scaling) then
+          this%ker(3,:,:,:) = RHO_SCALING *this%ker(1,:,:,:)
+          this%ker_density(3,:,:,:) = this%ker_density(1,:,:,:)
+        endif
+      else
+        do ip = 1, this%sr%nperiod
+          do i = 1, am%n_xyz(3)
+            this%ker(1,:,:,i) = this%ker(1,:,:,i) -this%adj_s(ip,:,:) * this%sen_vsRc(ip,:,:,i) &
+                                - this%adj_s(ip,:,:) * this%sen_vpRc(ip,:,:,i) * dalpha_dbeta(am%vs3d_opt(:,:,i)) &
+                                - this%adj_s(ip,:,:) * this%sen_rhoRc(ip,:,:,i) * &
+                                  drho_dalpha(am%vp3d_opt(:,:,i)) * dalpha_dbeta(am%vs3d_opt(:,:,i))
+            this%ker_density(1,:,:,i) = this%ker_density(1,:,:,i) - this%adj_density(ip,:,:) * this%sen_vsRc(ip,:,:,i) &
+                                    - this%adj_density(ip,:,:) * this%sen_vpRc(ip,:,:,i) * dalpha_dbeta(am%vs3d_opt(:,:,i)) &
+                                    - this%adj_density(ip,:,:) * this%sen_rhoRc(ip,:,:,i) * &
+                                      drho_dalpha(am%vp3d_opt(:,:,i)) * dalpha_dbeta(am%vs3d_opt(:,:,i))
+          enddo
+        enddo
+      endif
+      ! this%ker_beta = this%ker_beta / this%chi0
     endif
     call synchronize_all()
-    call sync_from_main_rank(this%ker_beta, am%n_xyz(1), am%n_xyz(2), am%n_xyz(3))
-    call sync_from_main_rank(this%ker_density, am%n_xyz(1), am%n_xyz(2), am%n_xyz(3))
+    call sync_from_main_rank(this%ker, nker, am%n_xyz(1), am%n_xyz(2), am%n_xyz(3))
+    call sync_from_main_rank(this%ker_density, nker, am%n_xyz(1), am%n_xyz(2), am%n_xyz(3))
   end subroutine combine_kernels
 
   subroutine post_processing_for_kernel(this)
@@ -230,10 +249,10 @@ subroutine post_processing_for_kernel(this)
     !! FILEPATH: /Users/xumijian/Codes/SurfATT/src/tomo.f90
     !! @param[inout] this an object of type att_acqui
     class(att_acqui), intent(inout) :: this
-    integer :: nxinv,nyinv,nzinv,nset
+    integer :: nxinv,nyinv,nzinv,nset,iker
     real(kind=dp) :: step, gkmax
     real(kind=dp), dimension(:), allocatable :: gk,gk_precond
-    real(kind=dp), dimension(:,:,:), allocatable :: update, precond
+    real(kind=dp), dimension(:,:,:), allocatable :: update, precond, this_ker
     character(len=MAX_STRING_LEN) :: fname
 
     if (myrank == 0) then
@@ -241,43 +260,43 @@ subroutine post_processing_for_kernel(this)
       nxinv = ap%inversion%n_inv_grid(1)
       nyinv = ap%inversion%n_inv_grid(2)
       nzinv = ap%inversion%n_inv_grid(3)
-      ! initial inversion grid kernel
       call write_log('This is post processing for '//trim(this%gr_name)//' kernel...',1,this%module)
-      gk = zeros(nset*nxinv*nyinv*nzinv)
-      ! gk_precond = zeros(nset*nxinv*nyinv*nzinv)
-      update = zeros(am%n_xyz(1),am%n_xyz(2),am%n_xyz(3))
-      ! to inversion grids
-      if (ap%inversion%kdensity_coe > 0) then
-        call this%regularize_ker_density(precond)
-        this%ker_beta = this%ker_beta*precond
-      endif
-      call inv_grid_iso(am%xinv,am%yinv,am%zinv, nxinv, nyinv, nzinv,&
-                        nset,this%ker_beta,am%n_xyz(1),am%n_xyz(2),am%n_xyz(3),&
-                        gk,am%xgrids,am%ygrids,am%zgrids)
-      ! call inv_grid_iso(am%xinv,am%yinv,am%zinv, nxinv, nyinv, nzinv,&
-      !                   nset,precond,am%n_xyz(1),am%n_xyz(2),am%n_xyz(3),&
-      !                   gk_precond,am%xgrids,am%ygrids,am%zgrids)
-      call inv2fwd_iso(am%xinv,am%yinv,am%zinv,nxinv,nyinv,nzinv,&
-                       ap%inversion%ncomponents,am%n_xyz(1),am%n_xyz(2),am%n_xyz(3), &
-                       gk,am%xgrids,am%ygrids,am%zgrids,update)
-      update = update / nset
-      gradient_s = gradient_s + update*ap%data%weights(this%itype)
+      do iker = 1, nker
+        ! initial inversion grid kernel
+        gk = zeros(nset*nxinv*nyinv*nzinv)
+        ! gk_precond = zeros(nset*nxinv*nyinv*nzinv)
+        update = zeros(am%n_xyz(1),am%n_xyz(2),am%n_xyz(3))
+        ! to inversion grids
+        if (ap%inversion%kdensity_coe > 0) then
+          call this%regularize_ker_density(iker, precond)
+          this%ker(iker,:,:,:) = this%ker(iker,:,:,:)*precond
+        endif
+        call inv_grid_iso(am%xinv,am%yinv,am%zinv, nxinv, nyinv, nzinv,&
+                          nset,this%ker(iker,:,:,:),am%n_xyz(1),am%n_xyz(2),am%n_xyz(3),&
+                          gk,am%xgrids,am%ygrids,am%zgrids)
+        call inv2fwd_iso(am%xinv,am%yinv,am%zinv,nxinv,nyinv,nzinv,&
+                        ap%inversion%ncomponents,am%n_xyz(1),am%n_xyz(2),am%n_xyz(3), &
+                        gk,am%xgrids,am%ygrids,am%zgrids,update)
+        update = update / nset / this%chi0
+        gradient_s(iker,:,:,:) = gradient_s(iker,:,:,:) + update*ap%data%weights(this%itype)
+      enddo
     endif
     call synchronize_all()
 
   end subroutine post_processing_for_kernel
 
-  subroutine regularize_ker_density(this, precond)
+  subroutine regularize_ker_density(this, kernel_id, precond)
     !> Regularize the kernel density.
     !! Inputs:
     !!   this: an object of type att_acqui
     !! Outputs:
     !!   precond: a 3D array of real numbers representing the preconditioned kernel density
     class(att_acqui), intent(inout) :: this
+    integer :: kernel_id
     real(kind=dp), dimension(:,:,:), allocatable, intent(out) :: precond
 
     precond = zeros(am%n_xyz(1), am%n_xyz(2), am%n_xyz(3))
-    precond = abs(this%ker_density)/maxval(abs(this%ker_density))
+    precond = abs(this%ker_density(kernel_id,:,:,:))/maxval(abs(this%ker_density(kernel_id,:,:,:)))
     where(precond<precond_thres)
       precond = 0
     elsewhere
@@ -307,8 +326,8 @@ subroutine allocate_shm_arrays(this)
       call prepare_shm_array_dp_4d(this%sen_vsRc, this%sr%nperiod, am%n_xyz(1), am%n_xyz(2), am%n_xyz(3), win_sen_vsRc)
       call prepare_shm_array_dp_4d(this%sen_vpRc, this%sr%nperiod, am%n_xyz(1), am%n_xyz(2), am%n_xyz(3), win_sen_vpRc)
       call prepare_shm_array_dp_4d(this%sen_rhoRc, this%sr%nperiod, am%n_xyz(1), am%n_xyz(2), am%n_xyz(3), win_sen_rhoRc)
-      call prepare_shm_array_dp_3d(this%ker_beta, am%n_xyz(1), am%n_xyz(2), am%n_xyz(3), win_ker_beta)
-      call prepare_shm_array_dp_3d(this%ker_density, am%n_xyz(1), am%n_xyz(2), am%n_xyz(3), win_ker_density)
+      call prepare_shm_array_dp_4d(this%ker, nker, am%n_xyz(1), am%n_xyz(2), am%n_xyz(3), win_ker)
+      call prepare_shm_array_dp_4d(this%ker_density, nker, am%n_xyz(1), am%n_xyz(2), am%n_xyz(3), win_ker_density)
     endif
     call synchronize_all()
   end subroutine allocate_shm_arrays
diff --git a/src/model.f90 b/src/model.f90
index 0353e62..0fedcc6 100644
--- a/src/model.f90
+++ b/src/model.f90
@@ -236,14 +236,16 @@ subroutine inv1d(this, vsinv, niter, misfits)
     niter = iter
   end subroutine inv1d
 
-  subroutine add_pert(this,nx,ny,nz,pert_vel,hmarg,anom_size)
+  subroutine add_pert(this,nx,ny,nz,pert_vel,hmarg,anom_size,only_vs)
     class(att_model), intent(inout) :: this
     integer,  intent(in) :: nx, ny, nz
     real(kind=dp), intent(in), optional :: pert_vel,hmarg,anom_size
+    logical, intent(in), optional :: only_vs
     real(kind=dp) :: ashmarg,aspert_vel,amp, asanom_size
     real(kind=dp), dimension(:), allocatable :: x_pert, y_pert, z_pert
     real(kind=dp), dimension(:,:,:), allocatable :: vs_pert
     real(kind=dp), dimension(3) :: para
+    logical :: asonly_vs
     integer :: ntaperx, ntapery, i, j, k
 
     if (present(pert_vel)) then 
@@ -264,6 +266,12 @@ subroutine add_pert(this,nx,ny,nz,pert_vel,hmarg,anom_size)
       asanom_size = 0.
     endif
 
+    if (present(only_vs)) then
+      asonly_vs = only_vs
+    else
+      asonly_vs = .false.
+    endif
+
     if (myrank == 0) then
       ntaperx = ashmarg/this%d_xyz(1)
       ntapery = ashmarg/this%d_xyz(2)
@@ -289,8 +297,10 @@ subroutine add_pert(this,nx,ny,nz,pert_vel,hmarg,anom_size)
         enddo
       enddo
       this%vs3d = this%vs3d*(1+vs_pert)
-      this%vp3d = empirical_vp(this%vs3d)
-      this%rho3d = empirical_rho(this%vp3d)
+      if (.not. asonly_vs) then
+        this%vp3d = empirical_vp(this%vs3d)
+        this%rho3d = empirical_rho(this%vp3d)
+      endif
     endif
     call synchronize_all()
     call sync_from_main_rank(this%vs3d, this%n_xyz(1), this%n_xyz(2), this%n_xyz(3))
@@ -314,6 +324,10 @@ subroutine write_model(this, subname)
       call h%add('/lat',this%ygrids)
       call h%add('/dep',this%zgrids)
       call h%add('/vs',transpose_3(this%vs3d))
+      if (ap%inversion%use_alpha_beta_rho) then
+        call h%add('/vp',transpose_3(this%vp3d))
+        call h%add('/rho',transpose_3(this%rho3d))
+      endif
       call h%close()
     endif
     call synchronize_all()    
diff --git a/src/optimize.f90 b/src/optimize.f90
index 9255104..777dc38 100644
--- a/src/optimize.f90
+++ b/src/optimize.f90
@@ -8,6 +8,9 @@
 !                           (c) October 2023
 !   
 !     Changing History: Apl 2024, Initialize Codes
+!                       May 2024, Add L-BFGS and CG optimization
+!                       Jun 2024, Add Hager-Zhang formula for CG
+!                       Jan 2025, change dims for Alpha_beta_rho
 !
 !=====================================================================
 module optimize
@@ -20,13 +23,13 @@ module optimize
 contains
   subroutine get_lbfgs_direction(iter, direction)
     integer, intent(in) :: iter
-    real(kind=dp), dimension(:,:,:), allocatable, intent(out) :: direction
-    real(kind=dp), dimension(:,:,:), allocatable :: gradient0,gradient1,model0,model1,&
+    real(kind=dp), dimension(:,:,:,:), allocatable, intent(out) :: direction
+    real(kind=dp), dimension(:,:,:,:), allocatable :: gradient0,gradient1,model0,model1,&
                                                     q_vector,r_vector
-    real(kind=dp), dimension(:,:,:,:), allocatable :: gradient_diff, model_diff
+    real(kind=dp), dimension(:,:,:,:,:), allocatable :: gradient_diff, model_diff
     real(kind=dp), dimension(:), allocatable :: p, a
     real(kind=dp) :: p_k_up_sum, p_k_down_sum, p_k, b
-    integer :: istore, i, dims(3), iter_store, nstore
+    integer :: istore, i, dims(4), iter_store, nstore
     integer, dimension(:), allocatable :: idx_iter
 
     iter_store = iter-m_store
@@ -35,8 +38,8 @@ subroutine get_lbfgs_direction(iter, direction)
 
     call get_gradient(iter, q_vector)
     dims = shape(q_vector)
-    allocate(gradient_diff(nstore, dims(1), dims(2), dims(3)))
-    allocate(model_diff(nstore, dims(1), dims(2), dims(3)))
+    allocate(gradient_diff(nstore, nker, dims(2), dims(3), dims(4)))
+    allocate(model_diff(nstore, nker, dims(2), dims(3), dims(4)))
     idx_iter = zeros(nstore)
     p = zeros(nstore)
     a = zeros(nstore)
@@ -49,22 +52,22 @@ subroutine get_lbfgs_direction(iter, direction)
       call get_model(istore, model0)
       call get_model(istore+1, model1)
 
-      model_diff(i,:,:,:) = model1 - model0
-      gradient_diff(i,:,:,:) = gradient1 - gradient0
+      model_diff(i,:,:,:,:) = model1 - model0
+      gradient_diff(i,:,:,:,:) = gradient1 - gradient0
 
-      p(i) = 1/sum(model_diff(i,:,:,:)*gradient_diff(i,:,:,:))
-      a(i) = sum(model_diff(i,:,:,:)*q_vector)*p(i)
-      q_vector = q_vector - a(i)*gradient_diff(i,:,:,:)
+      p(i) = 1/sum(model_diff(i,:,:,:,:)*gradient_diff(i,:,:,:,:))
+      a(i) = sum(model_diff(i,:,:,:,:)*q_vector)*p(i)
+      q_vector = q_vector - a(i)*gradient_diff(i,:,:,:,:)
     enddo
-    p_k_up_sum = sum(model_diff(1,:,:,:)*gradient_diff(1,:,:,:))
-    p_k_down_sum = sum(gradient_diff(1,:,:,:)*gradient_diff(1,:,:,:))
+    p_k_up_sum = sum(model_diff(1,:,:,:,:)*gradient_diff(1,:,:,:,:))
+    p_k_down_sum = sum(gradient_diff(1,:,:,:,:)*gradient_diff(1,:,:,:,:))
     p_k = p_k_up_sum/p_k_down_sum
     r_vector = p_k*q_vector
 
     do istore = iter_store,iter-1
       i = find_loc(idx_iter, istore)
-      b = sum(gradient_diff(i,:,:,:)*r_vector)*p(i)
-      r_vector = r_vector + model_diff(i,:,:,:)*(a(i)-b)
+      b = sum(gradient_diff(i,:,:,:,:)*r_vector)*p(i)
+      r_vector = r_vector + model_diff(i,:,:,:,:)*(a(i)-b)
     enddo
     direction = -1.0_dp * r_vector
 
@@ -73,8 +76,8 @@ end subroutine get_lbfgs_direction
   subroutine get_cg_direction(iter, direction)
     ! get the conjugate gradient direction using Hager-Zhang formula
     integer, intent(in) :: iter
-    real(kind=dp), dimension(:,:,:), allocatable, intent(out) :: direction
-    real(kind=dp), dimension(:,:,:), allocatable :: gradient0,gradient1,direction0,grad_diff
+    real(kind=dp), dimension(:,:,:,:), allocatable, intent(out) :: direction
+    real(kind=dp), dimension(:,:,:,:), allocatable :: gradient0,gradient1,direction0,grad_diff
     real(kind=dp) :: beta
 
     call get_gradient(iter-1, gradient0)
@@ -90,33 +93,50 @@ end subroutine get_cg_direction
 
   subroutine get_model(iter, model)
     integer, intent(in) :: iter
-    real(kind=dp), dimension(:,:,:), allocatable, intent(out) :: model
+    real(kind=dp), dimension(:,:,:,:), allocatable, intent(out) :: model
+    real(kind=dp), dimension(:,:,:), allocatable :: tmp_model
     character(len=MAX_NAME_LEN) :: key_name
 
     write(key_name, '("/vs_",I3.3)') iter
-    call h5read(model_fname, key_name, model)
-    model = transpose_3(model)
+    call h5read(model_fname, key_name, tmp_model)
+    tmp_model = transpose_3(tmp_model)
+    model = zeros( nker,size(tmp_model, 1), size(tmp_model, 2), size(tmp_model, 3))
+    model(1,:,:,:) = tmp_model
+    if (ap%inversion%use_alpha_beta_rho) then
+      ! read vp model
+      deallocate(tmp_model)
+      write(key_name, '("/vp_",I3.3)') iter
+      call h5read(model_fname, key_name, tmp_model)
+      model(2,:,:,:) = tmp_model
+
+      ! read rho model
+      deallocate(tmp_model)
+      write(key_name, '("/rho_",I3.3)') iter
+      call h5read(model_fname, key_name, tmp_model)
+      model(3,:,:,:) = tmp_model
+    endif
+
   end subroutine get_model
 
   subroutine get_gradient(iter, gradient)
     integer, intent(in) :: iter
-    real(kind=dp), dimension(:,:,:), allocatable, intent(out) :: gradient
+    real(kind=dp), dimension(:,:,:,:), allocatable, intent(out) :: gradient
     character(len=MAX_NAME_LEN) :: key_name
 
     write(key_name, '("/gradient_",I3.3)') iter
     call h5read(model_fname, key_name, gradient)
-    gradient = transpose_3(gradient)
+    gradient = transpose_4(gradient)
 
   end subroutine get_gradient
 
   subroutine get_direction(iter, direction)
     integer, intent(in) :: iter
-    real(kind=dp), dimension(:,:,:), allocatable, intent(out) :: direction
+    real(kind=dp), dimension(:,:,:,:), allocatable, intent(out) :: direction
     character(len=MAX_NAME_LEN) :: key_name
 
     write(key_name, '("/direction_",I3.3)') iter
     call h5read(model_fname, key_name, direction)
-    direction = transpose_3(direction)
+    direction = transpose_4(direction)
 
   end subroutine
 end module optimize
diff --git a/src/para.f90 b/src/para.f90
index b5fc9b9..99b73cb 100644
--- a/src/para.f90
+++ b/src/para.f90
@@ -55,6 +55,7 @@ module para
     character(len=MAX_NAME_LEN), dimension(2)              :: optim_name = ['CG   ', 'LBFGS']
     real(kind=dp)                                          :: vel_range(2), step_length, min_derr=0.001,&
                                                               kdensity_coe=0.5, maxshrink=0.8, sigma_2d=0.
+    logical                                                :: use_alpha_beta_rho, rho_scaling
   end type para_inversion
 
   type, public :: att_para
@@ -192,7 +193,13 @@ subroutine read_par_file(this, fname)
         this%inversion%kdensity_coe = inversion%get_real('kdensity_coe',default=0.0,error=io_err)
         this%inversion%sigma_2d = inversion%get_real('sigma_2d',default=0.0,error=io_err)
         m_store = inversion%get_integer('lbfgs_model_store',default=5,error=io_err)
-
+        this%inversion%use_alpha_beta_rho = inversion%get_logical( &
+          'use_alpha_beta_rho', default=.false., error=io_err &
+        )
+        if (this%inversion%use_alpha_beta_rho) nker = 3
+        this%inversion%rho_scaling = inversion%get_logical( &
+          'rho_scaling', default=.true., error=io_err &
+        )
       end select
       call root%finalize()
       deallocate(root)
@@ -238,6 +245,9 @@ subroutine read_par_file(this, fname)
     call bcast_all(this%inversion%optim_method)
     call bcast_all(this%inversion%sigma_2d)
     call bcast_all(m_store)
+    call bcast_all(this%inversion%use_alpha_beta_rho)
+    call bcast_all(this%inversion%rho_scaling)
+    call bcast_all(nker)
     
     call synchronize_all()
 
diff --git a/src/shared/argparse.f90 b/src/shared/argparse.f90
index ad2e842..c768e5c 100644
--- a/src/shared/argparse.f90
+++ b/src/shared/argparse.f90
@@ -61,18 +61,20 @@ subroutine argparse_tomo(fname, isfwd)
     endif
   end subroutine argparse_tomo
 
-  subroutine argparse_cb_fwd(fname, ncb, pert_vel, hmarg, anom_size, max_noise)
+  subroutine argparse_cb_fwd(fname, ncb, pert_vel, hmarg, anom_size, max_noise, only_vs)
     character(len=MAX_STRING_LEN),dimension(:), allocatable :: args
     character(len=MAX_STRING_LEN) :: arg, value
     character(len=MAX_STRING_LEN),intent(out) :: fname
     integer, dimension(3), intent(out) :: ncb
     real(kind=dp), intent(out) :: pert_vel, hmarg, anom_size, max_noise
-    integer :: i,nargs,m,ier, nopt=6
+    logical, intent(out) :: only_vs
+    integer :: i,nargs,m,ier, nopt=7
 
     pert_vel = 0.08
     hmarg = 0.
     anom_size = 0.
     max_noise = 0.
+    only_vs = .false.
     nargs = command_argument_count()
     allocate(args(nargs))
     do i = 1, nargs
@@ -91,6 +93,7 @@ subroutine argparse_cb_fwd(fname, ncb, pert_vel, hmarg, anom_size, max_noise)
       write(*,*)'optional arguments:'
       write(*,*)' -h                  Print help message'
       write(*,*)' -e tt_noise         Add random noise to travel time data, defaults to 0'
+      write(*,*)' -v                  Only add perturbations on Vs model'
       write(*,*)' -m margin_degree    Margin area in degree between anomaly as boundary, defaults to 0'
       write(*,*)' -p pert_vel         Magnitude of velocity perturbations, defaults to 0.08'
       write(*,*)' -s anom_size_km     size of anomalies at the top in km, default to uniform anomaly size in Z direction'
@@ -121,7 +124,10 @@ subroutine argparse_cb_fwd(fname, ncb, pert_vel, hmarg, anom_size, max_noise)
       elseif(arg(1:2) == '-e') then
         m = m+1
         read(args(i+1),*,iostat=ier) max_noise
-        if(ier/=0) stop 'Cannot parse argument -e'         
+        if(ier/=0) stop 'Cannot parse argument -e'
+      elseif(arg(1:2) == '-v') then
+        m = m+1
+        only_vs = .true.      
       endif
     enddo
     if (m<1) then
diff --git a/src/shared/hdf5_interface.f90 b/src/shared/hdf5_interface.f90
index 0788f26..97d66fe 100644
--- a/src/shared/hdf5_interface.f90
+++ b/src/shared/hdf5_interface.f90
@@ -6,12 +6,12 @@ module hdf5_interface
   public :: hdf5_file
 
   interface h5read
-    module procedure h5read_f_1d, h5read_f_2d, h5read_f_3d
+    module procedure h5read_f_1d, h5read_f_2d, h5read_f_3d, h5read_f_4d
   end interface h5read
   public :: h5read
 
   interface h5write
-    module procedure h5write_f_1d, h5write_f_2d, h5write_f_3d
+    module procedure h5write_f_1d, h5write_f_2d, h5write_f_3d, h5write_f_4d
   end interface h5write
   public :: h5write
 
@@ -52,6 +52,7 @@ module hdf5_interface
                         hdf_add_real1d,&
                         hdf_add_real2d,&
                         hdf_add_real3d,&
+                        hdf_add_real4d,&
                         hdf_add_string
 
     !> get dataset integer/real 0-3d
@@ -63,6 +64,7 @@ module hdf5_interface
                         hdf_get_real1d,&
                         hdf_get_real2d,&
                         hdf_get_real3d,&
+                        hdf_get_real4d,&
                         hdf_get_string
 
     !> add attribute
@@ -85,7 +87,9 @@ module hdf5_interface
     procedure,private :: hdf_add_real2d
     procedure,private :: hdf_get_real2d
     procedure,private :: hdf_add_real3d
+    procedure,private :: hdf_add_real4d
     procedure,private :: hdf_get_real3d
+    procedure,private :: hdf_get_real4d
     procedure,private :: hdf_get_string
     procedure,private :: hdf_add_string
     procedure,private :: hdf_adda_string
@@ -396,6 +400,21 @@ subroutine hdf_add_real3d(self,dname,value)
       rank(value), int(shape(value),HSIZE_T), h5kind_to_type(kind(value),H5_REAL_KIND), value, ierr)
 
   end subroutine hdf_add_real3d
+
+  !=============================================================================
+  subroutine hdf_add_real4d(self,dname,value)
+    class(hdf5_file), intent(in) :: self
+    character(*), intent(in) :: dname
+    real(kind=dp), intent(in)      :: value(:,:,:,:)
+
+    integer         :: ierr
+
+    call self%add(dname)
+
+    call h5ltmake_dataset_f(self%lid, dname, &
+      rank(value), int(shape(value),HSIZE_T), h5kind_to_type(kind(value),H5_REAL_KIND), value, ierr)
+
+  end subroutine hdf_add_real4d
   !=============================================================================
   subroutine hdf_add_string(self,dname,value)
     class(hdf5_file), intent(in) :: self
@@ -580,6 +599,25 @@ subroutine hdf_get_real3d(self, dname, value)
 
   end subroutine hdf_get_real3d
 
+  subroutine hdf_get_real4d(self, dname, value)
+
+    class(hdf5_file), intent(in)     :: self
+    character(*), intent(in)         :: dname
+    real(kind=dp), intent(out),allocatable :: value(:,:,:,:)
+
+    integer(HSIZE_T) :: dims(4)
+    integer(SIZE_T)  :: dsize
+    integer :: ierr, dtype
+
+    call h5ltget_dataset_info_f(self%lid, dname, dims, dtype, dsize, ierr)
+
+    allocate(value(dims(1),dims(2),dims(3),dims(4)))
+
+    call h5ltread_dataset_f(self%lid, dname, &
+         & h5kind_to_type(kind(value),H5_REAL_KIND), value, dims,  ierr)
+
+  end subroutine hdf_get_real4d
+
 !----- Helper functions
 
   elemental function to_lower(str)
@@ -642,6 +680,20 @@ subroutine h5read_f_3d(fname, dname, value)
     
   end subroutine h5read_f_3d
 
+  subroutine h5read_f_4d(fname, dname, value)
+    character(*), intent(in) :: fname
+    character(*), intent(in) :: dname
+    real(kind=dp), allocatable, intent(out) :: value(:,:,:,:)
+
+    type(hdf5_file) :: hdf
+    integer :: ierr
+
+    call hdf%open(fname)
+    call hdf%get(dname, value)
+    call hdf%close()
+    
+  end subroutine h5read_f_4d
+
   subroutine h5write_f_1d(fname, dname, value)
     character(*), intent(in) :: fname
     character(*), intent(in) :: dname
@@ -684,4 +736,17 @@ subroutine h5write_f_3d(fname, dname, value)
     
   end subroutine h5write_f_3d
 
+  subroutine h5write_f_4d(fname, dname, value)
+    character(*), intent(in) :: fname
+    character(*), intent(in) :: dname
+    real(kind=dp), intent(in) :: value(:,:,:,:)
+
+    type(hdf5_file) :: hdf
+    integer :: ierr
+
+    call hdf%open(fname)
+    call hdf%add(dname, value)
+    call hdf%close()
+    
+  end subroutine h5write_f_4d
 end module hdf5_interface
diff --git a/src/shared/shared_par.f90 b/src/shared/shared_par.f90
index 1c0c88e..2dbbada 100644
--- a/src/shared/shared_par.f90
+++ b/src/shared/shared_par.f90
@@ -13,6 +13,7 @@ module constants
   real(kind=dp), parameter :: deg2rad = pi/180.0d0
   real(kind=dp), parameter :: rad2deg = 180.0d0/pi
   real(kind=dp), parameter :: km2deg = 1/(6371.0d0*pi/180.0d0)
+  real(kind=dp), parameter :: rho_scaling = 0.33
 
   character(len=MAX_STRING_LEN),parameter :: srfile = 'src_rec_iter.h5'
   character(len=MAX_STRING_LEN),parameter :: modfile = 'model_iter.h5'
@@ -33,10 +34,11 @@ module shared_par
   integer, dimension(:,:), allocatable :: rank_map
   integer :: LID
   integer :: loglevel
-  real(kind=dp), dimension(:,:,:), allocatable :: gradient_s, direction
+  real(kind=dp), dimension(:,:,:,:), allocatable :: gradient_s, direction
   character(len=MAX_STRING_LEN) :: log_fname, model_fname
   character(len=MAX_STRING_LEN) :: VERSION
   integer :: m_store = 5
+  integer :: nker = 1
 
 contains
   subroutine get_version()
diff --git a/src/shared/utils.f90 b/src/shared/utils.f90
index 422baed..1ceb4f9 100644
--- a/src/shared/utils.f90
+++ b/src/shared/utils.f90
@@ -832,6 +832,23 @@ pure function transpose_3(x) result(y)
     return
   end function transpose_3
 
+  pure function transpose_4(x) result(y)
+    real(kind = DPRE), dimension(:,:,:,:), intent(in) :: x
+    real(kind = DPRE), dimension(size(x,4), size(x,3), size(x,2), size(x,1)) :: y
+    integer(kind = IPRE) :: i, j, k, l
+
+    do i = 1, size(x,1)
+      do j = 1, size(x,2)
+        do k = 1, size(x,3)
+          do l = 1, size(x,4)
+            y(l,k,j,i) = x(i,j,k,l)
+          enddo
+        end do
+      end do
+    end do
+    return
+  end function transpose_4
+
   pure subroutine meshgrid2_ij(ax, ay, x, y)
     real(kind = DPRE), dimension(:), intent(in) :: ax, ay
     real(kind = DPRE), dimension(:,:), allocatable, intent(out) :: x, y
diff --git a/src/surfatt_cb_fwd.f90 b/src/surfatt_cb_fwd.f90
index 9150943..c6d7c47 100644
--- a/src/surfatt_cb_fwd.f90
+++ b/src/surfatt_cb_fwd.f90
@@ -28,12 +28,13 @@ program surfatt_cb_fwd
   character(len=MAX_STRING_LEN) :: fname
   integer, dimension(3) :: ncb
   real(kind=dp) :: pert, hmarg, anom_size, max_noise
+  logical :: only_vs
 
   ! initialize MPI
   call init_mpi()
 
   ! read command line arguments
-  call argparse_cb_fwd(fname, ncb, pert, hmarg, anom_size, max_noise)
+  call argparse_cb_fwd(fname, ncb, pert, hmarg, anom_size, max_noise, only_vs)
 
   ! read parameter file
   call ap%read(fname)
@@ -54,7 +55,7 @@ program surfatt_cb_fwd
 
   ! add perturbations
   call am%add_pert(ncb(1),ncb(2),ncb(3), pert_vel=pert,&
-                   hmarg=hmarg, anom_size=anom_size)
+                   hmarg=hmarg, anom_size=anom_size,only_vs=only_vs)
 
   ! initialize grid
   if (ap%data%vel_type(1)) then
diff --git a/src/tomo.f90 b/src/tomo.f90
index 4c73db1..2be3642 100644
--- a/src/tomo.f90
+++ b/src/tomo.f90
@@ -8,6 +8,7 @@
 !                           (c) October 2023
 !   
 !     Changing History: Oct 2023, Initialize Codes
+!                       Jan 2025, change dims for Alpha_beta_rho
 !
 !=====================================================================
 module tomo
@@ -84,6 +85,10 @@ subroutine initialize_inv(this)
         call h%add('/lat', am%ygrids)
         call h%add('/dep', am%zgrids)
         call h%add('/vs_000', transpose_3(am%vs3d))
+        if (ap%inversion%use_alpha_beta_rho) then
+          call h%add('/vp_000', transpose_3(am%vp3d))
+          call h%add('/rho_000', transpose_3(am%rho3d))
+        endif
         call h%close(finalize=.true.)
       endif
     endif
@@ -220,9 +225,15 @@ subroutine steepest_descent(this)
         call write_log(this%message, 1, this%module)
       endif
       gradient_s = -updatemax * gradient_s / max_gk
-      am%vs3d = am%vs3d * (1 + gradient_s)
-      am%vp3d = empirical_vp(am%vs3d)
-      am%rho3d = empirical_rho(am%vp3d)
+      if (ap%inversion%use_alpha_beta_rho) then
+        am%vs3d = am%vs3d * (1 + gradient_s(1,:,:,:))
+        am%vp3d = empirical_vp(am%vs3d)
+        am%rho3d = empirical_rho(am%vp3d)
+      else
+        am%vs3d = am%vs3d * (1 + gradient_s(1,:,:,:))
+        am%vp3d = am%vp3d * (1 + gradient_s(2,:,:,:))
+        am%rho3d = am%rho3d * (1 + gradient_s(3,:,:,:))
+      endif
       call write_tmp_model()
     endif
     call synchronize_all()
@@ -254,7 +265,7 @@ subroutine line_search(this)
       endif
       if (ap%output%verbose_level > 0) then
         write(secname,'(a,i3.3)') '/direction_',iter-1
-        call h5write(model_fname, secname, transpose_3(direction))
+        call h5write(model_fname, secname, transpose_4(direction))
       endif
     endif
     call synchronize_all()
@@ -288,14 +299,20 @@ end subroutine line_search
   subroutine prepare_fwd_linesearch()
     ! class(att_tomo), intent(inout) :: this
     real(kind=dp) :: max_gk
-    real(kind=dp), dimension(:,:,:), allocatable :: gradient_ls
+    real(kind=dp), dimension(:,:,:,:), allocatable :: gradient_ls
 
     if (myrank == 0) then
       max_gk = maxval(abs(direction))
       gradient_ls = updatemax * direction / max_gk
-      am%vs3d_opt = am%vs3d * (1 + gradient_ls)
-      am%vp3d_opt = empirical_vp(am%vs3d_opt)
-      am%rho3d_opt = empirical_rho(am%vp3d_opt)    
+      if (ap%inversion%use_alpha_beta_rho) then
+        am%vs3d_opt = am%vs3d * (1 + gradient_ls(1,:,:,:))
+        am%vp3d_opt = am%vp3d * (1 + gradient_ls(2,:,:,:))
+        am%rho3d_opt = am%rho3d * (1 + gradient_ls(3,:,:,:))
+      else
+        am%vs3d_opt = am%vs3d * (1 + gradient_ls(1,:,:,:))
+        am%vp3d_opt = empirical_vp(am%vs3d_opt)
+        am%rho3d_opt = empirical_rho(am%vp3d_opt)
+      endif
     endif
     call synchronize_all()
     call sync_from_main_rank(am%vs3d_opt, am%n_xyz(1), am%n_xyz(2), am%n_xyz(3))
@@ -331,6 +348,12 @@ subroutine write_tmp_model()
     if (ap%output%verbose_level > 0) then
       write(secname,'(a,i3.3)') '/vs_',iter
       call h5write(model_fname, secname, transpose_3(am%vs3d))
+      if (ap%inversion%use_alpha_beta_rho) then
+        write(secname,'(a,i3.3)') '/vp_',iter
+        call h5write(model_fname, secname, transpose_3(am%vp3d))
+        write(secname,'(a,i3.3)') '/rho_',iter
+        call h5write(model_fname, secname, transpose_3(am%rho3d))
+      endif
     endif
 
   end subroutine write_tmp_model
@@ -340,12 +363,12 @@ subroutine write_gradient()
 
     if (ap%output%verbose_level > 0) then
       write(secname,'(a,i3.3)') '/gradient_',iter-1  
-      call h5write(model_fname, secname, transpose_3(gradient_s))
+      call h5write(model_fname, secname, transpose_4(gradient_s))
       do itype = 1, 2
         if (.not. ap%data%vel_type(itype)) cycle
         call select_type()
         write(secname,'(a,a,"_",i3.3)') '/kdensity_',trim(ap%data%gr_name(itype)),iter-1
-        call h5write(model_fname, secname, transpose_3(aq%ker_density))
+        call h5write(model_fname, secname, transpose_4(aq%ker_density))
       enddo
     endif
 
@@ -385,7 +408,7 @@ subroutine reset_opt(this)
     class(att_tomo), intent(inout) :: this
   
     if (myrank == 0) then
-      if (.not. this%is_fwd) gradient_s = zeros(am%n_xyz(1),am%n_xyz(2),am%n_xyz(3))
+      if (.not. this%is_fwd) gradient_s = zeros(nker, am%n_xyz(1),am%n_xyz(2),am%n_xyz(3))
       am%vs3d_opt = am%vs3d
       am%vp3d_opt = am%vp3d
       am%rho3d_opt = am%rho3d