! subroutine CalDdatSigma(dall,obst,cbst,sigmaT,meandeltaT,stddeltaT)
subroutine CalDdatSigma(dall, obst, cbst, sigmaT, meandeltaT)
implicit none
! abandoned
! used in Yunnan
! INPUT
integer dall
real obst(dall),cbst(dall)
! OUTPUT
real  sigmaT(dall)
real meandeltaT,stddeltaT
! PARAMETER
integer jj,kk,ii
integer i,j,k
real twostdratio
real :: deltaT(dall)
! sigma(i)=std(t_res/t_obs)*t_obs
deltaT=0
meandeltaT=0

do i = 1,dall
    deltaT(i)=abs(cbst(i)/obst(i))
    meandeltaT=meandeltaT+deltaT(i)
enddo
meandeltaT=meandeltaT/dall
! meandeltaT=sum(deltaT(1:dall))/dall
stddeltaT=0
do i=1,dall
    stddeltaT=stddeltaT+(deltaT(i)-meandeltaT)**2
enddo
stddeltaT=sqrt(stddeltaT/dall)
! stddeltaT=sqrt(sum((deltaT(1:dall)-meandeltaT)**2))
do i=1,dall
    twostdratio=abs(deltaT(i)/(1.5*stddeltaT))
    if (twostdratio.gt.1.0)then
        sigmaT(i)=stddeltaT*obst(i)*exp(twostdratio-1)
    else
        sigmaT(i)=stddeltaT*obst(i)
    endif
enddo
end subroutine

! Residual Norm !
subroutine CalVsReslNorm(maxvp,dall,GVs,VsInv,datweight,Tdata,fwdTvs,resbst)
    use  lsmrblasInterface, only : dnrm2
    implicit none
    ! INPUT
    integer,intent(in):: dall,maxvp
    real,intent(in):: GVs(dall,maxvp)
    ! real,intent(in):: dv(maxvp),datweight(dall),Tdata(dall)
    real,intent(in):: VsInv(maxvp),datweight(dall),Tdata(dall)
    ! OUTPUT
    real:: resbst(dall)

    real:: resW2Nm,res2Nm
    real:: PreRes
    ! PARAMETER

    real fwdTgc(dall),fwdTgs(dall), fwdTvs(dall),fwdt(dall)
    real fwdTaa(dall),res(dall),resW(dall),per_res(dall)
    integer jj,kk,ii
    integer i,j,k
    !--------------------------------------------------------------------!
    ! Target:
    !  Calculate traveltime resiudal between observe traveltime
    !  and forward traveltime which is generated by inversion result
    !  using G from ref  model.
    !   r=Gm-d
    !--------------------------------------------------------------------!
    ! VsInv=dv(1:maxvp)
    !--------------------------------------------------------------------!
    ! forward residual
    fwdTvs=matmul(GVs,VsInv)
    fwdt=0
    res=0
    resW=0
    per_res=0
    resbst=0
    do i=1,dall
        resbst(i)=Tdata(i)-fwdTvs(i)
        resW(i)=resbst(i)*datweight(i) ! here, 1/sigma is for
        per_res(i)=resbst(i)/Tdata(i)
    enddo
    resW2Nm=dnrm2(dall,resW,1)
    res2Nm=dnrm2(dall,resbst,1)
    PreRes=sum(abs(per_res(1:dall)))/dall

    write(6 ,'(a,2f12.3)') '  dVs:  ||(Gm-d)||^2  and ||W(Gm-d)||^2: ', res2Nm, resW2Nm
    write(66,'(a,2f12.3)') '  dVs:  ||(Gm-d)||^2  and ||W(Gm-d)||^2: ', res2Nm, resW2Nm
    ! write(6,*) ' abs mean Res(Vs)/Res  (%)',PreRes*100
    ! write(66,*)' abs mean Res(Vs)/Res  (%)',PreRes*100
end subroutine


