text:
!##########################################################################
        subroutine tecplot_T(ki)
!##########################################################################
        use multidata
        implicit none
        integer :: sn,i,j,k,toti,totj,totk
        integer :: is,ie,js,je,ks,ke
        integer :: ki,ib
        character*8 :: chb
        character*25 :: gf

        do ib=1,nbp

        write(chb,'(i4)') dom_id(ib)
        sn=len(trim(adjustl(chb)))
        chb=repeat('0',(4-sn))//trim(adjustl(chb))
        gf='tecout_T_'//trim(adjustl(chb))//'.plt'
        open (unit=88, file=gf)

        is=dom(ib)%isp; ie=dom(ib)%iep
        js=dom(ib)%jsp; je=dom(ib)%jep
        ks=dom(ib)%ksp; ke=dom(ib)%kep
        toti=(ie+1)-(is-1)+1
        totj=(je+1)-(js-1)+1
        totk=(ke+1)-(ks-1)+1

        write (88,*) 'title = Temperature field'
        write (88,*)'variables=x,y,z,T,Tm,Ttm'
        write (88,*)'zone ', ' i=',toti,', ',
     &  ' j=',totj,', k= ',totk,' f=point'
        do k=ks-1,ke+1
           do j=js-1,je+1
              do i=is-1,ie+1
                 write (88,88) dom(ib)%xc(i),dom(ib)%yc(j),dom(ib)%zc(k)
     & ,dom(ib)%T(i,j,k),dom(ib)%Tm(i,j,k),dom(ib)%Ttm(i,j,k)
              end do
           end do
        end do
        close (88)

        end do
88      format (15e25.8)

        end
!##########################################################################
        subroutine tecplot_S(ki)
!##########################################################################
        use multidata
        implicit none
        integer :: sn,i,j,k,toti,totj,totk
        integer :: is,ie,js,je,ks,ke
        integer :: ki,ib
        character*8 :: chb
        character*25 :: gf

        do ib=1,nbp

        write(chb,'(i4)') dom_id(ib)
        sn=len(trim(adjustl(chb)))
        chb=repeat('0',(4-sn))//trim(adjustl(chb))
        gf='tecout_S_'//trim(adjustl(chb))//'.plt'
        open (unit=88, file=gf)

        is=dom(ib)%isp; ie=dom(ib)%iep
        js=dom(ib)%jsp; je=dom(ib)%jep
        ks=dom(ib)%ksp; ke=dom(ib)%kep
        toti=(ie+1)-(is-1)+1
        totj=(je+1)-(js-1)+1
        totk=(ke+1)-(ks-1)+1

        write (88,*) 'title =  Scalar transport field'
        write (88,*)'variables=x,y,z,S,Sm,Stm,SUM,SVM,SWM'
        write (88,*)'zone  i=',toti,',  j=',totj,', k= ',totk,' f=point'

        do k=ks-1,ke+1
           do j=js-1,je+1
              do i=is-1,ie+1
                 write (88,88) dom(ib)%xc(i),dom(ib)%yc(j),dom(ib)%zc(k)
     & ,dom(ib)%S(i,j,k),dom(ib)%Sm(i,j,k),dom(ib)%Stm(i,j,k)
     & ,dom(ib)%SUtm(i,j,k),dom(ib)%SVtm(i,j,k),dom(ib)%SWtm(i,j,k)
              end do
           end do
        end do
        close (88)

        end do
88      format (15e25.8)
        end
!##########################################################################
        subroutine tecplot_p(ki)
!##########################################################################
        use multidata
        implicit none
        integer :: sn,i,j,k,toti,totj,totk
        integer :: is,ie,js,je,ks,ke
        integer :: ki,ib
        character*8 :: chb
        character*25 :: gf

        do ib=1,nbp

        write(chb,'(i4)') dom_id(ib)
        sn=len(trim(adjustl(chb)))
        chb=repeat('0',(4-sn))//trim(adjustl(chb))
        gf='tecout_p'//trim(adjustl(chb))//'.plt'
        open (unit=88, file=gf)

        is=dom(ib)%isp; ie=dom(ib)%iep
        js=dom(ib)%jsp; je=dom(ib)%jep
        ks=dom(ib)%ksp; ke=dom(ib)%kep
        toti=(ie+1)-(is-1)+1
        totj=(je+1)-(js-1)+1
        totk=(ke+1)-(ks-1)+1

        write (88,*) 'title = Pressure field'
        write (88,*)'variables=x,y,z,P,PM,ppm,vis,visM,ksgs,eps,epsm'
        write (88,*)'zone ', ' i=',toti,', ',
     &  ' j=',totj,', k= ',totk,' f=point'

        do k=ks-1,ke+1
           do j=js-1,je+1
              do i=is-1,ie+1
                 write (88,88) dom(ib)%xc(i),dom(ib)%yc(j),
     & dom(ib)%zc(k),dom(ib)%p(i,j,k),dom(ib)%pm(i,j,k),
     & dom(ib)%ppm(i,j,k),dom(ib)%vis(i,j,k),dom(ib)%vism(i,j,k),
     & dom(ib)%ksgs(i,j,k),dom(ib)%eps(i,j,k),dom(ib)%epsm(i,j,k)
              end do
           end do
        end do
        close (88)

        end do
88      format (15e25.8)

        end
!##########################################################################
        subroutine tec_turb(ki)
!##########################################################################
        use multidata
        use vars
        implicit none
        integer :: sn,i,j,k,toti,totj,totk
        integer :: is,ie,js,je,ks,ke
        integer :: ki,ib
        double precision u_cn,v_cn,w_cn,um_cn,vm_cn,wm_cn
        double precision uum_cn,vvm_cn,wwm_cn,uvml,uwml,vwml
        integer :: itoti,itotj,itotk,ntoti,ntotj,ntotk
        character*8 :: chb
        character*25 :: gf

        do ib=1,nbp

        write(chb,'(i4)') dom_id(ib)
        sn=len(trim(adjustl(chb)))
        chb=repeat('0',(4-sn))//trim(adjustl(chb))
        gf='tecturb'//trim(adjustl(chb))//'.plt'
        open (unit=88, file=gf)

        is=pl+1; ie=dom(ib)%ttc_i-pl
        js=pl+1; je=dom(ib)%ttc_j-pl
        ks=pl+1; ke=dom(ib)%ttc_k-pl
        toti=ie-(is-1)+1
        totj=je-(js-1)+1
        totk=ke-(ks-1)+1

!Plot every X points in every direction: change itotX respectively
	  itoti=1 	; ntoti=toti/itoti+1
	  itotj=1 	; ntotj=totj/itotj+1
	  itotk=1 	; ntotk=totk/itotk+1

	IF(.NOT. L_LSM) THEN
        write (88,*) 'title = turb'
        write (88,*)'variables=x,y,z,U,V,W,UM,VM,WM'
!     &, ,'uuM,vvM,wwM,uvM,uwM,vwM'
        write (88,*)'zone ', ' i=',ntoti,', ',
     &  ' j=',ntotj,', k= ',ntotk,' f=point'
	ELSE !WHEN LSM IS ALSO ACTIVATED
        write (88,*) 'title =  turb'
        write (88,*)'variables=x,y,z,U,V,W,phi,p'
        write (88,*)'zone ', ' i=',ntoti,', ',
     &  ' j=',ntotj,', k= ',ntotk,' f=point'
	ENDIF

        do k=ks-1,ke+1,itotk
           do j=js-1,je+1,itotj
              do i=is-1,ie+1,itoti

                 u_cn  =0.25*(dom(ib)%u(i,j,k)+
     &dom(ib)%u(i,j+1,k)+dom(ib)%u(i,j,k+1)+
     &dom(ib)%u(i,j+1,k+1))
                 um_cn  =0.25*(dom(ib)%um(i,j,k)+
     &dom(ib)%um(i,j+1,k)+dom(ib)%um(i,j,k+1)+
     &dom(ib)%um(i,j+1,k+1))
                 uum_cn  =0.25*(dom(ib)%uum(i,j,k)+
     &dom(ib)%uum(i,j+1,k)+dom(ib)%uum(i,j,k+1)+
     &dom(ib)%uum(i,j+1,k+1))

                 v_cn  =0.25*(dom(ib)%v(i,j,k)+
     &dom(ib)%v(i+1,j,k)+dom(ib)%v(i,j,k+1)+
     &dom(ib)%v(i+1,j,k+1))
                 vm_cn  =0.25*(dom(ib)%vm(i,j,k)+
     &dom(ib)%vm(i+1,j,k)+dom(ib)%vm(i,j,k+1)+
     &dom(ib)%vm(i+1,j,k+1))
                 vvm_cn  =0.25*(dom(ib)%vvm(i,j,k)+
     &dom(ib)%vvm(i+1,j,k)+dom(ib)%vvm(i,j,k+1)+
     &dom(ib)%vvm(i+1,j,k+1))

                 w_cn  =0.25*(dom(ib)%w(i,j,k)+
     &dom(ib)%w(i+1,j,k)+dom(ib)%w(i,j+1,k)+
     &dom(ib)%w(i+1,j+1,k))
                 wm_cn  =0.25*(dom(ib)%wm(i,j,k)+
     &dom(ib)%wm(i+1,j,k)+dom(ib)%wm(i,j+1,k)+
     &dom(ib)%wm(i+1,j+1,k))
                 wwm_cn  =0.25*(dom(ib)%wwm(i,j,k)+
     &dom(ib)%wwm(i+1,j,k)+dom(ib)%wwm(i,j+1,k)+
     &dom(ib)%wwm(i+1,j+1,k))

                 uvml  =0.125*(dom(ib)%uvm(i,j,k)+
     &dom(ib)%uvm(i+1,j,k)    +dom(ib)%uvm(i,j+1,k)+
     &dom(ib)%uvm(i+1,j+1,k)  +dom(ib)%uvm(i,j,k+1)+
     &dom(ib)%uvm(i+1,j,k+1)  +dom(ib)%uvm(i,j+1,k+1)+
     &dom(ib)%uvm(i+1,j+1,k+1))
                 uwml  =0.125*(dom(ib)%uwm(i,j,k)+
     &dom(ib)%uwm(i+1,j,k)    +dom(ib)%uwm(i,j+1,k)+
     &dom(ib)%uwm(i+1,j+1,k)  +dom(ib)%uwm(i,j,k+1)+
     &dom(ib)%uwm(i+1,j,k+1)  +dom(ib)%uwm(i,j+1,k+1)+
     &dom(ib)%uwm(i+1,j+1,k+1))
                 vwml  =0.125*(dom(ib)%vwm(i,j,k)+
     &dom(ib)%vwm(i+1,j,k)    +dom(ib)%vwm(i,j+1,k)+
     &dom(ib)%vwm(i+1,j+1,k)  +dom(ib)%vwm(i,j,k+1)+
     &dom(ib)%vwm(i+1,j,k+1)  +dom(ib)%vwm(i,j+1,k+1)+
     &dom(ib)%vwm(i+1,j+1,k+1))

	IF(.NOT. L_LSM) THEN
          write (88,88) dom(ib)%x(i),dom(ib)%y(j),
     &  dom(ib)%z(k),u_cn,v_cn,w_cn,um_cn,vm_cn,wm_cn