subroutine CalGcsReslNorm(maxvp,dall,GGc,GGs, GcInv, GsInv, datweight,Tdata,fwdTaa,resbst)
    use  lsmrblasInterface, only : dnrm2
    implicit none
    ! resW2Nm: predicated ||W(Gm-d)||2
    ! res2Nm: predicated ||(Gm-d)||2
    ! INPUT
    integer,intent(in):: dall,maxvp
    real,intent(in):: GGc(dall,maxvp),GGs(dall,maxvp)
    real,intent(in):: GcInv(maxvp), GsInv(maxvp), datweight(dall),Tdata(dall)
    ! OUTPUT
    real:: resbst(dall)
    real:: res2Nm,resW2Nm
    real:: PreRes
    ! PARAMETER
    ! real VsInv(maxvp),GcInv(maxvp),GsInv(maxvp)
    real fwdTgc(dall),fwdTgs(dall), fwdTvs(dall),fwdt(dall)
    real fwdTaa(dall),res(dall),resW(dall),per_res(dall)
    integer jj,kk,ii
    integer i,j,k
    real norm2,norm2W
    !--------------------------------------------------------------------!
    ! Target:
    !  Calculate traveltime resiudal between observe traveltime
    !  and forward traveltime which is generated by inversion result
    !  using G from ref  model.
    !   r=Gm-d
    !--------------------------------------------------------------------!
    ! GcInv=dv(1:maxvp)
    ! GsInv=dv(1+maxvp:2*maxvp)
    !--------------------------------------------------------------------!
    ! forward residual
    fwdTgc=matmul(GGc,GcInv)
    fwdTgs=matmul(GGs,GsInv)
    fwdTaa=0
    res=0
    resW=0
    per_res=0
    resbst=0
    norm2=0
    norm2W=0
    do i=1,dall
        fwdTaa(i)=fwdTgs(i)+fwdTgc(i)
        resbst(i)=Tdata(i)-fwdTaa(i)
        resW(i)=resbst(i)*datweight(i) ! here, 1/sigma is for data weight
        if (Tdata(i) .gt. 0)   per_res(i)=resbst(i)/Tdata(i)
        ! norm2=norm2+resbst(i)**2
        ! norm2W=norm2W+resW(i)**2
        ! write(6,*)'res=',resbst(i),'weight=',datweight(i),'resW=',resW(i),'norm2=',norm2,'norm2W=',norm2W
    enddo
    resW2Nm=dnrm2(dall,resW,1)
    res2Nm=dnrm2(dall,resbst,1)
    PreRes=sum(abs(per_res(1:dall)))/dall
    write(6 ,'(a,2f12.3)') '  Gcs:  ||(Gm-d)||^2  and ||W(Gm-d)||^2: ', res2Nm, resW2Nm
    write(66,'(a,2f12.3)') '  Gcs:  ||(Gm-d)||^2  and ||W(Gm-d)||^2: ', res2Nm, resW2Nm
    ! write(6,*) ' abs mean Res(aa)/Res  (%)',PreRes*100
    ! write(66,*)' abs mean Res(aa)/Res  (%)',PreRes*100
end subroutine