!     &  ,uum_cn,vvm_cn,wwm_cn,uvml,uwml,vwml
	else
          write (88,88) dom(ib)%x(i),dom(ib)%y(j), dom(ib)%z(k)
     & ,u_cn,v_cn,w_cn,dom(ib)%phi(i,j,k),dom(ib)%p(i,j,k)
	ENDIF
              end do
           end do
        end do

        close (88)

        end do
88      format (20e25.8)
        end
!##########################################################################
        subroutine tecplot_phi(ki)
!##########################################################################
        use vars
        use multidata
        implicit none

        integer :: sn,sn1,i,j,k,toti,totj,totk
        integer :: is,ie,js,je,ks,ke
        integer :: ki,ib
        character*8 :: chb,chb1
        character*27 :: gf

        do ib=1,nbp

	  if (L_anim_phi) then
          write(chb,'(i4)') dom_id(ib)
          write(chb1,'(i6)') ki
          sn=len(trim(adjustl(chb)))
          sn1=len(trim(adjustl(chb1)))
          chb=repeat('0',(4-sn))//trim(adjustl(chb))
          chb1=repeat('0',(6-sn1))//trim(adjustl(chb1))
          gf='tecout_phi_'//trim(adjustl(chb))//'_'//
     & trim(adjustl(chb1))//'.plt'
	  else
          write(chb,'(i4)') dom_id(ib)
          sn=len(trim(adjustl(chb)))
          chb=repeat('0',(4-sn))//trim(adjustl(chb))
          gf='tecout_phi_'//trim(adjustl(chb))//'.plt'
	  endif

	  open (unit=88, file=gf)

        is=dom(ib)%isp; ie=dom(ib)%iep
        js=dom(ib)%jsp; je=dom(ib)%jep
        ks=dom(ib)%ksp; ke=dom(ib)%kep
        toti=ie-(is-1)+1
        totj=je-(js-1)+1
        totk=ke-(ks-1)+1

        write (88,*)'title = Free-surface from LSM'
        write (88,*)'variables=x,y,z,phi,phim,dens,mu'
	  if (L_anim_phi) then
          write (88,*)'zone ','STRANDID=', 1, 'SOLUTIONTIME=', ctime,
     &    ' i=',toti,', ',' j=',totj,', k= ',totk,
     &    'zonetype=', 'ordered',', DATAPACKING=point'
	  else
          write(88,*)'zone  i=',toti,', j=',totj,', k= ',totk,' f=point'
	  endif

        do k=ks-1,ke
           do j=js-1,je
              do i=is-1,ie
                write (88,88) dom(ib)%xc(i),dom(ib)%yc(j),dom(ib)%zc(k)
     & ,dom(ib)%phi(i,j,k),dom(ib)%phim(i,j,k),dom(ib)%dens(i,j,k),
     & dom(ib)%mu(i,j,k)
              end do
           end do
        end do
        close (88)

        end do

88      format (8e20.8)

        end
!##########################################################################
        subroutine tecbin(ki)
!##########################################################################
        use multidata
        use vars
        implicit none
        integer :: sn,i,j,k,toti,totj,totk
        integer :: ki,ib,inind,jnind,knind
        character*8 :: chb
        character*25 :: gf

        do ib=1,nbp

        write(chb,'(i4)') dom_id(ib)
        sn=len(trim(adjustl(chb)))
        chb=repeat('0',(4-sn))//trim(adjustl(chb))
        gf='tecbin'//trim(adjustl(chb))//'.bin'
        open (unit=88, file=gf, form='unformatted')

        toti=dom(ib)%ttc_i
        totj=dom(ib)%ttc_j
        totk=dom(ib)%ttc_k

        write (88) toti,totj,totk
        write (88) pl
!====================================================================
        inind=0; jnind=0; knind=0
        if (dom(ib)%inext.lt.0 .and. dom(ib)%bc_east.ne.5) inind=-1
        if (dom(ib)%jnext.lt.0 .and. dom(ib)%bc_north.ne.5) jnind=-1
        if (dom(ib)%knext.lt.0 .and. dom(ib)%bc_top.ne.5) knind=-1
        write (88) inind,jnind,knind
!====================================================================

        do k=1,totk
           do j=1,totj
              do i=1,toti

        write (88) dom(ib)%x(i),dom(ib)%y(j),dom(ib)%z(k),
     & dom(ib)%p(i,j,k),dom(ib)%pm(i,j,k),dom(ib)%ppm(i,j,k),
     & dom(ib)%vis(i,j,k),dom(ib)%vism(i,j,k),
     & dom(ib)%u(i,j,k),dom(ib)%um(i,j,k),dom(ib)%uum(i,j,k),
     & dom(ib)%v(i,j,k),dom(ib)%vm(i,j,k),dom(ib)%vvm(i,j,k),
     & dom(ib)%w(i,j,k),dom(ib)%wm(i,j,k),dom(ib)%wwm(i,j,k),
     & dom(ib)%uvm(i,j,k),dom(ib)%uwm(i,j,k),dom(ib)%vwm(i,j,k),
     & dom(ib)%ksgs(i,j,k),dom(ib)%eps(i,j,k),dom(ib)%epsm(i,j,k)
              end do
           end do
        end do
        close (88)

        end do
        end
!##########################################################################
        subroutine tec_instant(ki)
!##########################################################################
        use multidata
        use vars
	use mpi
        implicit none
        integer :: sn,i,j,k,toti,totj,totk
        integer :: is,ie,js,je,ks,ke
        integer :: ki,ib,itoti,itotj,itotk,ntoti,ntotj,ntotk
        character*8 :: chb,chb2
        character*35 :: gf

!	if (nbp.gt.1) RETURN

        do ib=1,nbp

!	  if(rdiv(dom_id(ib)).le.2) RETURN

        write(chb,'(i4)') dom_id(ib)
        sn=len(trim(adjustl(chb)))
        chb=repeat('0',(4-sn))//trim(adjustl(chb))
        write(chb2,'(i6)') ki
        sn=len(trim(adjustl(chb2)))
        chb2=repeat('0',(6-sn))//trim(adjustl(chb2))
        gf='tecinst'//trim(adjustl(chb2))//
     &		'_'//trim(adjustl(chb))//'.plt'
        open (unit=95, file=gf)

        is=pl+1; ie=dom(ib)%ttc_i-pl
        js=pl+1; je=dom(ib)%ttc_j-pl
        ks=pl+1; ke=dom(ib)%ttc_k-pl
        toti=ie-(is-1)+1	!toti=nicell
        totj=je-(js-1)+1
        totk=ke-(ks-1)+1

!Plot every X points in every direction: change itotX respectively
	   itoti=1 	; ntoti=toti/itoti+1
	   itotj=1 	; ntotj=totj/itotj+1
	   itotk=1 	; ntotk=totk/itotk+1

        write(95,*)'variables=x,y,z,U,V,W'
        write(95,89)'zone i=',ntoti,', j=',ntotj,', k=',ntotk,' f=point'

        do k=ks-1,ke+1,itotk
           do j=js-1,je+1,itotj
              do i=is-1,ie+1,itoti
             write (95,88) dom(ib)%x(i),dom(ib)%y(j),dom(ib)%z(k),
     &  dom(ib)%u(i,j,k),dom(ib)%v(i,j,k),dom(ib)%w(i,j,k)
              end do
           end do
        end do

        close (95)
        end do

88      format (10e15.5)
89      format (a7,i4,a4,i4,a4,i4,a8)

        end
!##########################################################################
        subroutine tec_inst_plane(ki)
!##########################################################################
        use multidata
        use vars
	  use mpi
        implicit none
        integer :: sn,i,j,k,toti,totj,totk,is,ie,js,je,ks,ke,ki,ib
        integer :: itoti,itotj,itotk,ntoti,ntotj,ntotk
	  double precision :: xval,yval,zval
	  integer :: iplane
        character*8 :: chb,chb2
        character*35 :: gf

        do ib=1,nbp

        is=pl+1; ie=dom(ib)%ttc_i-pl
        js=pl+1; je=dom(ib)%ttc_j-pl
        ks=pl+1; ke=dom(ib)%ttc_k-pl
        toti=ie-(is-1)+1 ; itoti=1 	; ntoti=toti/itoti+1
        totj=je-(js-1)+1 ; itotj=1 	; ntotj=totj/itotj+1
        totk=ke-(ks-1)+1 ; itotk=1 	; ntotk=totk/itotk+1

!!!! Set the plane orientation and location at xval, yval or zval
	  iplane=2 			!1=x-plane, 2=y-plane, 3=z-plane
	  xval=(xen-xst)/2 ;  yval=(yen-yst)/2  ;  zval=(zen-zst)/2

	 if(iplane.eq.1) then
        do i=is-1,ie
	    if(dom(ib)%x(i).le.xval .and. dom(ib)%x(i+1).gt.xval) then

        write(chb,'(i4)') dom_id(ib)
        sn=len(trim(adjustl(chb)))
        chb=repeat('0',(4-sn))//trim(adjustl(chb))
        write(chb2,'(i6)') ki
        sn=len(trim(adjustl(chb2)))
        chb2=repeat('0',(6-sn))//trim(adjustl(chb2))
        gf='tecplane'//trim(adjustl(chb2))//
     &		'_'//trim(adjustl(chb))//'.plt'
           open (unit=95, file=gf)
        write(95,*)'variables=y,z,u,v,w'
        write(95,89)'zone i= 1, j=',ntotj,', k=',ntotk,' f=point'
           do k=ks-1,ke+1,itotk	!k=pl,nicell+pl
             do j=js-1,je+1,itotj
             write (95,88) dom(ib)%y(j),dom(ib)%z(k),
     &  dom(ib)%u(i,j,k),dom(ib)%v(i,j,k),dom(ib)%w(i,j,k)
             end do
           end do
           close (95)
	    endif ! within the interval
	  enddo
	 endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
	 if(iplane.eq.2) then
        do j=js-1,je
	    if(dom(ib)%y(j).le.yval .and. dom(ib)%y(j+1).gt.yval) then

        write(chb,'(i4)') dom_id(ib)
        sn=len(trim(adjustl(chb)))
        chb=repeat('0',(4-sn))//trim(adjustl(chb))
        write(chb2,'(i6)') ki
        sn=len(trim(adjustl(chb2)))
        chb2=repeat('0',(6-sn))//trim(adjustl(chb2))
        gf='tecplane'//trim(adjustl(chb2))//
     &		'_'//trim(adjustl(chb))//'.plt'
           open (unit=95, file=gf)
        write(95,*)'variables=x,z,u,v,w'
        write(95,89)'zone i=',ntoti,', j=1 , k=',ntotk,' f=point'
           do k=ks-1,ke+1,itotk
             do i=is-1,ie+1,itoti
             write (95,88) dom(ib)%x(i),dom(ib)%z(k),
     &  dom(ib)%u(i,j,k),dom(ib)%v(i,j,k),dom(ib)%w(i,j,k)
!     &  ,dom(ib)%phi(i,j,k)
             end do
           end do
           close (95)
	    endif ! within the interval
	  enddo
	 endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
	 if(iplane.eq.3) then
        do k=ks-1,ke
	    if(dom(ib)%z(k).le.zval .and. dom(ib)%z(k+1).gt.zval) then

        write(chb,'(i4)') dom_id(ib)
        sn=len(trim(adjustl(chb)))
        chb=repeat('0',(4-sn))//trim(adjustl(chb))
        write(chb2,'(i6)') ki
        sn=len(trim(adjustl(chb2)))
        chb2=repeat('0',(6-sn))//trim(adjustl(chb2))
        gf='tecplane'//trim(adjustl(chb2))//
     &		'_'//trim(adjustl(chb))//'.plt'
           open (unit=95, file=gf)
        write(95,*)'variables=x,y,u,v,w'
        write(95,89)'zone i=',ntoti,', j=',ntotj,', k=1 f=point'
           do j=js-1,je+1,itotj
            do i=is-1,ie+1,itoti
             write (95,88) dom(ib)%x(i),dom(ib)%y(j),
     &  dom(ib)%u(i,j,k),dom(ib)%v(i,j,k),dom(ib)%w(i,j,k)
             end do
           end do
           close (95)
	    endif ! within the interval
	  enddo
	 endif

	enddo !ib

88      format (10e15.5)
89      format (a7,i4,a4,i4,a4,i4,a8)

        end

ext: for


text:
      SUBROUTINE logsex( n, buf, lbuf, ibuf, scr)
c
c  ARGUMENTS:
c  buf   - The trace, with SEGY header as TYPE REAL
c  lbuf  - The trace, with SEGY header as TYPE INTEGER*4
c  ibuf  - The trace, with SEGY header as TYPE INTEGER*2
c  scr   - A scratch array.
c
c  COPYRIGHT (C) Seismic Reflection Processors, Solana Beach, CA. 92075
c  ALL RIGHTS RESERVED.
c  Written by Graham M. Kent, SRP, 16 September 1991
c  mod 1 Apr 09 - Use numdat in common for nsamps rather than the header.
c
      DIMENSION buf(1), lbuf(1), ibuf(1), scr(1)
      INTEGER*2 ibuf
      COMMON /readt/ itunit, numhdr, numdat, ihunit, ireeln, jntrcs,
     *               ifmt, nskip, secs, lrenum, isrcf, idtype,
     *               nfskip, jform, itxsi, itxdel, nfktrc, norigtr
      COMMON /sioap/ iasgnd, irelse, in, iout, nextad, lapsiz, ifree,
     *     iuseap
      COMMON /segyptr/ llsegptr, lrseqptr, lshotptr, lshtrptr, lrpnptr,
     *                 lrptrptr, itridptr, ldisptr,  lwbdptr,  lsxcoptr,
     *                 lrxcoptr, idelmptr, istmptr,  iendmptr, isampptr,
     *                 isiptr,   iyrptr,   idayptr,  ihrptr,   iminptr,
     *                 isecptr,  igmtptr,  ldelsptr, lsmusptr,lemusptr,
     *                 lsisptr,  lwbtsptr, lgatptr,  lssmsptr, lesmsptr,
     *                 lsbptr,   ifoldptr, icvleptr, lespnptr
      COMMON /apmem/ a(32767)
      DIMENSION params1(8), lparams1(8), xdata(10000), ydata(10000),
     *yderv(10000)
      EQUIVALENCE (params1(1),lparams1(1))
      INTEGER logtype
      COMMON /logcom/ ilogunit1, ilogunit2, nloglists, nlogwrds
      SAVE
      DATA mlists1/0/, mlists2/0/,xdata/10000*0./, ydata/10000*0./,
     *yderv/10000*0./
c
c
      IF( ibuf(itridptr) .EQ. 2 ) RETURN                                ! forget dead traces
c*GMK Have we read a list for logst1, if not get buf variables once - write over them
      IF (mlists1.eq.0) THEN
c           nsamps1 = ibuf(isampptr)
           nsamps1 = numdat
           delay1 = buf(ldelsptr)
           si1 = buf(lsisptr)                                                 ! sample interval in seconds
           wbt1 = buf(lwbtsptr)                                               ! water bottom time in seconds
           goto 10
      ELSE
c*GMK Have we read a list for logst2, if not get buf variables once - write over them
         IF (mlists2.eq.0) THEN
c           nsamps2 = ibuf(isampptr)
           nsamps2 = numdat
           delay2 = buf(ldelsptr)
           si2 = buf(lsisptr)                                                 ! sample interval in seconds
           wbt2 = buf(lwbtsptr)                                               ! water bottom time in seconds
           goto 10
         ELSE
           goto 11
         ENDIF
      ENDIF

c*GMK Checks logstretch and reads appropriate disk file
10    IF( n .EQ. 1 ) THEN
          CALL podisc( ilogunit1, 1, 0 )
          CALL rddisc( ilogunit1, params1, nlogwrds, istat )
          mlists1 = mlists1 + 1
          logtype1  = lparams1(1)
          rtsamp11  = params1(2)
          rtsamp21  = params1(3)
          rtcut1    = params1(4)
          rsltime1  = params1(5)
          reltime1  = params1(6)
          rloghz1   = params1(7)
          lprint1   = lparams1(8)
c*GMK Lprint statement n=1
         IF( lprint1 .EQ. 1 ) THEN
          PRINT *,'  type=',logtype1,' tsamp1=',rtsamp11,' tsamp2=',
     &rtsamp22,' tcut=',rtcut1,' sltime=',rsltime1,' eltime=',reltime1,
     &' loghz=', rloghz1
         ENDIF

      ELSE
          CALL podisc( ilogunit2, 1, 0 )
          CALL rddisc( ilogunit2, params1, nlogwrds, istat )
          mlists2 = mlists2 + 1
          logtype2  = lparams1(1)
          rtsamp12  = params1(2)
          rtsamp22  = params1(3)
          rtcut2    = params1(4)
          rsltime2  = params1(5)
          reltime2  = params1(6)
          rloghz2   = params1(7)
          lprint2   = lparams1(8)
c*GMK Lprint statement n=2
         IF( lprint2 .EQ. 1 ) THEN
          PRINT *,'  type=',logtype2,' tsamp1=',rtsamp12,' tsamp2=',
     &rtsamp22,' tcut=',rtcut2,' sltime=',rsltime2,' eltime=',reltime2,
     &' loghz=', rloghz2
         ENDIF

      ENDIF

c*GMK Determine what type of log stetch is needed
c*    And get appropriate variables which are generalized

11    if( n .eq. 1 ) then
          logtype =  logtype1
          rtsamp1 =  rtsamp11
          rtsamp2 =  rtsamp21
          rtcut   =  rtcut1
          rsltime =  rsltime1
          reltime =  reltime1
          rloghz  =  rloghz1
          lprint  =  lprint1
        goto 500
      else
          logtype =  logtype2
          rtsamp1 =  rtsamp12
          rtsamp2 =  rtsamp22
          rtcut   =  rtcut2
          rsltime =  rsltime2
          reltime =  reltime2
          rloghz  =  rloghz2
          lprint  =  lprint2
        goto 501
      endif

c*GMK Going to stretch the data prior to DMO
c*    *******************PROCESS LOGST1***********************************

500   continue

c*GMK Going to log stretch the data and ensure that TCUT is
c*    used so as to not take the LN of zero. Data are then
c*    place from array buf into ydata for spline routine and
c*    new stretched times are calculated for array xdata

      if (rtcut.lt.delay1) then
         do 100 i = 1, nsamps1
          t = delay1 + (i-1)*(si1/1.0)
          xdata(i) = log(t/rtcut)
c*GMK     Is data in AP
          if (in.eq.0) then
           ydata(i) = buf(i+numhdr)
          else
           ydata(i) = a(i)
          endif
100   continue
      nspline = nsamps1
      slog = xdata(1)
      elog = xdata(nspline)
      else

c*GMK IF tcut is greater, how much data do I need to skip before tcut - jump

          jump = nint( ( (rtcut-delay1)*1.0 )/si1) +1
         do 110 j = jump, nsamps1
          t = delay1 + (j-1)*(si1/1.0)
          xdata(j-jump+1) = log(t/rtcut)
c*GMK     Is data in AP
          if (in.eq.0) then
            ydata(j-jump+1) = buf(j+numhdr)
          else
            ydata(j-jump+1) = a(j)
          endif

110   continue
      nspline = (nsamps1-jump) + 1
      slog = xdata(1)
      elog = xdata(nspline)
      endif

c*GMK Big requires the natural spline to be taken (i= -1 and n+1 are 0)
      big = 1.0E+31

c*GMK Subroutine spline calculates the necessary 2nd derevitives
      call spline(xdata,ydata,nspline,big,big,yderv)

c*GMK Once the partials are calculated, we need to spline the data on
c*    to the new trace of sample interval tsamp1. New Trace will start a time
c*    t = 0 (going to FK anyways) and final time given by stretch function.

c*GMK Zero scratch array SCR to place final interpolated trace values
      do 120 k = 1, 10000
         scr(k) = 0.0
120   continue

c*GMK Find first sample greater than first value for log stretch, and do the same
c*    for largest value but firsat sample less than greatest log time value.

      islog = nint( (slog*1.0)/rtsamp1 ) + 1
      ielog = nint( (elog*1.0)/rtsamp1 )

c*GMK Calculate new value of y at every new evenly spaced x starting at time
c*    corresponding to sample islog

      do 200 n = islog, ielog
       x = (rtsamp1/1.0)*(n-1)
       call splint(xdata,ydata,yderv,nspline,x,y)
       scr(n) = y
200   continue

      do 210 i = 1, ielog
c*GMK     Is data in AP
       if (in.eq.0) then
        buf(i+numhdr) = scr(i)
       else
        a(i) = scr(i)
       endif
210   continue