subroutine CalReslNormJoint(maxvp, dall, GVs, GGc, GGs, VsInv, GcInv, GsInv, datweight, Tdata, fwdTvs, fwdTaa, resbst)
    use  lsmrblasInterface, only : dnrm2
    implicit none
    ! resW2Nm: predicated ||W(Gm-d)||2
    ! res2Nm: predicated ||(Gm-d)||2
    ! INPUT
    integer,intent(in):: dall,maxvp
    real,intent(in):: GVs(dall,maxvp), GGc(dall,maxvp), GGs(dall,maxvp)
    real,intent(in):: VsInv(maxvp), GcInv(maxvp), GsInv(maxvp), datweight(dall),Tdata(dall)
    ! OUTPUT
    real:: resbst(dall)
    real:: fwdTvs(dall),fwdTaa(dall)
    ! PARAMETER
    ! real VsInv(maxvp),GcInv(maxvp),GsInv(maxvp)
    real fwdTgc(dall),fwdTgs(dall),fwdt(dall)
    real res(dall),resW(dall),per_res(dall)
    integer jj,kk,ii
    integer i,j,k
    real norm2,norm2W
    real:: res2Nm,resW2Nm
    real:: PreRes
    real:: mean
    !--------------------------------------------------------------------!
    ! Target:
    !  Calculate traveltime resiudal between observe traveltime
    !  and forward traveltime which is generated by inversion result
    !  using G from ref  model.
    !   r=Gm-d
    !--------------------------------------------------------------------!
    ! GcInv=dv(1:maxvp)
    ! GsInv=dv(1+maxvp:2*maxvp)
    !--------------------------------------------------------------------!
    ! forward residual
    fwdTvs=matmul(GVs,VsInv)
    fwdTgc=matmul(GGc,GcInv)
    fwdTgs=matmul(GGs,GsInv)
    fwdTaa=0
    res=0
    resW=0
    per_res=0
    resbst=0
    norm2=0
    norm2W=0
    do i=1,dall
        fwdTaa(i)=fwdTgs(i)+fwdTgc(i)
        resbst(i)=Tdata(i)-fwdTaa(i)-fwdTvs(i)
        resW(i)=resbst(i)*datweight(i) ! here, 1/sigma is for

        ! norm2=norm2+resbst(i)**2
        ! norm2W=norm2W+resW(i)**2
        ! write(6,*)'res=',resbst(i),'weight=',datweight(i),'resW=',resW(i),'norm2=',norm2,'norm2W=',norm2W
    enddo
    resW2Nm=dnrm2(dall,resW,1)
    res2Nm=dnrm2(dall,resbst,1)


    write(6 ,'(a,2f12.3)') '  All:  ||(Gm-d)||^2  and ||W(Gm-d)||^2: ', res2Nm, resW2Nm
    write(66,'(a,2f12.3)') '  All:  ||(Gm-d)||^2  and ||W(Gm-d)||^2: ', res2Nm, resW2Nm
    write(66,'(a)')'  '
    write(6 ,'(a)')'  '
    ! analysis residual
    mean=sum(abs(fwdTaa(1:dall)))/dall
    write(66,'(a,f12.4,a)')'  ABS Mean T(AA): ',mean,'s'
    write(6 ,'(a,f12.4,a)')'  ABS Mean T(AA): ',mean,'s'
    mean=sum(abs(fwdTvs(1:dall)))/dall
    write(66,'(a,f12.4,a)')'  ABS Mean T(dVs):',mean,'s'
    write(6 ,'(a,f12.4,a)')'  ABS Mean T(dVs):',mean,'s'
end subroutine


! Model Norm !

subroutine Calmodel2Norm(nar1,nar,maxvp,count3,rw,col,dv,weight)
    use  lsmrblasInterface, only : dnrm2
    implicit none
    ! works for both Vs and Gcs inversion.
    integer,intent(in) ::nar1,nar,maxvp,count3
    real rw(*)
    integer col(*)
    real,intent(in):: dv(*)
    real,intent(in):: weight
    ! OUTPUT
    real Mnorm2,MwNorm2
    ! PARAMETER
    integer jj,kk,ii
    integer i,j,k
    real Lm(count3),LmWeight(count3)
    integer numC
    ! LOGICAL:: isTikh
    !--------------------------------------------------------------------!
    ! count3: row number of L
    ! Formula: Lm=[a12(m1-m2)
    !                              a13(m1-m3)]
    ! nar1: number of  nonzero value of G
    ! nar: number of nonzero value of Gbar-E
    ! dv: inversion result
    ! rw: value of Gbar
    ! iw2: row number of Gbar, not contain column
    ! col: column number of Gbar
    ! Mnorm2: 2 norm of Lm
    ! numC: count index of  nonzero L
    ! weight: in L, weight has added to Gbar, so divide weight
    ! Regularization Type: (T) 1st order Tikhonov ;(F) Gaussian
    Lm=0
    LmWeight=0

    do i=1,count3
      numC=nar1+i
      Lm(i)=rw(numC)*dv(col(numC))/weight
      LmWeight(i)=rw(numC)*dv(col(numC))
    enddo
    Mnorm2=dnrm2(count3,Lm,1)
    MwNorm2=dnrm2(count3,LmWeight,1)

    ! open(90,file='Lnorm.txt')
    ! do i=1,count3
    ! numC=nar1+i*2-1
    ! write(90,*)numC,col(numC),col(numC+1),dv(col(numC)),dv(col(numC+1)),rw(numC),rw(numC+1)
    ! enddo
    ! close(90)
    !
    ! open(90,file='dv_output.txt')
    ! do i=1,maxvp
    ! write(90,*)dv(i),dv(i+maxvp),dv(i+maxvp*2)
    ! enddo
    ! close(90)

    write(6 ,'(a,2f12.3)') '  dVs:  ||Lm||^2      and ||wLm||^2    : ', Mnorm2, MwNorm2
    write(66,'(a,2f12.3)') '  dVs:  ||Lm||^2      and ||wLm||^2    : ', Mnorm2, MwNorm2
end subroutine

subroutine Calmodel2NormJoint(nar1,nar,maxvp, NreVs, Nre, rw,col,dv, lameGcs, lameVs)
    use  lsmrblasInterface, only : dnrm2
    implicit none
    ! INPUT
    ! NreVs: number of Vs regularization in G
    ! Nre  : number of regularization in G
    integer,intent(in) ::nar1,nar,maxvp, Nre, NreVs
    real rw(*)
    integer col(*)
    real,intent(in):: dv(*)
    real,intent(in):: lameGcs,lameVs
    ! OUTPUT
    real Mnorm2,MwNorm2
    real VsNorm2,VswNorm2
    real GcsNorm2,GcswNorm2
    ! PARAMETER
    integer jj,kk,ii
    integer i,j,k
    real Lm(Nre), LmWeight(Nre)
    integer numC
    ! LOGICAL:: isTikh
    !--------------------------------------------------------------------!
    ! count3: row number of L
    ! Formula: Lm=[a12(m1-m2)
    !                              a13(m1-m3)]
    ! nar1: number of  nonzero value of G
    ! nar: number of nonzero value of Gbar-E
    ! dv: inversion result
    ! rw: value of Gbar
    ! iw2: row number of Gbar, not contain column
    ! col: column number of Gbar
    ! Mnorm2: 2 norm of Lm
    ! numC: count index of  nonzero L
    ! weight: in L, weight has added to Gbar, so divide weight
    ! Regularization Type: (T) 1st order Tikhonov ;(F) Gaussian
    Lm=0
    LmWeight=0

    do i=1, NreVs
      numC=nar1+i
      Lm(i)=rw(numC)*dv(col(numC))/lameVs
      LmWeight(i)=rw(numC)*dv(col(numC))
    enddo
    VsNorm2=dnrm2(NreVs, Lm(1:NreVs), 1)
    VswNorm2=dnrm2(NreVs, LmWeight(1:NreVs), 1)

    do i=NreVs+1, Nre
      numC=nar1+i
      Lm(i)=rw(numC)*dv(col(numC))/lameGcs
      LmWeight(i)=rw(numC)*dv(col(numC))
    enddo

    GcsNorm2 = dnrm2(Nre-NreVs, Lm(NreVs+1:Nre), 1)
    GcswNorm2= dnrm2(Nre-NreVs, LmWeight(NreVs+1:Nre), 1)

    Mnorm2=dnrm2(Nre, Lm, 1)
    MwNorm2=dnrm2(Nre, LmWeight, 1)


    write(6 ,'(a,2f12.3)') '  dVs:  ||Lm||^2   and   ||wLm||^2     : ', VsNorm2, VswNorm2
    write(66,'(a,2f12.3)') '  dVs:  ||Lm||^2   and   ||wLm||^2     : ', VsNorm2, VswNorm2
    write(6 ,'(a,2f12.3)') '  Gcs:  ||Lm||^2   and   ||wLm||^2     : ', GcsNorm2, GcswNorm2
    write(66,'(a,2f12.3)') '  Gcs:  ||Lm||^2   and   ||wLm||^2     : ', GcsNorm2, GcswNorm2
    write(6 ,'(a,2f12.3)') '  All:  ||Lm||^2   and   ||wLm||^2     : ', Mnorm2, MwNorm2
    write(66,'(a,2f12.3)') '  All:  ||Lm||^2   and   ||wLm||^2     : ', Mnorm2, MwNorm2
    write(66,'(a)')'  '
    write(6 ,'(a)')'  '
end subroutine