c*GMK Now change header values (Ouch that's dangerous) so the next PROCESS
c*    has a clue to the tweeked trace.

      ibuf(isampptr)  = ielog
      numdat = ielog
      buf(ldelsptr)   = 0.0
      ibuf(idelmptr)  = 0
      buf(lsisptr)    = rtsamp1
      ibuf(isiptr)    = rtsamp1*1.0E+06
c*GMK finished leave process logstr1
      goto 999

c*GMK Enter at 501 if data needs to be compressed subsequent to DMO
501   continue

c*GMK Going to log compress the data.  Data are then
c*    place from array buf into ydata for spline routine and
c*    new compressed times are calculated for array xdata

c*    *******************PROCESS LOGSTR2************************************

c*    This case should not happen because trace should have no delay after DMO
c*    but you never know!

      if (rtcut.lt.delay2) then
         do 1000 i = 1, nsamps2
          tau = delay2 + (i-1)*(si2/1.0)
          xdata(i) = rtcut*exp(tau)
c*GMK     Is data in AP
          if (in.eq.0) then
           ydata(i) = buf(i+numhdr)
          else
           ydata(i) = a(i)
          endif
1000  continue
      nspline = nsamps2
      slog = xdata(1)
      elog = xdata(nspline)
      else

c*GMK rtcut greater than delay
c*    This should always be true for this process

         do 1100 j = 1, nsamps2
C*GMK    t is really tau here because of compression
          tau = delay2 + (j-1)*(si2/1.0)
          xdata(j) = rtcut*exp(tau)
c*GMK     Is data in AP
          if (in.eq.0) then
           ydata(j) = buf(j+numhdr)
          else
           ydata(j) = a(j)
          endif
1100  continue
      nspline = nsamps2
      slog = xdata(1)
      elog = xdata(nspline)
      endif

c*GMK Big requires the natural spline to be taken (i= -1 and n+1 are 0)
      big = 1.0E+31
c*GMK Subroutine spline calculates the necessary 2nd derevitives
      call spline(xdata,ydata,nspline,big,big,yderv)

c*GMK Once the partials are calculated, we need to spline the data on
c*    to the new trace of sample interval tsamp1. New Trace will start a time
c*    t = 0 (going to FK anyways) and final time given by stretch function.

c*GMK Zero scratch array SCR to place final interpolated trace values
      do 1200 k = 1, 10000
         scr(k) = 0.0
1200  continue

c*GMK Find first sample greater than first value for log stretch, and do the same
c*    for largest value but firsat sample less than greatest log time value.

      islog = nint( (slog)/rtsamp2 ) + 1
      ielog = nint( (elog)/rtsamp2 )

c*GMK Calculate new value of y at every new evenly spaced x starting at time
c*    corresponding to sample islog

      do 2000 n = islog, ielog
       x = (rtsamp2/1.0)*(n-1)
       call splint(xdata,ydata,yderv,nspline,x,y)
       scr(n) = y
2000  continue

c*GMK Account for new tsamp rate, start and end times
c*    Going to have to tweeked the vales in common blks also to be consistent

      istart =  nint ( (rsltime-delay2)*(1.0/rtsamp2) ) + 1
      ifinish = nint ( (reltime-rsltime)*(1.0/rtsamp2) )
     &+istart
      do 2100 i = istart, ifinish
c*GMK     Is data in AP
          if (in.eq.0) then
           buf(i-istart+1+numhdr) = scr(i)
          else
           a(i-istart+1) = scr(i)
          endif
2100  continue

c*GMK Now change header values (Ouch that's dangerous) so the next PROCESS
c*    has a clue to the tweeked trace.

      numdat = (ifinish-istart) + 1
      ibuf(isampptr)  = numdat
      buf(ldelsptr)   = rsltime
      ibuf(idelmptr)  = rsltime*1000
      buf(lsisptr)    = rtsamp2
      ibuf(isiptr)    = rtsamp1*1.0E+06
c*GMK finished leave process logstr2
      goto 999

999   continue
      return
      end

      subroutine spline(x,y,n,yp1,ypn,y2)
      parameter (nmax=10000)
      dimension x(n),y(n),y2(n),u(nmax)
      if (yp1.gt..99e30) then
        y2(1)=0.
        u(1)=0.
      else
        y2(1)=-0.5
        u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
      endif
      do 11 i=2,n-1
        sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
        p=sig*y2(i-1)+2.
        y2(i)=(sig-1.)/p
        u(i)=(6.*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1))
     *      /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
11    continue
      if (ypn.gt..99e30) then
        qn=0.
        un=0.
      else
        qn=0.5
        un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
      endif
      y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
      do 12 k=n-1,1,-1
        y2(k)=y2(k)*y2(k+1)+u(k)
12    continue
      return
      end

      subroutine splint(xa,ya,y2a,n,x,y)
      dimension xa(n),ya(n),y2a(n)
      klo=1
      khi=n
1     if (khi-klo.gt.1) then
        k=(khi+klo)/2
        if(xa(k).gt.x)then
          khi=k
        else
          klo=k
        endif
      goto 1
      endif
      h=xa(khi)-xa(klo)
      if (h.eq.0.) THEN
         PRINT *,'****  ERROR  ****  bad xa input.'
         CALL EXIT
      ENDIF
      a=(xa(khi)-x)/h
      b=(x-xa(klo))/h
      y=a*ya(klo)+b*ya(khi)+
     *      ((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**2)/6.
      return
      end

ext: f

max_stars_repo_name: apthorpe/sioseis

max_issues_repo_name: henkart/sioseis

max_forks_repo_name: henkart/sioseis


text:
module Ref {

  @ Input port for receiving payload power on/off signal
  port receive_power_command_input(
                                    power: bool @< boolean val, true -> ON, false -> OFF
                                  ) -> cmdRecv

  @ Output port for sending success/failure in carrying out power operation
  port send_power_command_status_output(
                                        status: bool @< boolean val, true -> SUCCESS, false -> FAILURE
                                        )

  @ Output port for requesting current payload power status
  port request_power_status_output(
                                    req_status: bool @< boolean val, true -> get status
                                  )

  @ Input port for receiving the on/off status
  port receive_power_status_input(
                                  status: bool @< boolean val, true -> ON, false -> OFF
                                  ) -> cmdRecv

  @ Output port for sending a power on/off signal to EPS hardware
  port send_power_command_output(
                                  power: bool @< boolean val, true -> ON, false -> OFF
                                )

}

ext: fpp

max_stars_repo_name: MelbourneSpaceProgram/fprime

max_issues_repo_name: MelbourneSpaceProgram/fprime

max_forks_repo_name: MelbourneSpaceProgram/fprime


text:
! $Id$
! -----------------------------------------------------------
! Module McWeeny
! -----------------------------------------------------------
! Code area 4: density matrix
! -----------------------------------------------------------

!!****h* Conquest/McWeeny
!!  NAME
!!   McWeeny
!!  PURPOSE
!!   To handle the McWeeny/Manolopoulos initialisation of density matrix
!!  AUTHOR
!!   D.R.Bowler
!!  CREATION DATE
!!   21/10/98
!!  MODIFICATION HISTORY
!!   11/05/01 DRB
!!    ROBODoc headers and better output
!!   08/06/2001 dave
!!    Added RCS Id and Log tags; used GenComms throughout
!!    Changed my_barrier to come from GenComms
!!   11:09, 04/02/2003 drb
!!    Major change: stability fix in InitMcW to ensure mu lies between Hamiltonian limits
!!   10:09, 13/02/2006 drb
!!    Removed all explicit references to data_ variables and rewrote in terms of new
!!    matrix routines
!!   2006/11/21 16:53 dave
!!    Changed iprint statements
!!   2008/02/01 17:48 dave
!!    Changes for output to file not stdout
!!   2008/07/31 ast
!!    Added timers
!!   2012/03/01 L.Tong
!!    Added interface for McW_matrix_multiply
!!   2014/09/15 18:30 lat
!!    fixed call start/stop_timer to timer_module (not timer_stdlocks_module !)
!!   2019/10/24 10:30 dave
!!    Added inode and ionode as use-associated from GenComms in module header
!!    for efficiency and best practice; changed function calls
!!***
module McWeeny

  use global_module,          only: io_lun, area_DM
  use timer_module,           only: start_timer, stop_timer
  use timer_stdclocks_module, only: tmr_std_matrices
  use GenComms,               ONLY: inode, ionode

  implicit none

  !!****f* McWeeny/McW_matrix_multiply
  !! PURPOSE
  !!   Interface for McW_matrix_multiply_nondsw and
  !!   McW_matrix_multiply_ds
  !! AUTHOR
  !!   L.Tong
  !! CREATION DATE
  !!   2012/03/01
  !! MODIFICATION HISTORY
  !! SOURCE
  !!
  interface McW_matrix_multiply
     module procedure McW_matrix_multiply_nonds
     module procedure McW_matrix_multiply_ds
  end interface McW_matrix_multiply
  !!*****

contains

  !!****f* McWeeny/McWMin *
  !!
  !!  NAME
  !!   McWMin
  !!  USAGE
  !!
  !!  PURPOSE
  !!   Performs iterative minimisation using McWeeny algorithm
  !!  INPUTS
  !!
  !!
  !!  USES
  !!   datatypes, numbers, logicals, maxima_module, matrix_data
  !!   mult_module, multiply_module, GenBlas
  !!  AUTHOR
  !!   D.R.Bowler
  !!  CREATION DATE
  !!   21/10/98
  !!  MODIFICATION HISTORY
  !!   11/05/01 DRB
  !!    ROBODoc headers and better output also nsf from maxima_module
  !!   08/06/2001 dave
  !!    Tidied and indented
  !!   10:38, 16/05/2007 drb
  !!    Moved matrix_sum to stop error at end
  !!   2011/07/25  L.Tong
  !!     Implemented spin polarisation
  !!     added dependence on matL_dn, matH_dn, matLS, matSL, matL_dnS
  !!     and matSL_dn
  !!   2012/03/11 L.Tong
  !!   - Major rewrite for spin implementation
  !!   - implementation of spin non-polarised case and spin fixed-case
  !!     share same piece of code. Only implementation for the variable
  !!     spin case differs, due to need to use direct sum.
  !!  SOURCE
  !!
  subroutine McWMin(n_L_iterations, deltaE)

    use datatypes
    use numbers
    use matrix_data,   only: Lrange
    use mult_module,   only: allocate_temp_matrix, matrix_sum,         &
                             free_temp_matrix, matL, matH, matS, matT, &
                             symmetrise_L,  matLS, matSL
    use global_module, only: iprint_DM, IPRINT_TIME_THRES1,            &
                             IPRINT_TIME_THRES2,                       &
                             nspin, spin_factor,                       &
                             flag_fix_spin_population
    use timer_module,  only: cq_timer, start_timer, stop_print_timer,  &
                             WITH_LEVEL

    implicit none

    ! Passed variables
    integer,      intent(in)  :: n_L_iterations
    real(double), intent(out) :: deltaE
    ! Local variables
    integer                        :: n_iterations, spin
    integer,      dimension(nspin) :: matRhoNew, mat_oldL
    real(double), dimension(nspin) :: cn, oldE, omega1, c_old
    real(double)                   :: cn_tot, oldE_tot, omega1_tot, c_old_tot
    type(cq_timer)                 :: tmr_l_tmp1,tmr_l_tmp2
    logical,      dimension(nspin) :: done

    ! real(double) :: cn, oldE, omega1, c_old
    ! real(double) :: cn_up, cn_dn, oldE_up, oldE_dn, omega1_up, &
    !      omega1_dn, c_old_up, c_old_dn

    call start_timer(tmr_l_tmp1,WITH_LEVEL)

    done(:) = .false.
    oldE(:) = 1.0e30_double
    oldE_tot = spin_factor * sum(oldE(:))
    c_old(:) = zero
    c_old_tot = zero

    do spin = 1, nspin
       mat_oldL(spin) = allocate_temp_matrix(Lrange,0)
       matRhoNew(spin) = allocate_temp_matrix(Lrange,0)
       call matrix_sum(zero, mat_oldL(spin), one, matL(spin))
    end do

    if (nspin == 1 .or. flag_fix_spin_population) then

       iter: do n_iterations = 1, n_L_iterations

          call start_timer(tmr_l_tmp2, WITH_LEVEL)
          if (inode == ionode .and. iprint_DM >= 1) &
               write (io_lun, 1) n_iterations
          call McW_matrix_multiply_nonds(matH, matS, matL, matT, &
                                         matRhoNew, omega1, cn)

          do spin = 1, nspin
             if (omega1(spin) > oldE(spin) .or.                       &
                 abs(omega1(spin) - oldE(spin)) < 1.0e-08_double .or. &
                 cn(spin) > one .or.                                  &
                 cn(spin) < zero) then

                done(spin) = .true.

             end if
          end do
          if (done(1) .and. done(nspin)) then
             if (inode == ionode .and. iprint_DM > 0) then
                do spin = 1, nspin
                   write (io_lun, 2) spin, c_old(spin), oldE(spin),&
                                     cn(spin), omega1(spin)
                end do
             end if
             call stop_print_timer(tmr_l_tmp2, "a McWeeny iteration", &
                                   IPRINT_TIME_THRES2)
             exit iter
          end if

          deltaE = zero
          do spin = 1, nspin
              ! If new energy is lower, copy new rho into old rho and repeat
              deltaE = deltaE + spin_factor * (omega1(spin) - oldE(spin))
              call matrix_sum(zero, matL(spin), one, matRhoNew(spin))
              oldE(spin) = omega1(spin)
              c_old(spin) = cn(spin)
              if (inode == ionode .and. iprint_DM >= 0) &
                   write (io_lun, 3) spin, oldE(spin), c_old(spin)
           end do ! spin
           call symmetrise_L()
           call stop_print_timer(tmr_l_tmp2, "a McWeeny iteration",&
                                 IPRINT_TIME_THRES2)

        end do iter ! n_iterations

     else ! variable spin

        ! If the spin populations are allowed to vary there will be
        ! only one omega1_tot and cn_tot etc for the entire direct sum

        do n_iterations = 1, n_L_iterations
           call start_timer(tmr_l_tmp2, WITH_LEVEL)
           if (inode == ionode .and. iprint_DM >= 1) &
                write(io_lun, 1) n_iterations
           ! find new rho and energy
           call McW_matrix_multiply_ds(matH, matS, matL, matT, &
                                       matRhoNew, omega1_tot, cn_tot)
           if (omega1_tot > oldE_tot .or.                      &
               abs(omega1_tot - oldE_tot) < 1.0e-8_double .or. &
               cn_tot > one .or. cn_tot < zero) then

              if (inode == ionode .and. iprint_DM > 0) &
                   write (io_lun, 4) c_old_tot, oldE_tot, cn_tot, omega1_tot
              call stop_print_timer(tmr_l_tmp2, "a McWeeny iteration", &
                                    IPRINT_TIME_THRES2)
              exit

           end if
           ! If new energy is lower, copy new rho into old rho and repeat
           deltaE = omega1_tot - oldE_tot
           do spin = 1, nspin
              call matrix_sum(zero, matL(spin), one, matRhoNew(spin))
           end do
           oldE_tot = omega1_tot
           c_old_tot = cn_tot
           if (inode == ionode .and. iprint_DM >= 0) &
                write (io_lun, 5) oldE_tot, c_old_tot
           call symmetrise_L()
           call stop_print_timer(tmr_l_tmp2, "a McWeeny iteration", &
                                 IPRINT_TIME_THRES2)
        end do ! n_iterations

     end if ! (nspin == 1 .or. flag_fix_spin_population)

     do spin = nspin, 1, -1
        if (n_iterations < 3) &
             call matrix_sum(zero, matL(spin), one, mat_oldL(spin))
        call free_temp_matrix(matRhoNew(spin))
        call free_temp_matrix(mat_oldL(spin))
     end do

     if (inode == ionode .and. iprint_DM >= 1) then
        write (io_lun, 6) n_iterations
        if (nspin == 1 .or. flag_fix_spin_population) then
           do spin = 1, nspin
              write (io_lun, 7) spin, omega1(spin)
           end do
        else
           write (io_lun, 8) omega1_tot
        end if
     end if

     call stop_print_timer(tmr_l_tmp1,"MCWEENY",IPRINT_TIME_THRES1)

     return

1   format(/,20x,'McWeeny L iteration:',i5)
2   format(2x,'Rounding error found (spin=',i1,'): ',4f25.15)
3   format(2x,'Energy (spin=',i1,'): ',f25.15,' Mid-point: ',f25.15)
4   format(2x,'Rounding error found (direct sum): ',4f25.15)
5   format(2x,'Energy (direct sum): ',f25.15,' Mid-point: ',f25.15)
6   format(/,20x,'Functional value reached after ',i5,' L iterations')
7   format(20x,'Omega (spin=',i1,'): ',f25.15)
8   format(20x,'Omega (direct sum): ', f25.15)

  end subroutine McWMin
  !!***

  !!****f* McWeeny/InitMcW *
  !!
  !!  NAME
  !!   InitMcW
  !!  USAGE
  !!
  !!  PURPOSE
  !!   Initialises L to be a function of H with eigenvalues between 0 and 1
  !!  INPUTS
  !!
  !!
  !!  USES
  !!
  !!  AUTHOR
  !!   D.R.Bowler
  !!  CREATION DATE
  !!   End 10/98
  !!  MODIFICATION HISTORY
  !!   11/05/01 DRB
  !!    ROBODoc headers and better output also nsf from maxima_module
  !!   08/06/2001 dave
  !!    Tidying and added GenComms for gsum
  !!    Added my_barrier to use GenComms statement
  !!   15:25, 10/12/2002 drb
  !!    Bug fix - changed location of output for mubar and added gmin
  !!    and gmax calls on hmin and hmax
  !!   11:08, 04/02/2003 drb
  !!    Implemented quick fix for bracketing mu with Hamiltonian
  !!    limits (required for stability)
  !!   12:30, 2004/06/09 dave
  !!    Fixed assumption that T and S have same range
  !!   2004/10/29 drb
  !!    Removed gsum on SX: MAJOR bug found by TM which was messing up
  !!    initialisation; also removed loops to change hmin and hmax
  !!    (unnecessary, I think)
  !!   2011/07/25 L.Tong
  !!    Adding Spin polarisation
  !!   Friday, 2011/08/05 L.Tong
  !!    Removed the obsolete number_of_bands parameter, and the total
  !!    electron numbers are obtained directly from the global_module
  !!   2012/01/30 L.Tong
  !!    Rewrite for the variable spin case, works for the whole direct
  !!    sum.
  !!   2012/03/08 L.Tong
  !!    Rewrote spin implementation
  !!   2012/05/27 L.Tong
  !!   - Added check for the case ne(spin) > n_o(spin), i.e. the
  !!     number of electrons requested by the user at initialisation
  !!     in one spin channel exceeds the number of orbitals avaliable
  !!     for that spin channel. The end result of this is that
  !!     eigenvalues of the L will no longer be contained within 0 and
  !!     1, due to forced double occupancy of a same spin orbital. A
  !!     warning will be produced if this happens.
  !!  2013/03/06 L.Tong
  !!  - Added check for the case n_o == ne(spin) or n_e, in this case
  !!    mubar will be infinity, and cause the calculation to give
  !!    NaN. When n_o == ne(spin) or n_e, the initial L from
  !!    formulation should have been the identiy matrix (lambda = 0).
  !!  2018/11/15 nakata
  !!    Changed matS, matT, SX and A to be spin_SF dependent
  !!  2018/11/15 nakata
  !!    Removed direct_sum_factor and n_e and
  !!    introduced n_e_ds, n_o_ds and A_ds for variable spin case
  !!  SOURCE
  !!
  subroutine InitMcW

    use datatypes
    use numbers
    use global_module,  only: iprint_DM, ni_in_cell,                  &
                              IPRINT_TIME_THRES1, nspin, nspin_SF,    &
                              flag_fix_spin_population, ne_in_cell,   &
                              ne_spin_in_cell, spin_factor,           &
                              flag_SpinDependentSF
    use species_module, only: species, nsf_species
    use mult_module,    only: allocate_temp_matrix, matH, matS, matT, &
                              matL, mult, matrix_trace, T_S_TS,       &
                              matrix_product_trace, matrix_product,   &
                              free_temp_matrix
    use matrix_data,    only: Lrange, Srange, TSrange
    use GenBlas
    use GenComms,       only: gsum, my_barrier, gmin, gmax, cq_abort
    use timer_module,   only: cq_timer,start_timer, stop_print_timer, &
                              WITH_LEVEL

    implicit none

    ! Local variables
    real(double) :: n_e_ds, n_o, n_o_ds, hmax_ds, hmin_ds, SXHX_ds, A_ds

    real(double), dimension(nspin_SF) :: SX, A
    real(double), dimension(nspin) :: ne, mu1, mu2, mubar, lambda, &
                                      hmin, hmax, SXHX
    integer,      dimension(nspin) :: matXHX, matSXHX, mat_temp
    integer                        :: matTS, length, i, spin, spin_SF
    type(cq_timer)                 :: tmr_l_tmp1

    call start_timer(tmr_l_tmp1, WITH_LEVEL)
    if (inode == ionode .and. iprint_DM > 1) write (io_lun, 1)
    ! We must first initialise rho
    ! ne is the number of electrons per spin channel
    ne(1:nspin) = ne_spin_in_cell(1:nspin)
    ! n_e_ds is total number of electrons
    n_e_ds = spin_factor * sum(ne(:))

    ! n_o is number of support functions
    n_o = zero
    do i = 1, ni_in_cell
       n_o = n_o + nsf_species(species(i))
    end do
    ! we are working with the whole direct sum for variable spin case
    ! so we need to have twice the size of basis set
    n_o_ds = n_o + n_o

    ! check if we have enough orbitals for the given electron
    ! populations for the initialisation procedure to work properly.
    do spin = 1, nspin
       if (ne(spin) / n_o > one) then
          write (io_lun, '(/,2x,a,i1,a,i1,a,f12.6,a,f12.6,a,/)')      &
               '*** WARNING!!! You have too few number of support &
                &functions for the calculation. The number of &
                &electrons in spin = ', spin, ' is greater than the &
                &number of orbitals!. ne(', spin, ') = ', ne(spin),   &
                ', n_o for spin channel = ', n_o, &
                ' ***'
       end if
    end do

    if (nspin == 2) then
       if (inode == ionode .and. iprint_DM > 1) &
            write (io_lun, 3) ne(1), ne(2), n_o
    else
       if (inode == ionode .and. iprint_DM > 1) &
            write (io_lun, 2) n_e_ds, n_o_ds    ! check Miyazaki-san
    end if

    matTS = allocate_temp_matrix(TSrange,0)

    do spin_SF = 1, nspin_SF
       call matrix_product(matT(spin_SF), matS(spin_SF), matTS, mult(T_S_TS))
       SX(spin_SF) = matrix_trace(matTS)
    enddo
    if (inode == ionode .and. iprint_DM > 1) then
       if (nspin_SF == 1) then
          write (io_lun, fmt='(2x,"SX is    ",f25.15)') SX(1)
          write (io_lun, fmt='(2x,"SX/no is ",f25.15)') SX(1) / n_o
       else
          write (io_lun, fmt='(2x,"SX_up is    ",f25.15)') SX(1)
          write (io_lun, fmt='(2x,"SX_up/no is ",f25.15)') SX(1) / n_o
          write (io_lun, fmt='(2x,"SX_dn is    ",f25.15)') SX(2)
          write (io_lun, fmt='(2x,"SX_dn/no is ",f25.15)') SX(2) / n_o
       endif
    end if

    SXHX_ds = zero
    spin_SF = 1
    do spin = 1, nspin
       if (flag_SpinDependentSF) spin_SF = spin
       matXHX(spin) = allocate_temp_matrix(Lrange, 0)
       mat_temp(spin) = allocate_temp_matrix(Srange, 0)
       call McWXHX(matH(spin), matT(spin_SF), matXHX(spin), mat_temp(spin))
       SXHX(spin) = matrix_product_trace(matS(spin_SF), mat_temp(spin))
       SXHX_ds = SXHX_ds + spin_factor * SXHX(spin)
       if (inode == ionode .and. iprint_DM > 1) &
            write (io_lun, '(2x,"SXHX(spin=",i1,") = ",f25.15)') &
                  spin, SXHX(spin)
    end do
    if (inode == ionode .and. iprint_DM > 1) &
            write (io_lun, '(2x,"SXHX_ds = ",f25.15)') &
                  SXHX_ds

    call my_barrier()

    ! Get Gershgorin limits on H
    call GetHLimits(hmin, hmax)
    do spin = 1, nspin
       call gmin(hmin(spin))
       call gmax(hmax(spin))
       if (inode == ionode .and. iprint_DM > 1) &
          write (io_lun, 4) spin, hmin(spin), hmax(spin)
    end do
    ! for variable spin case find min and max of the direct sum
    if (nspin == 2 .and. (.not. flag_fix_spin_population)) then
       hmin_ds = minval(hmin(:))
       hmax_ds = maxval(hmax(:))
       if (inode == ionode .and. iprint_DM > 1) write (io_lun, 5) hmin, hmax
    end if

    do spin_SF = 1, nspin_SF
       A(spin_SF) = (one - SX(spin_SF)/n_o)
    enddo
    if (nspin == 2 .and. (.not. flag_fix_spin_population)) then
       if (nspin_SF==1) A_ds = (one - (SX(1)+SX(1))/n_o_ds)   ! SX(2) = SX(1)
       if (nspin_SF==2) A_ds = (one - (SX(1)+SX(2))/n_o_ds)
    endif

    if (nspin == 1 .or. flag_fix_spin_population) then
       spin_SF = 1
       do spin = 1, nspin
          if (flag_SpinDependentSF) spin_SF = spin
          ! for the rare case where n_o == ne(spin), we have lambda =
          ! 0, and mu_bar = infty. But initial rho will be identity
          ! matrix.
          if (n_o == ne(spin)) then
             mubar(spin) = zero
             lambda(spin) = zero
          else
             mu1(spin) = SXHX(spin) / n_o + hmax(spin) * A(spin_SF)
             mu2(spin) = (hmin(spin) * A(spin_SF) * ne(spin) / (n_o - ne(spin)) - &
                          SXHX(spin) / n_o) / &
                         (n_o * A(spin_SF) / (n_o - ne(spin)) - one)
             if (inode == ionode) &
                  write (io_lun, &
                         '(2x,"Mu1, Mu2, for spin = ",i1," are: ",2f25.15)') &
                         spin, mu1(spin), mu2(spin)
             if ((ne(spin) / (hmax(spin) - mu1(spin))) < &
                 ((n_o - ne(spin)) / (mu1(spin) - hmin(spin)))) then
                mubar(spin) = mu1(spin)
                lambda(spin) = ne(spin) / (hmax(spin) - mu1(spin))
             else if ((ne(spin) / (hmax(spin) - mu2(spin))) > &
                      ((n_o - ne(spin)) / (mu2(spin) - hmin(spin)))) then
                mubar(spin) = mu2(spin)
                lambda(spin) = (n_o - ne(spin)) / (mu2(spin) - hmin(spin))
             else
                call cq_abort('InitMcW: Cannot find mubar.')
             end if
          end if
          if (inode == ionode .and. iprint_DM > 1) &
               write (io_lun, &
                      '(2x,"mubar, lambda for spin = ",i1," are: ",2f25.15)') &
                     spin, mubar(spin), lambda(spin)
       end do
       ! Calculate L0.
       spin_SF = 1
       do spin = 1, nspin
          if (flag_SpinDependentSF) spin_SF = spin
          call McWRho0(matL(spin), matT(spin_SF), matXHX(spin), mubar(spin), &
               lambda(spin), ne(spin), n_o)
       end do
    else ! variable spin
       ! for the rare case where n_o_ds == n_e_ds, we have lambda =
       ! 0, and mu_bar = infty. But initial rho will be identity
       ! matrix.
       if (n_o_ds == n_e_ds) then
          mubar(:) = zero
          lambda(:) = zero
       else
          ! mu1(1) and mu1(2) always equal in this case. Same for mu2
          mu1(:) = SXHX_ds / n_o_ds + hmax_ds * A_ds
          mu2(:) = (hmin_ds * A_ds * n_e_ds / (n_o_ds - n_e_ds) - SXHX_ds / n_o_ds) / &
                   (n_o_ds * A_ds / (n_o_ds - n_e_ds) - one)
          if (inode == ionode) &
               write (io_lun, '(2x,"Mu1, Mu2: ",2f25.15)') mu1(1), mu2(1)
          if ((n_e_ds / (hmax_ds - mu1(1))) < &
              ((n_o_ds - n_e_ds) / (mu1(1) - hmin_ds))) then
             mubar(:) = mu1(:)
             lambda(:) = n_e_ds / (hmax_ds - mu1(:))
          else if ((n_e_ds / (hmax_ds - mu2(1))) < &
                   ((n_o_ds - n_e_ds) / (mu2(1) - hmin_ds))) then
             mubar(:) = mu2(:)
             lambda(:) = (n_o_ds - n_e_ds) / (mu2(:) - hmin_ds)
          else
             call cq_abort('InitMcW: Cannot find mubar.')
          end if
       end if
       if (inode == ionode .and. iprint_DM > 1) then
          write (io_lun, 6) mubar(1)
          write (io_lun, '(2x,"lambda is ",f25.15)') lambda(1)
       end if
       ! Calculate L0.
       spin_SF = 1
       do spin = 1, nspin
          if (flag_SpinDependentSF) spin_SF = spin
          call McWRho0(matL(spin), matT(spin_SF), matXHX(spin), mubar(spin), &
                       lambda(spin), n_e_ds, n_o_ds)
       end do
    end if

    ! Free the temporary matrices, order of deallocation is important
    do spin = nspin, 1, -1
       call free_temp_matrix(mat_temp(spin))
       call free_temp_matrix(matXHX(spin))
    end do
    call free_temp_matrix(matTS)

    call stop_print_timer(tmr_l_tmp1, "McWeeny initialisation", &
                          IPRINT_TIME_THRES1)

1   format(1x,'Welcome to InitMcW')
2   format(2x,'Electrons: ',f15.6,' Orbitals: ',f15.6)
3   format(2x,'Electrons_up: ',f15.6,' Electrons_dn: ',f15.6,' Orbitals for each spin: ',f15.6)
4   format(2x,'Minimum and maximum limites on H(spin=',i1,') are ',2f15.6)
5   format(2x,'Minimum and maximum limites on overall H are ',2f15.6)
6   format(2x,'Mubar is ',f15.6)

    return
  end subroutine InitMcW
  !!***

  !!****f* McWeeny/McW_matrix_multiply_nonds *
  !!
  !!  NAME
  !!   McW_matrix_multiply_nonds
  !!  USAGE
  !!
  !!  PURPOSE
  !!   Performs one iteration of McWeeny iterative scheme
  !!   Ref: PRB 58, 12704 (1998)
  !!  INPUTS
  !!
  !!  USES
  !!
  !!  AUTHOR
  !!   D.R.Bowler
  !!  CREATION DATE
  !!   21/10/98
  !!  MODIFICATION HISTORY
  !!   11/05/01 DRB
  !!    ROBODoc headers and better output also nsf from maxima_module
  !!   08/06/2001 dave
  !!    Added GenComms for gsum
  !!   2011/03/28 L.Tong
  !!    Removed local variable declaration matL2, redundant
  !!   2011/07/26 L.Tong
  !!    Implemented spin polarisation. For this subroutine, the spins
  !!    will be treated separately all together, and this subroutine
  !!    should be called twice involking different spin The matH, matL
  !!    and matRhoNew, and energy, c should be used accordingly for
  !!    different spin channels.
  !!    Decided for simplification, moving matLS and matSL refereces to
  !!    be part of the subroutine parameters instead of module
  !!    references. This tidies the code somewhat, in this case the
  !!    usual optional spin parameter is not required.
  !!    But still needs the flag_spin_polarisation flag to toggle the
  !!    factor of two for total energy calculation
  !!   2012/03/01 L.Tong
  !!    Changed name to McW_matrix_multiply_nonds, and to be called from
  !!    interface McW_matrix_multiply
  !!   2012/03/11 L.Tong
  !!    - Major rewrite of spin implementation
  !!    - Use matLS, matSL from mult module directly
  !!    - All inputs except matS and matInvS are arrays of dimension nspin
  !!   2018/11/15 nakata
  !!    Changed matS and matInvS to be spin_SF dependent
  !!  SOURCE
  !!
  subroutine McW_matrix_multiply_nonds(matH, matS, matL, matInvS, &
                                       matRhoNew, energy, c)
    use datatypes
    use numbers
    use primary_module, only: bundle
    use matrix_data,    only: LSLrange, Lrange, Hrange, Srange
    use mult_module,    only: matrix_product, matrix_sum,             &
                              matrix_product_trace, matrix_scale,     &
                              matrix_transpose, allocate_temp_matrix, &
                              free_temp_matrix, L_S_LS, LS_L_LSL,     &
                              LSL_SL_L, mult, matLS, matSL
    use global_module,  only: iprint_DM, nspin, spin_factor, flag_SpinDependentSF
    use GenComms,       only: gsum, cq_abort, my_barrier

    implicit none

    ! Passed variables
    integer,      dimension(:) :: matH, matS, matL, matInvS, matRhoNew
    real(double), dimension(:) :: energy, c
    ! Local Variables
    integer :: spin, spin_SF
    real(double), dimension(nspin) :: c1, c2, cn
    integer,      dimension(nspin) :: matLSL, matLSLSL, mat_top, mat_bottom

    spin_SF = 1

    do spin = 1, nspin
       if (flag_SpinDependentSF) spin_SF = spin
       call matrix_product(matL(spin), matS(spin_SF), matLS(spin), mult(L_S_LS))
       ! Generate LSL, LSLSL
       matLSL(spin) = allocate_temp_matrix(LSLrange,0)
       matLSLSL(spin) = allocate_temp_matrix(Lrange,0)
       mat_top(spin) = allocate_temp_matrix(Srange,0)
       mat_bottom(spin) = allocate_temp_matrix(Srange,0)
       call my_barrier
       call matrix_product(matLS(spin), matL(spin), matLSL(spin), &
                           mult(LS_L_LSL))
       call matrix_product(matLSL(spin), matLS(spin), matLSLSL(spin), &
                           mult(LSL_SL_L))
       !----------------------------------------------------------------
       ! Shorten L, LSL and LSLSL to range S and take traces (can do it
       ! with sdot because S is symmetric)
       !----------------------------------------------------------------
       call matrix_sum(zero, mat_top(spin),  one, matLSL(spin))
       call matrix_sum(one,  mat_top(spin), -one, matLSLSL(spin))
       c1(spin) = matrix_product_trace(matS(spin_SF), mat_top(spin))
       if (inode == ionode .and. iprint_DM >= 2) &
            write (io_lun, '(2x,"S.top (spin=",i1,") is  ",f25.15)') &
                  spin, c1(spin)
       call matrix_sum(zero, mat_bottom(spin), one, matL(spin))
       c1(spin) = matrix_product_trace(matS(spin_SF), mat_bottom(spin))
       if (inode == ionode .and. iprint_DM >= 2) &
            write (io_lun, '(2x,"N_e (spin=",i1,") is  ",f25.15)') &
                  spin, c1(spin)
       call matrix_sum(one, mat_bottom(spin), -one, matLSL(spin))
       c1(spin) = matrix_product_trace(matS(spin_SF), mat_top(spin))
       c2(spin) = matrix_product_trace(matS(spin_SF), mat_bottom(spin))
       if (c2(spin) /= zero) then
          cn(spin) = c1(spin) / c2(spin)
       else
          call cq_abort('McW_matrix_multiply: c2 is zero for spin', spin)
       endif
       c(spin) = cn(spin)
       if (inode == ionode .and. iprint_DM >= 2) &
            write (io_lun, '(2x,"c, c1, c2 (spin=",i1,") are ",3f25.15)') &
                  spin, cn(spin), c1(spin), c2(spin)
       ! Shorten LSL and LSLSL to range L
       call matrix_scale(zero, matRhoNew(spin))
       call matrix_sum(zero, matRhoNew(spin), -one, matLSLSL(spin))
       call matrix_sum(one, matRhoNew(spin), one + cn(spin), matLSL(spin))
       if(cn(spin) < half) then
          call matrix_sum(one, matRhoNew(spin), (one - two * &
                          cn(spin)), matL(spin))
          call matrix_scale((one / (one - cn(spin))), matRhoNew(spin))
       else
          call matrix_scale(one / cn(spin), matRhoNew(spin))
       endif

       ! Idempotency checks
       !%%! data_L3 = zero
       !%%! call axpy(NSF*NSF*mat(1,Lrange)%length,one,dat_L,1,data_L3,1)
       !%%! ! LSLSL -> 2LSLSL - 3LSL
       !%%! call matrix_add(2.0_double, data_LSLSL, mat(:,Lrange), &
       !%%!      -3.0_double, data_LSL, mat(:,LSLrange), bundle, inode)
       !%%! ! Lnew -> Lnew+(2LSLSL-3LSL) [= Lnew - McW(Lold)]
       !%%! call axpy(NSF*NSF*mat(1,Lrange)%length,one,data_LSLSL,1,data_L3,1)
       !%%! call mult_wrap(inode-1,data_L3,dat_S,data_LS,mult(L_S_LS))
       !%%! call mult_wrap(inode-1,dat_S,data_SL,data_XHX,mult(S_LS_L))
       !%%! cn = dot(NSF*NSF*mat(1,Lrange)%length,data_SLS2, 1, data_L3, 1)
       !%%! call gsum(cn)
       !%%! if(inode.eq.ionode.AND.iprint_DM>=2) &
       !%%!     write(io_lun,*) 'Idempot monitor is ',cn/256.0_double

       ! Shorten rhonew to range H into temp
       call free_temp_matrix(mat_bottom(spin))
       mat_bottom(spin) = allocate_temp_matrix(Hrange, 0)
       call matrix_sum(zero, mat_bottom(spin), one, matRhoNew(spin))
       energy(spin) = matrix_product_trace(mat_bottom(spin), matH(spin))
       if (inode == ionode .and. iprint_DM >= 3) &
            write (io_lun, '(2x,"energy (spin=",i1,") is ",f25.15)') &
                  spin, energy(spin)
       call matrix_sum(zero, mat_top(spin), one, matRhoNew(spin))
       c1(spin) = matrix_product_trace(mat_top(spin), matS(spin_SF))
       if (inode == ionode .and. iprint_DM >= 3) &
            write (io_lun, '(2x,"N_e(2) (spin=",i1,") is ",f25.15)') &
                  spin, c1(spin)
       ! note that N_e is printed without factor of two here even for
       ! spin un-polarised results.
       call free_temp_matrix(mat_bottom(spin))
       call free_temp_matrix(mat_top(spin))
       call free_temp_matrix(matLSLSL(spin))
       call free_temp_matrix(matLSL(spin))
       if (abs(c2(spin)) < 1.0e-6) c(spin) = -0.01_double
    end do ! spin

    return

  end subroutine McW_matrix_multiply_nonds
  !!***

  !!****f* McWeeny/McW_matrix_multiply_ds *
  !! PURPOSE
  !!   Performs one iteration of McWeeny iterative scheme treating the
  !!   overall matrices L, H as direct sums of their spin
  !!   components. This allows population migration from one spin
  !!   channel to another
  !!   Ref: PRB 58, 12704 (1998)
  !! INPUTS
  !! OUTPUT
  !! RETURN VALUE
  !! AUTHOR
  !!   L.Tong
  !! CREATION DATE
  !!   2011/07/26
  !! MODIFICATION HISTORY
  !!   Tuesday, 2011/08/02 L.Tong
  !!     Changed the order of input parameters, replaced energy_up and
  !!     energy_down with energy_total, and c with c_total
  !!   2012/03/01 L.Tong
  !!     Renamed to McW_matrix_multiply_ds, to be called by interface
  !!     McW_matrix_multiply
  !!   2012/03/11 L.Tong
  !!   - Major rewrite of the spin implementation
  !!   - all inputs except matS and matInvS are arrays of dimension nspin
  !!   2018/11/13 17:30 nakata
  !!    Changed matS and matInvS to be spin_SF dependent
  !! SOURCE
  !!
  subroutine McW_matrix_multiply_ds(matH, matS, matL, matInvS, &
                                    matRhoNew, energy_total, c_total)

    use datatypes
    use numbers
    use primary_module, only: bundle
    use matrix_data,    only: LSLrange, Lrange, Hrange, Srange
    use mult_module,    only: matrix_product, matrix_sum,             &
                              matrix_product_trace, matrix_scale,     &
                              matrix_transpose, allocate_temp_matrix, &
                              free_temp_matrix, matLS, matSL, L_S_LS, &
                              LS_L_LSL, LSL_SL_L, mult
    use global_module,  only: iprint_DM, nspin, spin_factor, flag_SpinDependentSF
    use GenComms,       only: gsum, cq_abort, my_barrier

    implicit none

    ! Passed variables
    integer, dimension(:) :: matH, matS, matL, matInvS, matRhoNew
    real(double)          :: energy_total, c_total

    ! Local Variables
    real(double)              :: c1, c2, tmp
    integer, dimension(nspin) :: matLSL, matLSLSL, mat_top, mat_bottom
    integer                   :: spin, spin_SF

    c1 = zero
    c2 = zero
    energy_total = zero
    spin_SF = 1

    do spin = 1, nspin
       if (flag_SpinDependentSF) spin_SF = spin
       call matrix_product(matL(spin), matS(spin_SF), matLS(spin), mult(L_S_LS))
       ! Generate LSL, LSLSL
       matLSL(spin) = allocate_temp_matrix(LSLrange, 0)
       matLSLSL(spin) = allocate_temp_matrix(Lrange, 0)
       mat_top(spin) = allocate_temp_matrix(Srange, 0)
       mat_bottom(spin) = allocate_temp_matrix(Srange, 0)
       call my_barrier
       call matrix_product(matLS(spin), matL(spin), matLSL(spin), &
                           mult(LS_L_LSL))
       call matrix_product(matLSL(spin), matLS(spin), matLSLSL(spin), &
                           mult(LSL_SL_L))
       !----------------------------------------------------------------
       ! Shorten L, LSL and LSLSL to range S and take traces (can do it
       ! with sdot because S is symmetric)
       !----------------------------------------------------------------
       call matrix_sum(zero, mat_top(spin),  one, matLSL(spin))
       call matrix_sum(one,  mat_top(spin), -one, matLSLSL(spin))
       tmp = matrix_product_trace(matS(spin_SF), mat_top(spin))
       if (inode == ionode .and. iprint_DM >= 2) &
            write (io_lun, '(2x,"S.top (spin=",i1,") is ",f25.15)') spin, tmp
       call matrix_sum(zero, mat_bottom(spin), one, matL(spin))
       tmp = matrix_product_trace(matS(spin_SF), mat_bottom(spin))
       if (inode == ionode .and. iprint_DM >= 2) &
            write (io_lun, '(2x,"N_e (spin=",i1,") is ",f25.15)') spin, tmp
       ! get the real mat_bottom for spin up
       call matrix_sum(one, mat_bottom(spin), -one, matLSL(spin))
       ! c1 and c2 calculated
       c1 = c1 + spin_factor * matrix_product_trace(matS(spin_SF), mat_top(spin))
       c2 = c2 + spin_factor * matrix_product_trace(matS(spin_SF), mat_bottom(spin))
       ! finished  mat_bottom here, change to Hrange for something else
       call free_temp_matrix(mat_bottom(spin))
       mat_bottom(spin) = allocate_temp_matrix(Hrange,0)
    end do ! spin
    if (c2 /= zero) then
       c_total = c1 / c2
    else
       call cq_abort('McW_matrix_multiply_ds: c2 is zero')
    end if
    if (inode == ionode .and. iprint_DM >= 2) &
         write (io_lun, '(2x,"c_total, c1, c2 are ",3f25.15)') &
               c_total, c1, c2
    ! Shorten LSL and LSLSL to range L
    do spin = 1, nspin
       call matrix_scale(zero, matRhoNew(spin))
       call matrix_sum(zero, matRhoNew(spin), -one, matLSLSL(spin))
       call matrix_sum(one, matRhoNew(spin), one + c_total, matLSL(spin))
       if (c_total < half) then
          call matrix_sum(one, matRhoNew(spin), (one - two * c_total),&
                          matL(spin))
          call matrix_scale((one/(one - c_total)), matRhoNew(spin))
       else
          call matrix_scale(one/c_total, matRhoNew(spin))
       end if
       call matrix_sum(zero, mat_bottom(spin), one, matRhoNew(spin))
       energy_total = energy_total + spin_factor * &
                      matrix_product_trace(mat_bottom(spin), matH(spin))
    end do ! spin
    if (inode == ionode .and. iprint_DM >= 3) &
         write (io_lun, '(2x,"energy_total: ",f25.15)') energy_total
    do spin = 1, nspin
       call matrix_sum(zero, mat_top(spin), one, matRhoNew(spin))
       tmp = matrix_product_trace(mat_top(spin), matS(spin_SF))
       if (inode == ionode .and. iprint_DM >= 3) &
            write (io_lun, '(2x,"N_e (spin=",i1,") is ",f25.15)') &
                  spin, tmp
    end do ! spin

    ! free temp matrices
    do spin = nspin, 1, -1
       call free_temp_matrix(mat_bottom(spin))
       call free_temp_matrix(mat_top(spin))
       call free_temp_matrix(matLSLSL(spin))
       call free_temp_matrix(matLSL(spin))
    end do

    return

  End Subroutine McW_matrix_multiply_ds
  !!***

  !!****f* MacWeeny/McWXHX
  !! PURPOSE
  !!  sbrt McWXHX: Pre- and post- multiplies H by S^-1, and returns
  !!  the result to both L range (bab) and S range (tm)
  !! INPUTS
  !! OUTPUT
  !! RETURN VALUE
  !! AUTHOR
  !!   drb
  !! CREATION DATE
  !!
  !! MODIFICATION HISTORY
  !! SOURCE
  !!
  subroutine McWXHX(matA, matB, matBAB, mat_temp)

    use datatypes
    use numbers
    use matrix_data,   only: THrange
    use mult_module,   only: matrix_product, matrix_sum, T_H_TH, &
                             TH_T_L, allocate_temp_matrix, mult, &
                             free_temp_matrix, matrix_trace
    use global_module, only: iprint_DM

    implicit none

    ! Passed variables
    integer :: matA, matB, matBAB, mat_temp

    ! Local variable
    integer :: matBA
    real(double) :: xx

    matBA = allocate_temp_matrix(THrange,0)
    call matrix_product(matB, matA, matBA, mult(T_H_TH))
    if (iprint_DM >= 2) xx = matrix_trace(matBA)
    if (inode == ionode .and. iprint_DM >= 2) &
         write (io_lun, '(2x,"Trace of BA:            ",f25.15)') xx
    call matrix_product(matBA, matB, matBAB, mult(TH_T_L))
    if (iprint_DM >= 2) xx = matrix_trace(matBAB)
    if (inode == ionode .and. iprint_DM >= 2) &
         write (io_lun, '(2x,"Trace of BAB (L range): ",f25.15)') xx
    call matrix_sum(zero, mat_temp, one, matBAB)
    if (iprint_DM >= 2) xx = matrix_trace(mat_temp)
    if (inode == ionode .and. iprint_DM >= 2) &
         write (io_lun, '(2x,"Trace of BAB (S range): ",f25.15)') xx
    call free_temp_matrix(matBA)
    return
  end Subroutine McWXHX
  !!***

  !!****f*
  !! PURPOSE
  !! INPUTS
  !! OUTPUT
  !! RETURN VALUE
  !! AUTHOR
  !!   drb
  !! CREATION DATE
  !! MODIFICATION HISTORY
  !! SOURCE
  !!
  subroutine McWRho0(matA, matB, matC, m, l, n_e, n_o)

    use datatypes
    use numbers
    use mult_module,   only: matrix_sum
    use global_module, only: iprint_DM

    implicit none

    integer :: matA, matB, matC
    Real( double ) :: m, l, n_e, n_o, tmp

    tmp = (l*m+n_e)/n_o
    if (inode == ionode .and. iprint_DM >= 2) &
         write (io_lun, '(2x,"l, m, n_e, n_o and tmp are ",5f25.15)') &
               l, m, n_e, n_o, tmp
    call matrix_sum(zero, matA, tmp, matB)
    tmp = -l/n_o
    if (inode == ionode .and. iprint_DM >= 2) &
         write (io_lun, '(2x,"tmp is ",f25.15)') tmp
    call matrix_sum(one, matA, tmp, matC)
    return
  end Subroutine McWRho0
  !!***

  !!****f*
  !! PURPOSE
  !!   Finds the Gershgorin limits on the eigenvalues of the Hamiltonian
  !! INPUTS
  !! OUTPUT
  !! RETURN VALUE
  !! AUTHOR
  !!   drb
  !! CREATION DATE
  !! MODIFICATION HISTORY
  !!   2011/07/25 L.Tong
  !!     Implemented Spin Polarisation support: optional variable spin, that
  !!     makes the subroutine to calculate spin up of spin down components.
  !!   Wednesday, 2011/08/03 L.Tong
  !!     Fixed the bug caused by the conditional if (present (spin)
  !!     .and. spin == 2), on some compilers which does not short circuit
  !!     and evaluates both sides of .and. and when the subroutine is
  !!     called without the optional spin parameter, then spin == 2 gives
  !!     a segmentation fault. Fixed this by introducing logical variable
  !!     do_spin_down, we work out this variable correctly first, and
  !!     then apply conditional on it in place of the previous present
  !!     (spin) .and. spin == 2 conditional. Also moved the conditional
  !!     outside of the inner loops to improve efficiency
  !!   2012/03/11 L.Tong
  !!   - Major rewrite of spin implementation
  !!   - removed optional parameter spin
  !!   - MinH and MaxH are now arrays of spin, each contain min and max
  !!     of H in each spin channels.
  !! SOURCE
  !!
  subroutine GetHLimits(MinH, MaxH)

    use datatypes
    use numbers
    use matrix_module,  only: matrix, matrix_halo
    use matrix_data,    only: mat, Hrange, halo
    use cover_module,   only: BCS_parts
    use primary_module, only: bundle
    use mult_module,    only: matH, return_matrix_value, matrix_pos
    use GenComms,       only: cq_abort
    use maxima_module,  only: maxnsf
    use global_module,  only: nspin
    use memory_module,  only: reg_alloc_mem, reg_dealloc_mem, type_dbl

    implicit none

    ! Passed variables
    real(double), dimension(:) :: MinH, MaxH

    ! Local variables
    integer :: np, nat, nb, ist, i, j, gcspart, wheremat, iprim, nsfi,&
               nsfj, stat, spin
    real(double), dimension(:,:), allocatable :: HMin, HMax
    real(double) :: Hval

    allocate(HMin(maxnsf,nspin), HMax(maxnsf,nspin), STAT=stat)
    if (stat /= 0) &
         call cq_abort("GetHLimits: Error allocating mem: ", maxnsf, nspin)
    call reg_alloc_mem(area_DM, 2*maxnsf*nspin, type_dbl)

    MinH(:) = 1.0e30_double
    MaxH(:) = -1.0e30_double

    call start_timer(tmr_std_matrices)
    do spin = 1, nspin
       iprim = 0
       do np = 1, bundle%groups_on_node
          if (bundle%nm_nodgroup(np) > 0) then
             do nat = 1,bundle%nm_nodgroup(np)
                iprim=iprim+1
                nsfi = mat(np,Hrange)%ndimi(nat)
                HMin(:,spin) = zero
                HMax(:,spin) = zero
                do nb = 1, mat(np,Hrange)%n_nab(nat)
                   ist = mat(np,Hrange)%i_acc(nat) + nb - 1
                   nsfj = mat(np,Hrange)%ndimj(ist)
                   gcspart = &
                        BCS_parts%icover_ibeg(mat(np,Hrange)%i_part(ist)) +&
                        mat(np,Hrange)%i_seq(ist) - 1
                   wheremat = matrix_pos(matH(spin), iprim, &
                                         halo(Hrange)%i_halo(gcspart))
                   if (wheremat == mat(np, Hrange)%onsite(nat)) then
                      do i = 1, nsfj
                         do j = 1, nsfi
                            Hval = return_matrix_value(matH(spin), np, &
                                                       nat, iprim, nb, j, i)
                            HMin(i,spin) = HMin(i,spin) + abs(Hval)
                            HMax(i,spin) = HMax(i,spin) + abs(Hval)
                         end do
                      end do
                   else
                      do i = 1,nsfj
                         do j= 1, nsfi
                            Hval = return_matrix_value(matH(spin), np, &
                                                       nat, iprim, nb, j, i)
                            HMin(i,spin) = HMin(i,spin) - abs(Hval)
                            HMax(i,spin) = HMax(i,spin) + abs(Hval)
                         end do
                      end do
                   end if
                end do ! do nb = 1, mat(np,Hrange)%n_nab(nat)
                ! do a simple search for the min amd max among HMin(i) and HMax(i)
                do i = 1, nsfi
                   if (HMin(i,spin) < MinH(spin)) MinH(spin) = HMin(i,spin)
                   if (HMax(i,spin) > MaxH(spin)) MaxH(spin) = HMax(i,spin)
                enddo
             end do ! nat = bundle%nm_nodgroup
          end if ! (bundle%nm_nodgroup(np) > 0)
       end do ! np
    end do ! spin

    call stop_timer(tmr_std_matrices)

    deallocate(HMin, HMax, STAT=stat)
    if (stat /= 0) &
         call cq_abort("GetHLimits: Error deallocating mem")
    call reg_dealloc_mem(area_DM, 2*maxnsf*nspin, type_dbl)

    return
  end subroutine GetHLimits
  !!***

end module McWeeny

max_stars_repo_name: negf/CONQUEST-release

max_issues_repo_name: negf/CONQUEST-release

max_forks_repo_name: negf/CONQUEST-release


text:
! { dg-do run }
! Tests the fix for PR40591 in which the interface 'sub2'
! for 'pptr2' was not resolved.
!
! Contributed by Tobias Burnus  <burnus@gcc.gnu.org>
!
program main
  call test
contains
  subroutine sub1(arg)
    integer arg
    arg = arg + 1
  end subroutine sub1
  subroutine test()
    procedure(sub1), pointer :: pptr1
    procedure(sub2), pointer :: pptr2
    integer i
    i = 0
    pptr1 => sub1
    call pptr1 (i)
    pptr1 => sub2
    call pptr1 (i)
    pptr2 => sub1
    call pptr2 (i)
    pptr2 => sub2
    call pptr2 (i)
    if (i .ne. 22) call abort
  end subroutine test
  subroutine sub2(arg)
    integer arg
    arg = arg + 10
  end subroutine sub2
end program main

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab