module interpolation
use phys_consts, only: wp

implicit none (type, external)
public

interface ! interp2d.f90
module pure function interp2_plaid(x1,x2,f,x1i,x2i)
real(wp), dimension(:), intent(in) :: x1, x2, x1i, x2i
real(wp), dimension(:,:), intent(in) :: f
real(wp) :: interp2_plaid(1:size(x1i,1),1:size(x2i,1))
end function interp2_plaid

module pure function interp2(x1,x2,f,x1i,x2i)
real(wp), dimension(:), intent(in) :: x1, x2, x1i, x2i
real(wp), dimension(:,:), intent(in) :: f
real(wp) :: interp2(1:size(x1i,1))
end function interp2
end interface

contains

pure real(wp) function interp1(x1,f,x1i)

!------------------------------------------------------------
!-------A 1D LINEAR INTERPOLATION FUNCTION.  THE INDEPENDENT
!-------VARIABLE FOR THE GIVEN DATA GRID MUST BE MONOTONICALLY
!-------INCREASING, BUT THE LOCATIONS FOR THE INTERPOLATION DATA
!-------CAN BE IN ANY ORDER.  NOTE THAT EXTRAPOLATION IS NOT DONE
!-------AT ALL BY THIS FUNCTION - VALUES OUTSIDE DOMAIN OF
!-------ORIGINAL DATA ARE SIMPLY SET TO ZERO.
!------------------------------------------------------------

real(wp), dimension(:), intent(in) :: x1,f, x1i
dimension :: interp1(1:size(x1i,1))

integer :: lx1,lx1i,ix1,ix1i
real(wp) :: slope

integer :: ix10,ix1fin

lx1=size(x1,1)
lx1i=size(x1i,1)

do ix1i=1,lx1i
!      !find the 'bin' for this point; i.e. find ix1 s.t. xi(ix1i) is between x(ix1-1) and x(ix1)
  ix10=1
  ix1=lx1/2
  ix1fin=lx1
  if (x1i(ix1i)>=x1(1) .and. x1i(ix1i)<=x1(lx1)) then    !in bounds
    do while(.not.(x1i(ix1i)>=x1(ix1-1) .and. x1i(ix1i)<=x1(ix1)))    !keep going until we are in the interval we want
      if (x1i(ix1i)>=x1(ix10) .and. x1i(ix1i)<=x1(ix1)) then    !left half (correct guess)
        ix1fin=ix1
      else    !wrong take the "right" half (har har)
        ix10=ix1
      end if
      ix1=(ix1fin+ix10)/2
      if (ix10==ix1) then
        ix1=lx1
      end if
    end do
  else if (x1i(ix1i)<x1(1)) then
    ix1=1
  else
    ix1=lx1
  end if

  !execute interpolation for this point
  if (ix1>1 .and. ix1<=lx1) then   !interpolation
    slope=(f(ix1)-f(ix1-1))/(x1(ix1)-x1(ix1-1))
    interp1(ix1i)=f(ix1-1)+slope*(x1i(ix1i)-x1(ix1-1))
  else
    interp1(ix1i)=0
  end if
end do

!THERE IS SOME ISSUE WITH POINTS OUTSIDE INTERPOLANT DOMAIN - THIS IS A
!WORKAROUND UNTIL I CAN PIN DOWN THE EXACT PROBLEM
do ix1i=1,lx1i
  if(x1i(ix1i)<x1(1) .or. x1i(ix1i)>x1(lx1)) then
    interp1(ix1i)=0
  end if
end do

end function interp1

pure real(wp) function interp3(x1,x2,x3,f,x1i,x2i,x3i)

!------------------------------------------------------------
!-------A 2D BILINEAR INTERPOLATION FUNCTION.  THIS VERSION ASSUMES
!-------THAT THE LIST OF OUTPUT POINTS IS A 'FLAT LIST' RATHER THAN
!-------DESCRIPTIVE OF A 2D MESHGRID.
!------------------------------------------------------------

real(wp), dimension(:), intent(in) :: x1,x2,x3,x1i,x2i,x3i
real(wp), dimension(:,:,:), intent(in) :: f
dimension :: interp3(1:size(x1i,1))     !interpolated points are a flat list

real(wp) :: fx1ix2pix3p,fx1ix2nix3p,fx1ix2pix3n,fx1ix2nix3n    !function estimates at x1i point  vs. at x2 interfaces
real(wp) :: fx2ix3p,fx2ix3n
real(wp) :: slope                     !temp value for slope for interpolations

integer :: lx1,lx2,lx3,lxi,ix1,ix2,ix3,ixi
integer :: ix10,ix1fin,ix20,ix2fin,ix30,ix3fin

lx1=size(x1,1)
lx2=size(x2,1)
lx3=size(x3,1)
lxi=size(x1i,1)    !only one size since this a flat list of grid points

do ixi=1,lxi
  !find the x1 'bin' for this point; i.e. find ix1 s.t. xi(ix1i) is between x(ix1-1) and x(ix1)
  ix10=1
  ix1=lx1/2
  ix1fin=lx1
  if (x1i(ixi)>=x1(1) .and. x1i(ixi)<=x1(lx1)) then    !in bounds
    do while(.not.(x1i(ixi)>=x1(ix1-1) .and. x1i(ixi)<=x1(ix1)))    !keep going until we are in the interval we want
      if (x1i(ixi)>=x1(ix10) .and. x1i(ixi)<=x1(ix1)) then    !left half (correct guess)
        ix1fin=ix1
      else    !wrong take the "right" half (har har)
        ix10=ix1
      end if
      ix1=(ix1fin+ix10)/2
      if (ix10==ix1) then
        ix1=lx1
      end if
    end do
  else if (x1i(ixi)<x1(1)) then
    ix1=1
  else
    ix1=lx1
  end if

  !find the x2 'bin' for this point; i.e. find ix2 s.t. x2i(ix2i) is between x2(ix2-1) and x2(ix2)
  ix20=1
  ix2=lx2/2
  ix2fin=lx2
  if (x2i(ixi)>=x2(1) .and. x2i(ixi)<=x2(lx2)) then    !in bounds
    do while(.not.(x2i(ixi)>=x2(ix2-1) .and. x2i(ixi)<=x2(ix2)))    !keep going until we are in the interval we want
      if (x2i(ixi)>=x2(ix20) .and. x2i(ixi)<=x2(ix2)) then    !left half (correct guess)
        ix2fin=ix2
      else    !wrong take the "right" half (har har)
        ix20=ix2
      end if
      ix2=(ix2fin+ix20)/2
      if (ix20==ix2) then
        ix2=lx2
      end if
    end do
  else if (x2i(ixi)<x2(1)) then
    ix2=1
  else
    ix2=lx2
  end if

  !find the x3 'bin' for this point; i.e. find ix3 s.t. x3i(ix3i) is between x3(ix3-1) and x3(ix3)
  ix30=1
  ix3=lx3/2
  ix3fin=lx3
  if (x3i(ixi)>=x3(1) .and. x3i(ixi)<=x3(lx3)) then    !in bounds
    do while(.not.(x3i(ixi)>=x3(ix3-1) .and. x3i(ixi)<=x3(ix3)))    !keep going until we are in the interval we want
      if (x3i(ixi)>=x3(ix30) .and. x3i(ixi)<=x3(ix3)) then    !left half (correct guess)
        ix3fin=ix3
      else    !wrong take the "right" half (har har)
        ix30=ix3
      end if
      ix3=(ix3fin+ix30)/2
      if (ix30==ix3) then
        ix3=lx3
      end if
    end do
  else if (x3i(ixi)<x3(1)) then
    ix3=1
  else
    ix3=lx3
  end if

  if (ix1>1 .and. ix1<=lx1 .and. ix2>1 .and. ix2<=lx2 .and. ix3>1 .and. ix3<=lx3) then   !interpolation
    !interpolate x1 for fixed values of x2,x3 (four separate interps)
    !first the "prev" x2 value, "prev" x2 value
    slope=(f(ix1,ix2-1,ix3-1)-f(ix1-1,ix2-1,ix3-1))/(x1(ix1)-x1(ix1-1))
    fx1ix2pix3p=f(ix1-1,ix2-1,ix3-1)+slope*(x1i(ixi)-x1(ix1-1))

    !now the "next" x2 value, "prev" x3
    slope=(f(ix1,ix2,ix3-1)-f(ix1-1,ix2,ix3-1))/(x1(ix1)-x1(ix1-1))
    fx1ix2nix3p=f(ix1-1,ix2,ix3-1)+slope*(x1i(ixi)-x1(ix1-1))

    !prev x2, next x3
    slope=(f(ix1,ix2-1,ix3)-f(ix1-1,ix2-1,ix3))/(x1(ix1)-x1(ix1-1))
    fx1ix2pix3n=f(ix1-1,ix2-1,ix3)+slope*(x1i(ixi)-x1(ix1-1))

    !next x3, next x3
    slope=(f(ix1,ix2,ix3)-f(ix1-1,ix2,ix3))/(x1(ix1)-x1(ix1-1))
    fx1ix2nix3n=f(ix1-1,ix2,ix3)+slope*(x1i(ixi)-x1(ix1-1))

    !interpolate between each x2 value (two separate interps)
    !interp in x2 for the x3 prev points
    slope=(fx1ix2nix3p-fx1ix2pix3p)/(x2(ix2)-x2(ix2-1))
    fx2ix3p=fx1ix2pix3p+slope*(x2i(ixi)-x2(ix2-1))

    !interp in 2 for the next x3 points
    slope=(fx1ix2nix3n-fx1ix2pix3n)/(x2(ix2)-x2(ix2-1))
    fx2ix3n=fx1ix2pix3n+slope*(x2i(ixi)-x2(ix2-1))

    !finally an interpolation in x2 to finish things off (single interp)
    slope=(fx2ix3n-fx2ix3p)/(x3(ix3)-x3(ix3-1))
    interp3(ixi)=fx2ix3p+slope*(x3i(ixi)-x3(ix3-1))
  else
    interp3(ixi)=0._wp
  end if
end do

!THERE IS SOME ISSUE WITH POINTS OUTSIDE INTERPOLANT DOMAIN - THIS IS A WORKAROUND UNTIL I CAN PIN DOWN THE EXACT PROBLEM
do ixi=1,lxi
if(x1i(ixi)<x1(1) .or. x1i(ixi)>x1(lx1) .or. x2i(ixi)<x2(1) .or. x2i(ixi)>x2(lx2) .or. x3i(ixi)<x3(1) .or. &
   x3i(ixi)>x3(lx3) .or. lx1==1 .or. lx2==1 .or. lx3==1) then     !also cover the case where there is a singleton dimension...
  interp3(ixi)=0._wp
end if
end do

end function interp3

end module interpolation


text:
      subroutine invcpfft2L2(psi,nr1L2,nr2L2,nr3L2)

c     Written by A. Canning (CRAY-EPFL) 25th July 94
c
c     output = psi(z,y,x)  wavefunction in real space z columns
c     input  = psi(x,y,z)  x columns of g vectors in load balalanced distr.
c
c     computes inverse fft specifically for CP algo
c     ie taking sphere and going to cylinder then two (one for gamma)
c     slices then a complete cube.
c     see fftprep.f for more details
c     fftprep must be called once before this routine
c     for setups.
c
c     last revised ADV 7/7/95 to parametrise the no. of
c     chunks on each PE
c     AMC revised 21/6/97 hardwired for gamma point complex to real fft

      use load_data
      use fft_data
      use data

      implicit none

      include "mpif.h"
c
c
      real*8,allocatable,dimension(:) :: worknr1,worknr2,worknr3
      real*8,allocatable,dimension(:) :: psiy,combuf1,combuf2

      integer inode,nnodes
      integer inode_tot,nnodes_tot,icolor,num_group,MPI_COMM_K,
     &       MPI_COMM_N

      integer ireq(nnodes)

      integer jj

      common /mpi_data/inode,nnodes,inode_tot,nnodes_tot,
     &  icolor,num_group,MPI_COMM_K,MPI_COMM_N

      integer mpistatus(mpi_status_size)

c     scalars used

      integer i,ib,ic,idum,ii,ilocadd,isc,isign,itar,itaradd,
     c     itnode,iw,ix,iy,izb,j,jcol,ngy,ngyadd,ngz,nr1L2,
     c     nr2L2,nr3L2,nr3uL2,iloc_dum,ierr
      integer nworknr1,nworknr2,nworknr3

      real*8 psi(mr_nL2)
      real*8 psidum(nr3L2),sum

      nworknr1 = 20000+2.28*nr1xL2
      nworknr2 = 20000+2.28*nr2xL2
      nworknr3 = 20000+2.28*nr3xL2

cccc

      allocate(worknr1(nworknr1))
      allocate(worknr2(nworknr2))
      allocate(worknr3(nworknr3))
      allocate(psiy(mr_nL2))
      allocate(combuf1(mr_nL2))
      allocate(combuf2(mr_nL2))

      call mpi_barrier(MPI_COMM_K,ierr)

c
c     for gamma point fft
c
      nr3uL2 = nr3L2/2 + 1
c
c     memory allocation for local arrays used with SHMEM routines
c
c
      isign = -1
c
c     do the FFT's in place on the non-zero columns
c     of psi
c     so this is the x dir FFT
c

      do jcol = 1,ncol2L2(inode)
         ilocadd = 1+(jcol-1)*(2*nr1L2)

         call system_ccfft(0,isign,nr1L2,1.0d0,psi(ilocadd),
     &    psi(ilocadd),tabnr1linL2,worknr1,0,ntabnr1L2,nworknr1)

      enddo                     !i

c
c     transpose to the slice  y col mode  in psiy
c
      do ii = 1,ncoly2L2*2*nr2L2
         psiy(ii) = 0.0
      enddo

      idum = 1
      do i = 1,nnodes
        do j = 1,ivpacn1lL2(i)
          combuf1(idum) = psi(ivpac1lL2(idum))
          idum = idum + 1
        enddo
      enddo

      call mpi_barrier(mpi_comm_k,ierr)

      idum = 1
      do i = 1,nnodes
       call mpi_isend(combuf1(idum),ivpacn1lL2(i),mpi_real8,i-1,
     &             inode,mpi_comm_k,ireq(i),ierr)
       idum = idum + ivpacn1lL2(i)
      enddo

      idum = 1
      do i = 1,nnodes
       call mpi_recv(combuf2(idum),ivunpn1lL2(i),mpi_real8,i-1,i,
     &               mpi_comm_k,mpistatus,ierr)
       idum = idum + ivunpn1lL2(i)
      enddo

      do i = 1, nnodes
        call mpi_wait(ireq(i), mpistatus, ierr)
      end do

      call mpi_barrier(mpi_comm_k,ierr)

      idum = 1
      do i = 1,nnodes
        do j = 1,ivunpn1lL2(i)
          psiy(ivunp1lL2(idum)) = combuf2(idum)
          idum = idum + 1
        enddo
      enddo
c
c     do FFT on the y direction on the two (one for gamma) slices
c     this is the y dir FFT
c     each PE should have ncoly columns for the FFT
c

      do i = 1,ncoly2L2
         ilocadd = 1+(i-1)*(2*nr2L2)
         call system_ccfft(0,isign,nr2L2,1.0d0,psiy(ilocadd),
     c      psiy(ilocadd),tabnr2linL2,worknr2,0,ntabnr2L2,nworknr2)
      enddo

c
c     now use strided put to put the z strips to the correct
c     processors. Put them into psi.
c     zero values of psi which we do not put into

      do i = 1,ncolz2L2
         idum = (i-1)*2*nr3uL2
         do j = 2*mgz2L2+1,2*nr3uL2
            psi(idum+j)=0.0
         enddo
      enddo

      idum = 1
      do i = 1,nnodes
        do j = 1,ivpacn2lL2(i)
          combuf1(idum) = psiy(ivpac2lL2(idum))
          idum  = idum + 1
        enddo
      enddo

      call mpi_barrier(mpi_comm_k,ierr)

      idum = 1
      do i = 1,nnodes
       call mpi_isend(combuf1(idum),ivpacn2lL2(i),mpi_real8,i-1,inode,
     &                mpi_comm_k,ireq(i),ierr)
       idum = idum + ivpacn2lL2(i)
      enddo

      idum = 1
      do i = 1,nnodes
       call mpi_recv(combuf2(idum),ivunpn2lL2(i),mpi_real8,i-1,i,
     &               mpi_comm_k,mpistatus,ierr)
       idum = idum + ivunpn2lL2(i)
      enddo

      do i = 1, nnodes
        call mpi_wait(ireq(i), mpistatus, ierr)
      end do

      call mpi_barrier(mpi_comm_k,ierr)

      idum = 1
      do i = 1,nnodes
        do j = 1,ivunpn2lL2(i)
          psi(ivunp2lL2(idum)) = combuf2(idum)
          idum = idum + 1
        enddo
      enddo
c

c
c     now do FFT's on z direction so each proc has nr1*nr2/nnodes
c     columns of height nr3 complex numbers (real for gamma point)
c

      do i = 1,ncolzL2
         ilocadd = 1+(i-1)*(2*nr3uL2)

         call system_csfft(0,isign,nr3L2,1.0d0,psi(ilocadd),
     &    psidum,tabnr3lcrL2,worknr3,0,ntabnr3L2,nworknr3)

         iloc_dum = 1+(i-1)*nr3L2

         do ii = 1,nr3L2
            psi(iloc_dum-1+ii)=psidum(ii)
         enddo

      enddo

      call mpi_barrier(mpi_comm_k,ierr)

      deallocate(worknr1)
      deallocate(worknr2)
      deallocate(worknr3)
      deallocate(psiy)
      deallocate(combuf1)
      deallocate(combuf2)

      end

ext: f

max_stars_repo_name: qsnake/petot

max_issues_repo_name: qsnake/petot

max_forks_repo_name: qsnake/petot


text:
      SUBROUTINE XSFADD
C
CXSFABD
C
C     REVISED  8/89 BY G.C./UNISYS
C          1.  THE ORDER OF COMM AND XFIAT IN /XSFA1/ ARE REVERSED IN
C              THIS ROUTINE AND IN THE FOLLOWING 7 SUBROUTINES -
C              XCLEAN, XDPH, XPOLCK, XPUNP, XPURGE, XSFA AND XSOSGN.
C              ANY INCREASE IN SIZE OF XFIAT CAN THEREFORE BE MADE
C              EASILY THROUGH OUT THESE GROUP OF ROUTINES BY JUST
C              CHANGING THE XFIAT DIMENSION HERE.
C          2.  IN THIS GROUP OF ROUTINES, THE ARRAY XFIAT IN /XSFA1/ IS
C              RENAMED TO XFIAT, NOT TO BE CONFUSED WITH THE XFIAT ARRAY
C              IN /XFIAT/
C          3.  ENTN1 MUST EQUAL ICFIAT, THE 24TH WORD OF /SYSTEM/
C              HOWEVER, XSFA AND XPURGE ROUTINES INITIALIZE ENTN1 AGAIN
C              TO ICFIAT, JUST TO BE SURE.
C          4.  THE DIMENSION OF XFIAT SHOULD BE 800 WHEN ENTN1 = 8, OR
C              1100 WHEN ENTN1 IS 11
C
      INTEGER         ALMSK,APNDMK,COMM,CURSNO,ENTN1,ENTN2,ENTN3,
     1                ENTN4,FLAG,FNX,RMSK,RXMSK,S,SCORNT,SOS,TAPMSK,
     2                THCRMK,XFIAT,ZAP
CWKBR COMMON /XSFA1 / MF(401),SOS(1501),COMM(20),XFIAT(1100)
      COMMON /XSFA1 / MF(401),SOS(1501),COMM(20),XFIAT(1320)
      EQUIVALENCE            (COMM (1),ALMSK ),(COMM (2),APNDMK),
     1     (COMM (3),CURSNO),(COMM (4),ENTN1 ),(COMM (5),ENTN2 ),
     2     (COMM (6),ENTN3 ),(COMM (7),ENTN4 ),(COMM (8),FLAG  ),
     3     (COMM (9),FNX   ),(COMM(10),LMSK  ),(COMM(11),LXMSK ),
     4     (COMM(12),MACSFT),(COMM(13),RMSK  ),(COMM(14),RXMSK ),
     5     (COMM(15),S     ),(COMM(16),SCORNT),(COMM(17),TAPMSK),
     6     (COMM(18),THCRMK),(COMM(19),ZAP   )
      ENTN1 = 11
      ENTN2 = 3
      ENTN3 = 4
      ENTN4 = 3
      FLAG  = 0
      DO 10 I = 1, 1320
   10 XFIAT(I) = 0
      TAPMSK = 32768
C            TAPMSK = O 000000100000  = Z 00008000
      APNDMK = 1073741824
C            APNDMK = O 010000000000  = Z 40000000
      RMSK   = 32767
C            RMSK   = O 000000077777  = Z 00007FFF
      RXMSK  = 65535
C            RXMSK  = O 000000177777  = Z 0000FFFF
      LMSK   = 1073676288
C            LMSK   = O 007777600000  = Z 3FFF0000
      LXMSK  = 2147418112
C            LXMSK  = O 017777600000  = Z 7FFF0000
      SCORNT = 1073708992
C            SCORNT = O 007777677700  = Z 3FFF7FC0
      ZAP    = 32767
C            ZAP    = O 000000077777  = Z 00007FFF
      END


text:
!====================================================================
!
! MODULE: ChanType
!
!   Function: Get_ChanType()
!
! Created by Y.Tahara in Nov,02
!
!====================================================================

module ChanType

  !--- modules

  use type_kinds, only : fp_kind
  use ParametersGenCoef

  !--- implicit

  implicit none

  !--- private & public

  private

  public Get_ChanType

 contains

  !====================================================================
  !
  ! Function: Get_ChanType()
  !
  !   Get_ChanType detects type of a particular channel from
  !   its frequency.
  !
  !====================================================================

  function Get_ChanType( wavenum, freq ) result ( chtype )

    !--- interface

    real(fp_kind),intent(in) :: wavenum! wave number (cm^-1)
    real(fp_kind),intent(in) :: freq   ! frequency (GHz)
    integer                  :: chtype ! channel type
                                                ! 1:dry 2:dominated by H2O abs line 3:ozone
                                                ! 4:typically called dry but dominated by H2O continuum

    !--- infrared

    if( wavenum > 100._fp_kind )then

      !--- default is dry

      chtype = 1

      !--- Ozone dominates absorption

      if( wavenum >= 1000._fp_kind .and. wavenum < 1070._fp_kind )then
        chtype = 3
      endif

      !--- H2O absorber lines dominate absorption

      if( wavenum >= 1200._fp_kind .and. wavenum < 2100._fp_kind )then
        chtype = 2
      endif

      !--- H2O continuum affected

      if( wavenum >=  710._fp_kind .and. wavenum < 2100._fp_kind )then
        chtype = chtype + 10
      endif

      return

    endif

    !--- microwave

    !--- default is dry

    chtype = 1

    !--- H2O absorber lines dominate absorption

    if( abs( freq - 23._fp_kind ) < 5._fp_kind )then
      chtype = 2
    endif

    !--- H2O absorbes more than oxigen

    if( (freq >  80._fp_kind .and. freq < 110._fp_kind) .or. &
         freq > 130._fp_kind                            )then
      chtype = 2
    endif

  end function Get_ChanType

end module ChanType

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm


text: >http://www.geocities.com/fish0uttawater/Pictures007.jpg< >Kati was born on Mars. She likes cantaloupes. and Town Wildlife bugs.< >Why cant we all just get along?<

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      module e3if_lhs_m
c
        use e3if_defs_m
c
        implicit none
c
      contains
c
      subroutine set_lhs_matrices
c
        integer :: n,p,q,r
c
        do q = 1,nflow
          do p = 1,nflow
c
            do n = 1,nshl0
c
              AiNa0(:,1,p,q) = Ai0(:,1,p,q) * shp0(:,n)
              AiNa0(:,2,p,q) = Ai0(:,2,p,q) * shp0(:,n)
              AiNa0(:,3,p,q) = Ai0(:,3,p,q) * shp0(:,n)
c
              KijNaj0(:,1,p,q) = Kij0(:,1,1,p,q)*shg0(:,n,1)
     &                         + Kij0(:,1,2,p,q)*shg0(:,n,2)
     &                         + Kij0(:,1,3,p,q)*shg0(:,n,3)
              KijNaj0(:,2,p,q) = Kij0(:,2,1,p,q)*shg0(:,n,1)
     &                         + Kij0(:,2,2,p,q)*shg0(:,n,2)
     &                         + Kij0(:,2,3,p,q)*shg0(:,n,3)
              KijNaj0(:,3,p,q) = Kij0(:,3,1,p,q)*shg0(:,n,1)
     &                         + Kij0(:,3,2,p,q)*shg0(:,n,2)
     &                         + Kij0(:,3,3,p,q)*shg0(:,n,3)
c
            enddo
c
            do n = 1,nshl1
c
              AiNa1(:,1,p,q) = Ai1(:,1,p,q) * shp1(:,n)
              AiNa1(:,2,p,q) = Ai1(:,2,p,q) * shp1(:,n)
              AiNa1(:,3,p,q) = Ai1(:,3,p,q) * shp1(:,n)
c
              KijNaj1(:,1,p,q) = Kij1(:,1,1,p,q)*shg1(:,n,1)
     &                         + Kij1(:,1,2,p,q)*shg1(:,n,2)
     &                         + Kij1(:,1,3,p,q)*shg1(:,n,3)
              KijNaj1(:,2,p,q) = Kij1(:,2,1,p,q)*shg1(:,n,1)
     &                         + Kij1(:,2,2,p,q)*shg1(:,n,2)
     &                         + Kij1(:,2,3,p,q)*shg1(:,n,3)
              KijNaj1(:,3,p,q) = Kij1(:,3,1,p,q)*shg1(:,n,1)
     &                         + Kij1(:,3,2,p,q)*shg1(:,n,2)
     &                         + Kij1(:,3,3,p,q)*shg1(:,n,3)
c
            enddo
c
          enddo
        enddo
c
        KijNajC0 = zero
        KijNajC1 = zero
c
        do q = 1,nflow
          do p = 1,nflow
            do r = 1,nflow
              KijNajC0(:,1,p,q) = KijNajC0(:,1,p,q) + KijNaj0(:,1,p,r)*cmtrx(:,r,q)
              KijNajC0(:,2,p,q) = KijNajC0(:,2,p,q) + KijNaj0(:,2,p,r)*cmtrx(:,r,q)
              KijNajC0(:,3,p,q) = KijNajC0(:,3,p,q) + KijNaj0(:,3,p,r)*cmtrx(:,r,q)
              KijNajC1(:,1,p,q) = KijNajC1(:,1,p,q) + KijNaj1(:,1,p,r)*cmtrx(:,r,q)
              KijNajC1(:,2,p,q) = KijNajC1(:,2,p,q) + KijNaj1(:,2,p,r)*cmtrx(:,r,q)
              KijNajC1(:,3,p,q) = KijNajC1(:,3,p,q) + KijNaj1(:,3,p,r)*cmtrx(:,r,q)
            enddo
          enddo
        enddo
c
      end subroutine set_lhs_matrices
c
      subroutine calc_egmass(egmass,AiNa1,KijNaj0,KijNaj1,KijNajC0,shp0,shp1,n0,n1,WdetJ,nshl0,nshl1)
c
        real*8, dimension(:,:,:), intent(out) :: egmass
        real*8, dimension(:,:,:,:), intent(in) :: AiNa1,KijNaj0,KijNaj1,KijNajC0
        real*8, dimension(:,:), intent(in) :: shp0,shp1,n0,n1
        real*8, dimension(:), intent(in) :: WdetJ
        integer, intent(in) :: nshl0,nshl1
c
        integer :: i,j,p,q,inode0,inode1,isd
c
c... loop through the columns
c
        do inode1 = 1, nshl1
          q = nflow * (inode1 - 1)
c
          do inode0 = 1, nshl0
            p = nflow * (inode0 - 1)
c
            do isd = 1,nsd
              egmass(:,p,q) = egmass(:,p,q) + (
     &              pt50 * shp0(:,inode0) * (AiNa1(:,isd,p,q)+KijNaj1(:,isd,p,q))*n0(:,isd)
c     &            + pt50 * s * KijNajC0(:,isd,p,q)*n1(:,isd)
c     &            + e*mu/h * ctc(:,p,q)*shp0(:,inode0)*n0(:,isd)*shp1(:,inode1)*n1(:,isd)
     &            )*WdetJ
            enddo
c
          enddo
c
        enddo
c
      end subroutine calc_egmass
c
      end module e3if_lhs_m


module stream_module

  use amrex_fort_module, only : amrex_real, dim=>amrex_spacedim

  implicit none

  public

contains

  subroutine pushvtog(lo, hi, dlo, dhi, U, U_lo, U_hi, nc) bind(c,name='pushvtog')
    implicit none
    integer, intent(in) :: nc, lo(3),  hi(3), dlo(3), dhi(3)
    integer, intent(in) :: U_lo(3), U_hi(3)
    real(amrex_real), intent(inout) :: U(U_lo(1):U_hi(1),U_lo(2):U_hi(2),U_lo(3):U_hi(3),nc)

    integer :: n
    real(amrex_real) :: xlo(3), dx(3)

    ! Make up something for these that gets what we want
    dx(1:3)  = 1._amrex_real
    xlo(1:3) = 0._amrex_real

    do n = 1,nc
       call hoextraptocc(U(:,:,:,n),U_lo(1),U_lo(2),U_lo(3),U_hi(1),U_hi(2),U_hi(3),lo,hi,dx,xlo)
    enddo

  end subroutine pushvtog

  subroutine pushvtog2d(lo, hi, dlo, dhi, U, U_lo, U_hi, nc) bind(c,name='pushvtog2d')
    implicit none
    integer, intent(in) :: nc, lo(2),  hi(2), dlo(2), dhi(2)
    integer, intent(in) :: U_lo(2), U_hi(2)
    real(amrex_real), intent(inout) :: U(U_lo(1):U_hi(1),U_lo(2):U_hi(2),nc)

    integer :: n
    real(amrex_real) :: xlo(2), dx(2)

    ! Make up something for these that gets what we want
    dx(1:2)  = 1._amrex_real
    xlo(1:2) = 0._amrex_real

    do n = 1,nc
       call hoextraptocc(U(:,:,n),U_lo(1),U_lo(2),U_hi(1),U_hi(2),lo,hi,dx,xlo)
    enddo

  end subroutine pushvtog2d

end module stream_module


text:
      subroutine zero
*
*
*       initialization of global scalars.
*       ---------------------------------
*
      include 'common.h'
*
      integer i,j
*
*      integer idum2,iv,iy,i
*
c      common /randx/ idum2,iy,iv(32)
*
*      data idum2/123456789/, iv/32*0/, iy/0/
*
*
*       initialize parameters, arrays, counters & set useful constants.
*
      time    = 0.0d0
      timet   = 0.0d0
      ttp     = 0.0d0
      tte     = 0.0d0
      rtid    = 0.0d0
      ekickt  = 0.0d0
      ekicktbs  = 0.0d0
      ekicktbd  = 0.0d0
      ekicktb2  = 0.0d0
      escsta  = 0.0d0
      escbi2  = 0.0d0
      escbi3  = 0.0d0
      escbi4  = 0.0d0
      escbip  = 0.0d0
      escbb2  = 0.0d0
      escbb3  = 0.0d0
      escbb4  = 0.0d0
      escbbp  = 0.0d0
      escb2s  = 0.0d0
      escb3s  = 0.0d0
      escb4s  = 0.0d0
      escbps  = 0.0d0
      ehbin2  = 0.0d0
      ehbin3  = 0.0d0
      ehbi3p  = 0.0d0
      ehb3b3  = 0.0d0
      ehbin4  = 0.0d0
      ehbinp  = 0.0d0
      ecbin2  = 0.0d0
      ehmlev  = 0.0d0
      erelb3  = 0.0d0
      erb3in  = 0.0d0
      enepot  = 0.0d0
      enekin  = 0.0d0
      error   = 0.0d0
      enrad   = 0.0d0
      sloses  = 0.0d0
      slob3s  = 0.0d0
      slob3b  = 0.0d0
      slosev  = 0.0d0
      slosco  = 0.0d0
      sescrt  = 0.0d0
      eccoll  = 0.0d0
      ehcoll  = 0.0d0
      pbin2   = 0.0d0
      pbin3   = 0.0d0
      pbin4   = 0.0d0
      pbin2s  = 0.0d0
      pbin3s  = 0.0d0
      pbin4s  = 0.0d0
      pbinps  = 0.0d0
      pb2b2   = 0.0d0
      pb2b3   = 0.0d0
      pb2b4   = 0.0d0
      pb2bp   = 0.0d0
      pb3b3   = 0.0d0
      pb3b4   = 0.0d0
      pb3bp   = 0.0d0
      pb4b4   = 0.0d0
      pb4bp   = 0.0d0
      pbpbp   = 0.0d0
      pcoll   = 0.0d0
*
      mnsbh   = 0
      ikickt  = 0
      ikicktbs  = 0
      ikicktbd  = 0
      nexchang = 0
      ntsn1 = 0
      ntsn2 = 0
      ntsnb = 0
      nmerge  = 0
      nescst  = 0
      nescb2  = 0
      nescb3  = 0
      nescb4  = 0
      nescrt  = 0
      nescbp  = 0
      nesb2s  = 0
      nesb3s  = 0
      nesb4s  = 0
      nesbps  = 0
      nmloev  = 0
      nmloco  = 0
      iobt = 0
      ncoll   = 0
      necb2   = 0
      nehb2   = 0
      necb3   = 0
      nehb3   = 0
      necb4   = 0
      nehb4   = 0
      necbp   = 0
      nehbp   = 0
      nbin2   = 0
      nbin3   = 0
      nb3b3   = 0
      nb3fin  = 0
      nbin4   = 0
      nbinp   = 0
      ndist3  = 0
      ndist4  = 0
      iesc1   = 0
      ivnewg  = 0
      ivrr    = 0
      ibsm    = 0
      ibsc    = 0
      ibs3    = 0
      ibs4    = 0
*
      do 20 i=1,32
         iv(i) = 0
 20   continue
*
      idum2 = 123456789
      iy = 0
*
*       set fractional constants pi & two pi
*
      one2  = 1.0d0/2.0d0
      one3  = 1.0d0/3.0d0
      one4  = 1.0d0/4.0d0
      one5  = 1.0d0/5.0d0
      one6  = 1.0d0/6.0d0
      one9  = 1.0d0/9.0d0
      one12 = 1.0d0/12.0d0
      pi = 4.0d0*atan(1.0d0)
      twopi = 2.0d0*pi
      rmglob = 1.8d0/float(nt)
      rsuntopc = 6.9599e10/3.085678e18
*
      do 30 i = 1,nmax
         nkick(i) = 0
         inexch(i) = 0
         iname(i) = 0
         ikind(i) = 0
         names(i) = 0
         nameb(i) = 0
         nwhich(i) = 0
         nbinar(i) = 0
         ibstra(i) = 0
         body(i)= 0.d0
         vr(i) = 0.d0
         vt(i) = 0.d0
         r(i) = 0.d0
         u(i) = 0.d0
         xescp(i) = 0.d0
         xesct(i) = 0.d0
         vrr(i) = 0.d0
         ro(i) = 0.d0
         vro(i) = 0.d0
         vto(i) = 0.d0
         uo(i) = 0.d0
         rn(i) = 0.d0
         xmin(i) = 0.d0
         xmax(i) = 0.d0
         xmaxz(i) = 0.d0
         xgmin(i) = 0.d0
         xgmax(i) = 0.d0
         x(i,1) = 0.d0
         x(i,2) = 0.d0
         x(i,3) = 0.d0
         xdot(i,1) = 0.d0
         xdot(i,2) = 0.d0
         xdot(i,3) = 0.d0
         uptime(i) = 0.d0
         oldtime(i) = 0.d0
         vrp(i) = 0.d0
         vtp(i) = 0.d0
 30   continue
*
      do 40 i = 1,nbmax3
         iinte3(i) = 0
         do 41 j = 1,8
             bin(i,j) = 0.d0
 41      continue
 40   continue
*
      do 50 i =1,50*nbmax3
         do 51 j = 1,7
            binin(i,j) = 0.d0
 51      continue
 50   continue
*
      return
*
      end
*
*
*
*

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
! *****************************COPYRIGHT*******************************
! (C) Crown copyright Met Office. All rights reserved.
! For further details please refer to the file COPYRIGHT.txt
! which you should have received as part of this distribution.
! *****************************COPYRIGHT*******************************
! Contains values of all physical constants as parameters
!      CP, R, RATIO_MOL_WTS, G, RLVAP, RLSUB,
! and functions of the above
!      RLVAP_ON_CP, RLSUB_ON_CP, R_ON_CP, G_ON_2
Module physconst

  Use typeKind

  Implicit None

  REAL(WP), PARAMETER :: CP = 1005.
  ! Specific heat of gas at constant pressure
  REAL(WP), PARAMETER :: R = 287.05
  ! Universal gas constant
  REAL(WP), PARAMETER :: RATIO_MOL_WTS = 1.608
  ! Ratio of molecular weights of air and water vapour
  REAL(WP), PARAMETER :: G = 9.81
  ! Acceleration due to gravity
  REAL(WP), PARAMETER :: RLVAP = 2.501e6
  ! Latent heat of vapourisation
  REAL(WP), PARAMETER :: RLSUB = 2.834e6
  ! Latent heat of sublimation
  REAL(WP), PARAMETER :: RLFUS = RLSUB - RLVAP
  ! Latent heat of fusion
  REAL(WP), PARAMETER :: RLVAP_ON_CP = RLVAP/CP
  REAL(WP), PARAMETER :: RLSUB_ON_CP = RLSUB/CP
  REAL(WP), PARAMETER :: RLFUS_ON_CP = RLSUB_ON_CP-RLVAP_ON_CP
  REAL(WP), PARAMETER :: R_ON_CP = R/CP
  REAL(WP), PARAMETER :: G_ON_2 = 0.5*G
  ! Time data which is useful to keep constant
  REAL(WP), PARAMETER ::      &
  &   secinhr   = 3600.   &   ! second in hour
  & , hrinday   = 24.     &   ! hours in day
  & , secsinday = 86400.  &   ! seconds in day
  & , rsecsinday = 1./secsinday

  real(wp), parameter :: p0=100000. ! reference surface pressure
  real(wp), parameter :: pi=3.141592653589793239_wp ! pi

  real(wp), parameter :: tk0c=273.15 ! temperature in Kelvin at
                                     ! 0 degrees C

  real(wp), parameter :: rhoW   = 1000.   ! density of water
  real(wp), parameter :: Rw     = 461.5   ! gas conbstant for water vapour
  real(wp), parameter :: THcond = 0.0243  ! thermal conductivity of air
  real(wp), parameter :: diffWV = 2.26E-5 ! diffusivity of water vapour in air
  real(wp), parameter :: visair = 1.44E-5 ! kinematic viscosity of air
  real(wp), parameter :: Cwater = 4187.   ! specific heat capacity of liquid water
  real(wp), parameter :: Cice   = 2093.   ! specific heat capacity of ice

end Module physconst

max_stars_repo_name: Adehill/KiD-A

max_issues_repo_name: Adehill/KiD-A

max_forks_repo_name: Adehill/KiD-A


 FUNCTION JONSWAP_spectrum(FREQ, Hs, Tp, gamma)
!
! This function returns the value of the JONSWAP spectrum
! for the frequency FREQ, significant wave height Hs, and peak period
! Tp.

   IMPLICIT none
   integer, parameter :: long=selected_real_kind(12,99)
   REAL(kind=long) :: jonswap_spectrum, Hs, Tp, A, B, a_exp, alpha, gamma, &
        sigma, freq, f_p, four=4._long, pi, zero=0._long, one=1._long,     &
        two=2._long, half=.5_long, five=5._long, grav=9.82_long
   pi = acos (-one)
   f_p=one/Tp
   alpha=3.32285_long * Hs**2 * f_p**4
   A=alpha*grav**2/(two*pi)**4
   B=five/four * f_p**4

   IF (freq.le.zero) then
      jonswap_spectrum = zero
   ELSEIF(freq<=f_p)THEN
      sigma=.07_long
      a_exp=exp(-half*((freq-f_p)/(sigma*f_p))**2)
      jonswap_spectrum = A/freq**5 * exp(-B/freq**4) * gamma**a_exp
   ELSE
      sigma=.09_long
      a_exp=exp(-half*((freq-f_p)/(sigma*f_p))**2)
      jonswap_spectrum = A/freq**5 * exp(-B/freq**4) * gamma**a_exp
   ENDIF
   RETURN
 END FUNCTION JONSWAP_spectrum


text:
module strength
use hamiltonian_types
use maths
use build_hamiltonian_routines
implicit none
contains

!------------------------------------------------
! Locates a given value in a vector
!------------------------------------------------
	function locate(a,b)
	integer :: locate
	real(kind=8), INTENT(IN) :: a(:), b
	integer :: location(1)
		locate = -1

		location = minloc(abs(a-b))
		if (location(1) /= 0) then
			locate = location(1)
		endif

	end function locate

!------------------------------------------------
! Returns the kronecker's delta
!------------------------------------------------
	function kronecker(a,b)
	real*8 :: kronecker
	real*8, INTENT(IN) :: a, b
		kronecker = 0.d0
		if (a == b) then
			kronecker = 1.d0
		endif
	end function kronecker

! ------------------------------------------------------------
! This function returns the matrix elements of the dipole moment in the Hund's a)
! eigenfunctions
! ------------------------------------------------------------
	function q_o_m(Lambda_up,Omega_up,J_up,M_up,Lambda_low,Omega_low,J_low,M_low)
	real*8 :: q_o_m
	real*8, INTENT(IN) :: Lambda_up, Omega_up, J_up, M_up, Lambda_low, Omega_low, J_low, M_low
	real*8 :: q_m, q_o

		q_m = qm(J_up,M_up,J_low,M_low)
		q_o = qo(Lambda_up,Omega_up,J_up,Lambda_low,Omega_low,J_low)

		q_o_m = q_o * q_m

	end function q_o_m

! ------------------------------------------------------------
! This function returns the matrix elements of q_o
! ------------------------------------------------------------
	function qm(J_up,M_up,J_low,M_low)
	real*8 :: qm
	real*8, INTENT(IN) :: J_up, M_up, J_low, M_low
	real*8 :: deltam, deltaj, J, M

		deltaj = J_up - J_low
		deltam = M_up - M_low

		J = J_low
		M = M_low
		qm = 0.d0

		if (deltam == 0.d0) then
			if (deltaj == 1.d0) then
				qm = dsqrt((J+1.d0)**2 - M**2)
			else if (deltaj == 0) then
				qm = M
			else if (deltaj == -1.d0) then
				qm = dsqrt(J**2 - M**2)
			endif
		endif

		if (dabs(deltam) == 1.d0) then
			if (deltaj == 1.d0) then
				qm = -deltam * dsqrt(0.5d0*(J+1.d0+M*deltam)*(J+2.d0+M*deltam))
			else if (deltaj == 0) then
				qm = dsqrt(0.5d0*(J-M*deltam)*(J+1.d0+M*deltam))
			else if (deltaj == -1.d0) then
				qm = deltam * dsqrt(0.5d0*(J-M*deltam)*(J-1.d0-M*deltam))
			endif
		endif

	end function qm

! ------------------------------------------------------------
! This function returns the matrix elements of q_o
! ------------------------------------------------------------
	function qo(Lambda_up,Omega_up,J_up,Lambda_low,Omega_low,J_low)
	real*8 :: qo
	real*8, INTENT(IN) :: Lambda_up, Omega_up, J_up, Lambda_low, Omega_low, J_low
	real*8 :: deltao, deltaj, Omega, J, pi_eps

		deltao = Omega_up - Omega_low
		deltaj = J_up - J_low
		pi_eps = 1.d0 + kronecker(1.d0,Lambda_up+Lambda_low)

		J = J_low
		Omega = Omega_low
		qo = 0.d0

		if (deltao == 0.d0) then
			if (deltaj == 1.d0) then
				qo = 1.d0/(J+1.d0) * dsqrt(( (J+1.d0)**2 - Omega**2)/( (2.d0*J+1.d0)*(2.d0*J+3.d0)) )
			else if (deltaj == 0) then
				qo = Omega / (J*(J+1.d0))
			else if (deltaj == -1.d0) then
				qo = 1.d0/J * dsqrt(( J**2 - Omega**2)/((2.d0*J-1.d0)*(2.d0*J+1.d0)))
			endif
		endif

		if (dabs(deltao) == 1.d0) then
			if (deltaj == 1.d0) then
				qo = -deltao/(J+1.d0) * dsqrt(pi_eps*(J+1.d0+Omega*deltao)*(J+2.d0+Omega*deltao) /&
					(2.d0*(2.d0*J+1.d0)*(2.d0*J+3.d0)))
			else if (deltaj == 0) then
				qo = 1.d0/(J*(J+1.d0)) * dsqrt(0.5d0*pi_eps*(J-Omega*deltao)*(J+1.d0+Omega*deltao))
			else if (deltaj == -1.d0) then
				qo = deltao/J * dsqrt(pi_eps*(J-Omega*deltao)*(J-1.d0-Omega*deltao) /&
					(2.d0*(2.d0*J-1.d0)*(2.d0*J+1.d0)))
			endif
		endif

	end function qo

!----------------------------------------------------------------------
!----------------------------------------------------------------------
! Strength of the Zeeman components in the Hund's case (a) coupling
!----------------------------------------------------------------------
!----------------------------------------------------------------------
	function case_a_strength(Ju,Jl,Mu,Ml,Ou,Ol)
	real(kind=8) :: case_a_strength
	real(kind=8) :: Mu, Ml, Ju, Jl, Ou, Ol, qm, qo

		qm = Mu - Ml
		qo = Ou - Ol

! Falta un factor sqrt((2Ju+1)*(2Jl+1)) ??

		case_a_strength = w3js(int(2.d0*Jl),int(2.d0*Ju),2,int(2.d0*Ml),-int(2.d0*Mu),int(2.d0*qm)) * &
			w3js(int(2.d0*Jl),int(2.d0*Ju),2,int(2.d0*Ol),-int(2.d0*Ou),int(2.d0*qo)) * &
			sqrt((2.d0*Ju+1.d0)*(2.d0*Jl+1.d0))

	end function case_a_strength

!----------------------------------------------------------------------
!----------------------------------------------------------------------
! Count the number of Zeeman components and generate the index arrays
!----------------------------------------------------------------------
!----------------------------------------------------------------------
	subroutine count_zeeman_components(trans)
	type(type_transition) :: trans
	integer :: i, j, i2, j2, k, nJup, nJlow, ind_up, ind_low
	real(kind=8) :: Ju, Jl, Mu, Ml, Ou, Ol, case_a_str, spli_up, spli_low, spli, m_up, m_low
	real(kind=8), allocatable :: direction_cosine(:,:), strength_matrix(:,:), caca(:,:)

		nJup = int(2.d0*trans%upper%J+1.d0)
		nJlow = int(2.d0*trans%lower%J+1.d0)

		k = 0
		do i = 1, trans%upper%hamiltonian_size
			do j = 1, trans%lower%hamiltonian_size
				if ( abs(trans%upper%ket(i)%M-trans%lower%ket(j)%M) <= 1.d0 .and. &
					trans%upper%active(i) /= -9999 .and. trans%lower%active(j) /= -9999) then
						k = k + 1
				endif
			enddo
		enddo
		print *, 'Number of Zeeman transitions : ', k

		allocate(direction_cosine(trans%upper%hamiltonian_size,trans%lower%hamiltonian_size))
		allocate(strength_matrix(trans%upper%hamiltonian_size,trans%lower%hamiltonian_size))

		do i = 1, trans%upper%hamiltonian_size
			do j = 1, trans%lower%hamiltonian_size
				Ju = trans%upper%ket_m_ordered(i)%J
				Jl = trans%lower%ket_m_ordered(j)%J
				Mu = trans%upper%ket_m_ordered(i)%M
				Ml = trans%lower%ket_m_ordered(j)%M
				Ou = trans%upper%ket_m_ordered(i)%Omega
				Ol = trans%lower%ket_m_ordered(j)%Omega

				direction_cosine(i,j) = kronecker(trans%upper%ket_m_ordered(i)%Sigma,trans%lower%ket_m_ordered(j)%Sigma) * &
					q_o_m(trans%upper%Lambda,Ou,Ju,Mu,trans%lower%Lambda,Ol,Jl,Ml)
			enddo
		enddo

		strength_matrix = matmul(matmul(transpose(trans%upper%eigenvec),direction_cosine),trans%lower%eigenvec)

		allocate(trans%ind_up(k))
		allocate(trans%ind_low(k))
		allocate(trans%str(k))
		trans%str = 0.d0

		print *, trans%upper%hamiltonian_size, trans%lower%hamiltonian_size
		write(18,*) strength_matrix
		write(17,*) direction_cosine

		open(unit=12,file='pattern.dat',action='write',status='replace')
		write(12,*) k

		k = 0
		do i = 1, 2*trans%upper%J+1
			m_up = -trans%upper%J + (i-1)
			do j = 1, 2*trans%lower%J+1
				m_low = -trans%lower%J + (j-1)

				if (dabs(m_up-m_low) <= 1.d0) then

					ind_up = locate(trans%upper%active,m_up)
					ind_low = locate(trans%lower%active,m_low)

					spli_up = trans%upper%eigenval(ind_up)-trans%upper%eigenval_zerofield(ind_up)
					spli_low = trans%lower%eigenval(ind_low)-trans%lower%eigenval_zerofield(ind_low)
					spli = -(spli_up-spli_low)*(trans%wavelength*1.d-8)**2 * 1.d8
					write(*,FMT='(F5.1,1X,F5.1,1X,F13.8,3X,F13.8,3X,E13.7)') trans%upper%ket(ind_up)%M, &
							trans%lower%ket(ind_low)%M, spli, strength_matrix(ind_up,ind_low)**2
					write(12,FMT='(F5.1,1X,F5.1,1X,F13.8,3X,F13.8,3X,E13.7)') trans%upper%ket(ind_up)%M, &
							trans%lower%ket(ind_low)%M, spli, strength_matrix(ind_up,ind_low)**2

				endif
			enddo
		enddo
		close(12)

	end subroutine count_zeeman_components

!----------------------------------------------------------------------
!----------------------------------------------------------------------
! Strength of the Zeeman components
!----------------------------------------------------------------------
!----------------------------------------------------------------------
	subroutine zeeman_strength(trans)
	type(type_transition) :: trans

! First count the number of components of the Zeeman transition
		call count_zeeman_components(trans)

	end subroutine zeeman_strength

end module strength

max_stars_repo_name: aasensio/diatomic

max_issues_repo_name: aasensio/diatomic

max_forks_repo_name: aasensio/diatomic


text:
module mod_band_types
  use mod_prec
  use mod_rotation
  use mod_utility
  implicit none
  ! Routines for working with compact representations of banded
  ! matrices.  A banded matrix $A=[a_{jk}]$ with $a_{jk}=0$ for
  ! $k-j>s$ and $k-j < -r$ is stored in one of two ways: In the first
  ! format, elements of $A$ are stored in a matrix $B$ with columns
  ! aligned with those of $A$ and $a_{j,k}=b_{j-k+s+1,k}$. Thus each
  ! row of $B$ corresponds to a diagonal of $A$ and $B$ can be chosen
  ! to be $(r+s+1)\times n$.  Routines working with this format are
  ! given the tag `_bc`. In the second format, $B$ has rows aligned
  ! with those of $A$ and $a_{j,k}=b_{j,k-j+r+1}$.  In this case, each
  ! column of $B$ corresponds to a diagonal of $A$ and $B$ can be
  ! chosen to be $n\times (r+s+1)$.  Routines working with this format
  ! are given the tag `_br`.  The module provides routines for getting
  ! and setting elements $a_{jk}$ of a matrix $A$ stored in this
  ! format; for applying applying rotations in a truncated manner
  ! corresponding to $j$-leading or $j$-trailing transformations; for
  ! printing band matrices in a readable format; and for converting
  ! between different storage formats.

  private

  public :: get_el_bc, d_get_el_bc, z_get_el_bc, &
       set_el_bc, d_set_el_bc, z_set_el_bc, &
       get_el_br, d_get_el_br, z_get_el_br, &
       set_el_br, d_set_el_br, z_set_el_br

  public :: bc_to_general, d_bc_to_general, z_bc_to_general, &
       br_to_general, d_br_to_general, z_br_to_general

  public :: print_bc, d_print_bc, z_print_bc, &
       print_br, d_print_br, z_print_br, &
       print_abs_bc, d_print_abs_bc, z_print_abs_bc, &
       print_abs_br, d_print_abs_br, z_print_abs_br

  public :: f_rotation_times_tbc, f_d_rotation_times_tbc, &
       f_z_rotation_times_tbc, f_tbc_times_rotation, f_d_tbc_times_rotation, &
       f_z_tbc_times_rotation, f_tbr_times_rotation, f_d_tbr_times_rotation, &
       f_z_tbr_times_rotation, f_rotation_times_tbr, f_d_rotation_times_tbr, &
       f_z_rotation_times_tbr

  public :: rotation_times_tbc, d_rotation_times_tbc, z_rotation_times_tbc, &
       tbc_times_rotation, d_tbc_times_rotation, z_tbc_times_rotation, &
       tbr_times_rotation, d_tbr_times_rotation, z_tbr_times_rotation, &
       rotation_times_tbr, d_rotation_times_tbr, z_rotation_times_tbr

  public :: bc_to_br, d_bc_to_br, z_bc_to_br, &
       br_to_bc, d_br_to_bc, z_br_to_bc

  public :: extract_diagonals_bc, d_extract_diagonals_bc, &
       z_extract_diagonals_bc, extract_diagonals_br, d_extract_diagonals_br, &
       z_extract_diagonals_br

  public :: truncate_profile_bc, d_truncate_profile_bc, z_truncate_profile_bc, &
       truncate_profile_br, d_truncate_profile_br, z_truncate_profile_br

  public :: d_band_norm_bc, z_band_norm_bc, band_norm_bc, &
       d_band_norm_br, z_band_norm_br, band_norm_br

  interface get_el_bc
     module procedure d_get_el_bc, z_get_el_bc
  end interface get_el_bc

  interface set_el_bc
     module procedure d_set_el_bc, z_set_el_bc
  end interface set_el_bc

  interface get_el_br
     module procedure d_get_el_br, z_get_el_br
  end interface get_el_br

  interface set_el_br
     module procedure d_set_el_br, z_set_el_br
  end interface set_el_br

  ! to general

  interface bc_to_general
     module procedure d_bc_to_general, z_bc_to_general
  end interface bc_to_general

  interface br_to_general
     module procedure d_br_to_general, z_br_to_general
  end interface br_to_general

  ! printing

  interface print_bc
     module procedure d_print_bc, z_print_bc
  end interface print_bc

  interface print_br
     module procedure d_print_br, z_print_br
  end interface print_br

  interface print_abs_bc
     module procedure d_print_abs_bc, z_print_abs_bc
  end interface print_abs_bc

  interface print_abs_br
     module procedure d_print_abs_br, z_print_abs_br
  end interface print_abs_br

  ! rotations

  interface f_rotation_times_tbc
     module procedure f_d_rotation_times_tbc, f_z_rotation_times_tbc
  end interface f_rotation_times_tbc

  interface f_tbc_times_rotation
     module procedure f_d_tbc_times_rotation, f_z_tbc_times_rotation
  end interface f_tbc_times_rotation

  interface f_tbr_times_rotation
     module procedure f_d_tbr_times_rotation, f_z_tbr_times_rotation
  end interface f_tbr_times_rotation

  interface f_rotation_times_tbr
     module procedure f_d_rotation_times_tbr, f_z_rotation_times_tbr
  end interface f_rotation_times_tbr

  interface rotation_times_tbc
     module procedure d_rotation_times_tbc, z_rotation_times_tbc
  end interface rotation_times_tbc

  interface tbc_times_rotation
     module procedure d_tbc_times_rotation, z_tbc_times_rotation
  end interface tbc_times_rotation

  interface tbr_times_rotation
     module procedure d_tbr_times_rotation, z_tbr_times_rotation
  end interface tbr_times_rotation

  interface rotation_times_tbr
     module procedure d_rotation_times_tbr, z_rotation_times_tbr
  end interface rotation_times_tbr

  ! conversions

  interface bc_to_br
     module procedure d_bc_to_br, z_bc_to_br
  end interface bc_to_br

  interface br_to_bc
     module procedure d_br_to_bc, z_br_to_bc
  end interface br_to_bc

  interface extract_diagonals_br
     module procedure d_extract_diagonals_br, z_extract_diagonals_br
  end interface extract_diagonals_br

  interface extract_diagonals_bc
     module procedure d_extract_diagonals_bc, z_extract_diagonals_bc
  end interface extract_diagonals_bc

  interface truncate_profile_bc
     module procedure d_truncate_profile_bc, z_truncate_profile_bc
  end interface truncate_profile_bc

  interface truncate_profile_br
     module procedure d_truncate_profile_br, z_truncate_profile_br
  end interface truncate_profile_br

  interface band_norm_bc
     module procedure d_band_norm_bc, z_band_norm_bc
  end interface band_norm_bc

  interface band_norm_br
     module procedure d_band_norm_br, z_band_norm_br
  end interface band_norm_br

contains

  ! Get and set functions for band row and band column.

  real(kind=dp) function d_get_el_bc(bc,ubw,j,k) result(x)
    real(kind=dp), dimension(:,:), intent(in) :: bc
    integer(kind=int32), intent(in) :: ubw,j,k
    !
    x=bc(j-k+ubw+1,k)
  end function d_get_el_bc

  complex(kind=dp) function z_get_el_bc(bc,ubw,j,k) result(x)
    complex(kind=dp), dimension(:,:), intent(in) :: bc
    integer(kind=int32), intent(in) :: ubw,j,k
    !
    x=bc(j-k+ubw+1,k)
  end function z_get_el_bc

  subroutine d_set_el_bc(bc,ubw,j,k,x)
    real(kind=dp), dimension(:,:), intent(inout) :: bc
    integer(kind=int32), intent(in) :: ubw,j,k
    real(kind=dp), intent(in) :: x
    !
    bc(j-k+ubw+1,k) = x
  end subroutine d_set_el_bc

  subroutine z_set_el_bc(bc,ubw,j,k,x)
    complex(kind=dp), dimension(:,:), intent(inout) :: bc
    integer(kind=int32), intent(in) :: ubw,j,k
    complex(kind=dp), intent(in) :: x
    !
    bc(j-k+ubw+1,k) = x
  end subroutine z_set_el_bc

  real(kind=dp) function d_get_el_br(br,lbw,j,k) result(x)
    real(kind=dp), dimension(:,:), intent(in) :: br
    integer(kind=int32), intent(in) :: lbw,j,k
    !
    x=br(j,k-j+lbw+1)
  end function d_get_el_br

  complex(kind=dp) function z_get_el_br(br,lbw,j,k) result(x)
    complex(kind=dp), dimension(:,:), intent(in) :: br
    integer(kind=int32), intent(in) :: lbw,j,k
    !
    x=br(j,k-j+lbw+1)
  end function z_get_el_br

  subroutine d_set_el_br(br,lbw,j,k,x)
    real(kind=dp), dimension(:,:), intent(inout) :: br
    integer(kind=int32), intent(in) :: lbw,j,k
    real(kind=dp), intent(in) :: x
    !
    br(j,k-j+lbw+1)=x
  end subroutine d_set_el_br

  subroutine z_set_el_br(br,lbw,j,k,x)
    complex(kind=dp), dimension(:,:), intent(inout) :: br
    integer(kind=int32), intent(in) :: lbw,j,k
    complex(kind=dp), intent(in) :: x
    !
    br(j,k-j+lbw+1)=x
  end subroutine z_set_el_br

  ! band to general

  subroutine d_bc_to_general(bc,lbw,ubw,a)
    real(kind=dp), dimension(:,:), intent(in) :: bc
    integer(kind=int32), intent(in) :: ubw, lbw
    real(kind=dp), dimension(:,:), intent(out) :: a
    !
    integer(kind=int32) :: d, n, k
    a=0.0_dp
    n=size(bc,2)
    do d=1,ubw+1
       do k=ubw-d+2,n
          a(d+k-ubw-1,k)=bc(d,k)
       end do
    end do
    do d=ubw+2, ubw+lbw+1
       do k=1,n-d+ubw+1
          a(d+k-ubw-1,k)=bc(d,k)
       end do
    end do
  end subroutine d_bc_to_general

  subroutine z_bc_to_general(bc,lbw,ubw,a)
    complex(kind=dp), dimension(:,:), intent(in) :: bc
    integer(kind=int32), intent(in) :: ubw, lbw
    complex(kind=dp), dimension(:,:), intent(out) :: a
    !
    integer(kind=int32) :: d, n, k
    a=(0.0_dp,0.0_dp)
    n=size(bc,2)
    do d=1,ubw+1
       do k=ubw-d+2,n
          a(d+k-ubw-1,k)=bc(d,k)
       end do
    end do
    do d=ubw+2, ubw+lbw+1
       do k=1,n-d+ubw+1
          a(d+k-ubw-1,k)=bc(d,k)
       end do
    end do
  end subroutine z_bc_to_general

  subroutine d_br_to_general(br,lbw,ubw,a)
    real(kind=dp), dimension(:,:), intent(in) :: br
    integer(kind=int32), intent(in) :: ubw, lbw
    real(kind=dp), dimension(:,:), intent(out) :: a
    !
    integer(kind=int32) :: d, n, j
    a=0.0_dp
    n=size(br,1)
    do d=1,lbw+1
       do j=lbw-d+2, n
          a(j,d+j-lbw-1)=br(j,d)
       end do
    end do
    do d=lbw+2,ubw+lbw+1
       do j=1, n-d+lbw+1
          a(j,d+j-lbw-1)=br(j,d)
       end do
    end do
  end subroutine d_br_to_general

  subroutine z_br_to_general(br,lbw,ubw,a)
    complex(kind=dp), dimension(:,:), intent(in) :: br
    integer(kind=int32), intent(in) :: ubw, lbw
    complex(kind=dp), dimension(:,:), intent(out) :: a
    !
    integer(kind=int32) :: d, n, j
    a=(0.0_dp,0.0_dp)
    n=size(br,1)
    do d=1,lbw+1
       do j=lbw-d+2, n
          a(j,d+j-lbw-1)=br(j,d)
       end do
    end do
    do d=lbw+2,ubw+lbw+1
       do j=1, n-d+lbw+1
          a(j,d+j-lbw-1)=br(j,d)
       end do
    end do
  end subroutine z_br_to_general

  ! printing

  subroutine d_print_bc(bc,lbw,ubw)
    real(kind=dp), dimension(:,:), intent(in) :: bc
    integer(kind=int32), intent(in) :: ubw, lbw
    !
    real(kind=dp), dimension(size(bc,2),size(bc,2)) :: a
    call d_bc_to_general(bc,lbw,ubw,a)
    call print_matrix(a)
  end subroutine d_print_bc

  subroutine z_print_bc(bc,lbw,ubw)
    complex(kind=dp), dimension(:,:), intent(in) :: bc
    integer(kind=int32), intent(in) :: ubw, lbw
    !
    complex(kind=dp), dimension(size(bc,2),size(bc,2)) :: a
    call z_bc_to_general(bc,lbw,ubw,a)
    call print_matrix(a)
  end subroutine z_print_bc

  subroutine d_print_br(br,lbw,ubw)
    real(kind=dp), dimension(:,:), intent(in) :: br
    integer(kind=int32), intent(in) :: ubw, lbw
    !
    real(kind=dp), dimension(size(br,1),size(br,1)) :: a
    call d_br_to_general(br,lbw,ubw,a)
    call print_matrix(a)
  end subroutine d_print_br

  subroutine z_print_br(br,lbw,ubw)
    complex(kind=dp), dimension(:,:), intent(in) :: br
    integer(kind=int32), intent(in) :: ubw, lbw
    !
    complex(kind=dp), dimension(size(br,1),size(br,1)) :: a
    call z_br_to_general(br,lbw,ubw,a)
    call print_matrix(a)
  end subroutine z_print_br

  ! abs print

  subroutine d_print_abs_bc(bc,lbw,ubw)
    real(kind=dp), dimension(:,:), intent(in) :: bc
    integer(kind=int32), intent(in) :: ubw, lbw
    !
    real(kind=dp), dimension(size(bc,2),size(bc,2)) :: a
    call d_bc_to_general(bc,lbw,ubw,a)
    call print_matrix(abs(a))
  end subroutine d_print_abs_bc

  subroutine z_print_abs_bc(bc,lbw,ubw)
    complex(kind=dp), dimension(:,:), intent(in) :: bc
    integer(kind=int32), intent(in) :: ubw, lbw
    !
    complex(kind=dp), dimension(size(bc,2),size(bc,2)) :: a
    call z_bc_to_general(bc,lbw,ubw,a)
    call print_matrix(abs(a))
  end subroutine z_print_abs_bc

  subroutine d_print_abs_br(br,lbw,ubw)
    real(kind=dp), dimension(:,:), intent(in) :: br
    integer(kind=int32), intent(in) :: ubw, lbw
    !
    real(kind=dp), dimension(size(br,1),size(br,1)) :: a
    call d_br_to_general(br,lbw,ubw,a)
    call print_matrix(abs(a))
  end subroutine d_print_abs_br

  subroutine z_print_abs_br(br,lbw,ubw)
    complex(kind=dp), dimension(:,:), intent(in) :: br
    integer(kind=int32), intent(in) :: ubw, lbw
    !
    complex(kind=dp), dimension(size(br,1),size(br,1)) :: a
    call z_br_to_general(br,lbw,ubw,a)
    call print_matrix(abs(a))
  end subroutine z_print_abs_br

  ! Rotations applied to band matrices.
  !
  ! In applying rotations to a band matrix:
  !
  ! Column k has nonzeros in rows [max(k-ubw,1),min(k+lbw,m)]
  ! Column k+1 has nonzeros in rows [max(k-ubw+1,1),min(k+1+lbw,m)]
  ! A rotation acting from the right in columns k and k+1
  ! acts on the intersection of these row sets: [max(k-ubw+1,1),min(k+1+lbw,m)].
  !
  ! Row j has nonzeros in columns [max(1,j-lbw), min(j+ubw,n)]
  ! Row j+1 has nonzeros in columns [max(1,j+1-lbw), min(j+1+ubw,n)]
  ! A rotation acting on rows j and j+1 acts on the intersection:
  ! [max(1,j+1-lbw),min(n,j+ubw)]

  ! Rotations for band matrices with aligned columns
  !
  ! Rotation times truncated band matrix with aligned columns.
  ! First l and last p columns are not modified.

  subroutine f_d_rotation_times_tbc(c,s,b,n,lbw,ubw,l,p,j)
    real(kind=dp), dimension(:,:), intent(inout) :: b
    real(kind=dp), intent(in) :: c, s
    integer(kind=int32), intent(in) :: j,n, ubw,lbw, l, p
    !
    real(kind=dp) :: tmp
    integer(kind=int32) :: k, k0, k1, d
    ! first and last columns operated on, given by the
    ! intersection of [l+1,n-p] and [max(1,j+1-lbw),min(n,j+ubw)]
    ! which is [maxl+1,max(1,j+1-lbw)), min(n-p,min(n,j+ubw))]
    k0=max(l+1,j+1-lbw)
    k1=min(n-p,j+ubw)
    ! loop over relevant columns in b.
    do k=k0,k1
       d=j-k+ubw+1
       tmp = b(d,k)
       b(d,k)=c*tmp - s*b(d+1,k)
       b(d+1,k) = s*tmp + c*b(d+1,k)
    end do
  end subroutine f_d_rotation_times_tbc

  subroutine d_rotation_times_tbc(r,b,n,lbw,ubw,l,p,j)
    real(kind=dp), dimension(:,:), intent(inout) :: b
    type(d_rotation), intent(in) :: r
    integer(kind=int32), intent(in) :: j,n, ubw,lbw, l, p
    !
    real(kind=dp) :: c, s, tmp
    integer(kind=int32) :: k, k0, k1, d
    c=r%cosine; s=r%sine
    ! first and last columns operated on, given by the
    ! intersection of [l+1,n-p] and [max(1,j+1-lbw),min(n,j+ubw)]
    ! which is [maxl+1,max(1,j+1-lbw)), min(n-p,min(n,j+ubw))]
    k0=max(l+1,j+1-lbw)
    k1=min(n-p,j+ubw)
    ! loop over relevant columns in b.
    do k=k0,k1
       d=j-k+ubw+1
       tmp = b(d,k)
       b(d,k)=c*tmp - s*b(d+1,k)
       b(d+1,k) = s*tmp + c*b(d+1,k)
    end do
  end subroutine d_rotation_times_tbc

  ! Rotation times truncated band matrix with aligned columns.
  ! First l and last p columns are not modified.
  subroutine f_z_rotation_times_tbc(c,s,b,n,lbw,ubw,l,p,j)
    complex(kind=dp), dimension(:,:), intent(inout) :: b
    complex(kind=dp), intent(in) :: s
    real(kind=dp), intent(in) :: c
    integer(kind=int32), intent(in) :: j,n, ubw,lbw, l, p
    !
    complex(kind=dp) :: tmp
    integer(kind=int32) :: k, k0, k1, d
    ! first and last columns operated on, given by the
    ! intersection of [l+1,n-p] and [max(1,j+1-lbw),min(n,j+ubw)]
    ! which is [maxl+1,max(1,j+1-lbw)), min(n-p,min(n,j+ubw))]
    k0=max(l+1,j+1-lbw)
    k1=min(n-p,j+ubw)
    ! loop over relevant columns in b.
    do k=k0,k1
       d=j-k+ubw+1
       tmp = b(d,k)
       b(d,k)=c*tmp - conjg(s)*b(d+1,k)
       b(d+1,k) = s*tmp + c*b(d+1,k)
    end do
  end subroutine f_z_rotation_times_tbc

  ! Rotation times truncated band matrix with aligned columns.
  ! First l and last p columns are not modified.
  subroutine z_rotation_times_tbc(r,b,n,lbw,ubw,l,p,j)
    complex(kind=dp), dimension(:,:), intent(inout) :: b
    type(z_rotation), intent(in) :: r
    integer(kind=int32), intent(in) :: j,n, ubw,lbw, l, p
    !
    complex(kind=dp) :: s, tmp
    real(kind=dp) :: c
    integer(kind=int32) :: k, k0, k1, d
    c=r%cosine; s=r%sine
    ! first and last columns operated on, given by the
    ! intersection of [l+1,n-p] and [max(1,j+1-lbw),min(n,j+ubw)]
    ! which is [maxl+1,max(1,j+1-lbw)), min(n-p,min(n,j+ubw))]
    k0=max(l+1,j+1-lbw)
    k1=min(n-p,j+ubw)
    ! loop over relevant columns in b.
    do k=k0,k1
       d=j-k+ubw+1
       tmp = b(d,k)
       b(d,k)=c*tmp - conjg(s)*b(d+1,k)
       b(d+1,k) = s*tmp + c*b(d+1,k)
    end do
  end subroutine z_rotation_times_tbc

  ! apply a rotation to columns k and k+1 of a truncated band matrix with
  ! aligned columns. first l and last p rows are not modified.
  subroutine f_d_tbc_times_rotation(b,n,lbw,ubw,l,p,c,s,k)
    real(kind=dp), dimension(:,:), intent(inout) :: b
    real(kind=dp), intent(in) :: c, s
    integer(kind=int32), intent(in) :: k,n, ubw,lbw, l, p
    !
    real(kind=dp) :: tmp
    integer(kind=int32) :: d, j0, j1, j
    ! first and last rows operated on are given by the intersection
    ! of [l+1,n-p] with [max(k-ubw+1,1),min(k+lbw,n)] which is
    ! [max(l+1,k-ubw+1),min(k+lbw,n-p)]
    j0=max(l+1,k-ubw+1)
    j1=min(k+lbw,n-p)
    do j=j0,j1
       d=j-k+ubw+1
       tmp=b(d,k)
       b(d,k)=c*tmp+s*b(d-1,k+1)
       b(d-1,k+1)=-s*tmp+c*b(d-1,k+1)
    end do
  end subroutine f_d_tbc_times_rotation

  ! apply a rotation to columns k and k+1 of a truncated band matrix with
  ! aligned columns. first l and last p rows are not modified.
  subroutine d_tbc_times_rotation(b,n,lbw,ubw,l,p,r,k)
    real(kind=dp), dimension(:,:), intent(inout) :: b
    type(d_rotation), intent(in) :: r
    integer(kind=int32), intent(in) :: k,n, ubw,lbw, l, p
    !
    real(kind=dp) :: c, s, tmp
    integer(kind=int32) :: d, j0, j1, j
    c=r%cosine; s=r%sine
    ! first and last rows operated on are given by the intersection
    ! of [l+1,n-p] with [max(k-ubw+1,1),min(k+lbw,n)] which is
    ! [max(l+1,k-ubw+1),min(k+lbw,n-p)]
    j0=max(l+1,k-ubw+1)
    j1=min(k+lbw,n-p)
    do j=j0,j1
       d=j-k+ubw+1
       tmp=b(d,k)
       b(d,k)=c*tmp+s*b(d-1,k+1)
       b(d-1,k+1)=-s*tmp+c*b(d-1,k+1)
    end do
  end subroutine d_tbc_times_rotation

  subroutine f_z_tbc_times_rotation(b,n,lbw,ubw,l,p,c,s,k)
    complex(kind=dp), dimension(:,:), intent(inout) :: b
    complex(kind=dp), intent(in) :: s
    real(kind=dp), intent(in) :: c
    integer(kind=int32), intent(in) :: k,n, ubw,lbw, l, p
    !
    complex(kind=dp) :: tmp
    integer(kind=int32) :: d, j0, j1, j
    ! first and last rows operated on are given by the intersection
    ! of [l+1,n-p] with [max(k-ubw+1,1),min(k+lbw,n)] which is
    ! [max(l+1,k-ubw+1),min(k+lbw,n-p)]
    j0=max(l+1,k-ubw+1)
    j1=min(k+lbw,n-p)
    do j=j0,j1
       d=j-k+ubw+1
       tmp=b(d,k)
       b(d,k)=c*tmp+s*b(d-1,k+1)
       b(d-1,k+1)=-conjg(s)*tmp+c*b(d-1,k+1)
    end do
  end subroutine f_z_tbc_times_rotation

  subroutine z_tbc_times_rotation(b,n,lbw,ubw,l,p,r,k)
    complex(kind=dp), dimension(:,:), intent(inout) :: b
    type(z_rotation), intent(in) :: r
    integer(kind=int32), intent(in) :: k,n, ubw,lbw, l, p
    !
    complex(kind=dp) :: s, tmp
    real(kind=dp) :: c
    integer(kind=int32) :: d, j0, j1, j
    c=r%cosine; s=r%sine
    ! first and last rows operated on are given by the intersection
    ! of [l+1,n-p] with [max(k-ubw+1,1),min(k+lbw,n)] which is
    ! [max(l+1,k-ubw+1),min(k+lbw,n-p)]
    j0=max(l+1,k-ubw+1)
    j1=min(k+lbw,n-p)
    do j=j0,j1
       d=j-k+ubw+1
       tmp=b(d,k)
       b(d,k)=c*tmp+s*b(d-1,k+1)
       b(d-1,k+1)=-conjg(s)*tmp+c*b(d-1,k+1)
    end do
  end subroutine z_tbc_times_rotation

  ! rotations for band matrices with aligned rows.
  !
  ! apply a rotation to columns k and k+1 of a truncated band matrix with
  ! aligned rows. The rotation is not applied to the first l or last p rows.
  subroutine f_d_tbr_times_rotation(b,m,lbw,ubw,l,p,c,s,k)
    real(kind=dp), dimension(:, :), intent(inout) :: b
    real(kind=dp), intent(in) :: c, s
    integer(kind=int32), intent(in) :: k,m, ubw,lbw, l,p
    !
    real(kind=dp) :: tmp
    integer(kind=int32) :: j, j0, j1, d
    ! first and last rows operated on are given by
    ! the intersection of [l+1,m-p] and [max(k-ubw+1,1),min(k+lbw,m)]
    ! which is [max(l+1,k-ubw+1), min(k+lbw,m-p)]
    j0=max(l+1,k-ubw+1)
    j1=min(k+lbw,m-p)
    do j=j0,j1
       d=k-j+lbw+1
       tmp=b(j,d)
       b(j,d)=c*tmp+s*b(j,d+1)
       b(j,d+1)=-s*tmp+c*b(j,d+1)
    end do
  end subroutine f_d_tbr_times_rotation

  ! apply a rotation to columns k and k+1 of a truncated band matrix with
  ! aligned rows. The rotation is not applied to the first l or last p rows.
  subroutine d_tbr_times_rotation(b,m,lbw,ubw,l,p,r,k)
    real(kind=dp), dimension(:, :), intent(inout) :: b
    type(d_rotation), intent(in) :: r
    integer(kind=int32), intent(in) :: k,m, ubw,lbw, l,p
    !
    real(kind=dp) :: c, s, tmp
    integer(kind=int32) :: j, j0, j1, d
    c=r%cosine; s=r%sine
    ! first and last rows operated on are given by
    ! the intersection of [l+1,m-p] and [max(k-ubw+1,1),min(k+lbw,m)]
    ! which is [max(l+1,k-ubw+1), min(k+lbw,m-p)]
    j0=max(l+1,k-ubw+1)
    j1=min(k+lbw,m-p)
    do j=j0,j1
       d=k-j+lbw+1
       tmp=b(j,d)
       b(j,d)=c*tmp+s*b(j,d+1)
       b(j,d+1)=-s*tmp+c*b(j,d+1)
    end do
  end subroutine d_tbr_times_rotation

  subroutine f_z_tbr_times_rotation(b,m,lbw,ubw,l,p,c,s,k)
    complex(kind=dp), dimension(:,:), intent(inout) :: b
    complex(kind=dp), intent(in) :: s
    real(kind=dp), intent(in) :: c
    integer(kind=int32), intent(in) :: k,m, ubw,lbw, l, p
    !
    complex(kind=dp) :: tmp
    integer(kind=int32) :: j, j0, j1, d
    ! first and last rows operated on are given by
    ! the intersection of [l+1,m-p] and [max(k-ubw+1,1),min(k+lbw,m)]
    ! which is [max(l+1,k-ubw+1), min(k+lbw,m-p)]
    j0=max(l+1,k-ubw+1)
    j1=min(k+lbw,m-p)
    do j=j0,j1
       d=k-j+lbw+1
       tmp=b(j,d)
       b(j,d)=c*tmp+s*b(j,d+1)
       b(j,d+1)=-conjg(s)*tmp+c*b(j,d+1)
    end do
  end subroutine f_z_tbr_times_rotation

  subroutine z_tbr_times_rotation(b,m,lbw,ubw,l,p,r,k)
    complex(kind=dp), dimension(:,:), intent(inout) :: b
    type(z_rotation), intent(in) :: r
    integer(kind=int32), intent(in) :: k,m, ubw,lbw, l, p
    !
    complex(kind=dp) :: s, tmp
    real(kind=dp) :: c
    integer(kind=int32) :: j, j0, j1, d
    c=r%cosine; s=r%sine
    ! first and last rows operated on are given by
    ! the intersection of [l+1,m-p] and [max(k-ubw+1,1),min(k+lbw,m)]
    ! which is [max(l+1,k-ubw+1), min(k+lbw,m-p)]
    j0=max(l+1,k-ubw+1)
    j1=min(k+lbw,m-p)
    do j=j0,j1
       d=k-j+lbw+1
       tmp=b(j,d)
       b(j,d)=c*tmp+s*b(j,d+1)
       b(j,d+1)=-conjg(s)*tmp+c*b(j,d+1)
    end do
  end subroutine z_tbr_times_rotation

  ! A rotation times a truncated band matrix with aligned rows.  The
  ! rotation is not applied to the first l or last p columns.
  subroutine f_d_rotation_times_tbr(c,s,b,m,lbw,ubw,l,p,j)
    real(kind=dp), dimension(:,:), intent(inout) :: b
    real(kind=dp), intent(in) :: c, s
    integer(kind=int32), intent(in) :: j,m, ubw,lbw, l, p
    !
    real(kind=dp) :: tmp
    integer(kind=int32) :: d, k0,k1,k
    ! first and last columns operated on, given by the
    ! intersection of [l+1,m-p] and [max(1,j+1-lbw),min(m,j+ubw)]
    ! which is [maxl+1,max(1,j+1-lbw)), min(m-p,min(n,j+ubw))]
    k0=max(l+1,j+1-lbw)
    k1=min(m-p,j+ubw)
    do k=k0,k1
       d=k-j+lbw+1
       tmp=b(j,d)
       b(j,d)=c*tmp-s*b(j+1,d-1)
       b(j+1,d-1)=s*tmp+c*b(j+1,d-1)
    end do
  end subroutine f_d_rotation_times_tbr

  ! A rotation times a truncated band matrix with aligned rows.  The
  ! rotation is not applied to the first l or last p columns.
  subroutine d_rotation_times_tbr(r,b,m,lbw,ubw,l,p,j)
    real(kind=dp), dimension(:,:), intent(inout) :: b
    type(d_rotation), intent(in) :: r
    integer(kind=int32), intent(in) :: j,m, ubw,lbw, l, p
    !
    real(kind=dp) :: c, s, tmp
    integer(kind=int32) :: d, k0,k1,k
    c=r%cosine; s=r%sine
    ! first and last columns operated on, given by the
    ! intersection of [l+1,m-p] and [max(1,j+1-lbw),min(m,j+ubw)]
    ! which is [maxl+1,max(1,j+1-lbw)), min(m-p,min(n,j+ubw))]
    k0=max(l+1,j+1-lbw)
    k1=min(m-p,j+ubw)
    do k=k0,k1
       d=k-j+lbw+1
       tmp=b(j,d)
       b(j,d)=c*tmp-s*b(j+1,d-1)
       b(j+1,d-1)=s*tmp+c*b(j+1,d-1)
    end do
  end subroutine d_rotation_times_tbr

  subroutine f_z_rotation_times_tbr(c,s,b,m,lbw,ubw,l,p,j)
    complex(kind=dp), dimension(:,:), intent(inout) :: b
    complex(kind=dp), intent(in) :: s
    real(kind=dp), intent(in) :: c
    integer(kind=int32), intent(in) :: j,m, ubw,lbw, l,p
    !
    complex(kind=dp) :: tmp
    integer(kind=int32) :: d, k0,k1,k
    ! first and last columns operated on, given by the
    ! intersection of [l+1,m-p] and [max(1,j+1-lbw),min(m,j+ubw)]
    ! which is [maxl+1,max(1,j+1-lbw)), min(m-p,min(n,j+ubw))]
    k0=max(l+1,j+1-lbw)
    k1=min(m-p,j+ubw)
    do k=k0,k1
       d=k-j+lbw+1
       tmp=b(j,d)
       b(j,d)=c*tmp-conjg(s)*b(j+1,d-1)
       b(j+1,d-1)=s*tmp+c*b(j+1,d-1)
    end do
  end subroutine f_z_rotation_times_tbr

  subroutine z_rotation_times_tbr(r,b,m,lbw,ubw,l,p,j)
    complex(kind=dp), dimension(:,:), intent(inout) :: b
    type(z_rotation), intent(in) :: r
    integer(kind=int32), intent(in) :: j,m, ubw,lbw, l,p
    !
    complex(kind=dp) :: s, tmp
    real(kind=dp) :: c
    integer(kind=int32) :: d, k0,k1,k
    c=r%cosine; s=r%sine
    ! first and last columns operated on, given by the
    ! intersection of [l+1,m-p] and [max(1,j+1-lbw),min(m,j+ubw)]
    ! which is [maxl+1,max(1,j+1-lbw)), min(m-p,min(n,j+ubw))]
    k0=max(l+1,j+1-lbw)
    k1=min(m-p,j+ubw)
    do k=k0,k1
       d=k-j+lbw+1
       tmp=b(j,d)
       b(j,d)=c*tmp-conjg(s)*b(j+1,d-1)
       b(j+1,d-1)=s*tmp+c*b(j+1,d-1)
    end do
  end subroutine z_rotation_times_tbr

  subroutine d_bc_to_br(bc,br,lbw,ubw)
    real(kind=dp), dimension(:,:), intent(in) :: bc
    real(kind=dp), dimension(:,:), intent(out) :: br
    integer(kind=int32), intent(in) :: lbw, ubw

    integer(kind=int32) :: d, j, n
    n=size(bc,2)
    do d=1,lbw+1
       do j=1,n-lbw+d-1
          br(lbw-d+j+1,d)=bc(ubw+lbw+2-d,j)
       end do
    end do
    do d=lbw+2,lbw+ubw+1
       do j=1,n-d+lbw+1
          br(j,d)=bc(ubw+lbw+2-d,j+d-lbw-1)
       end do
    end do
  end subroutine d_bc_to_br

  subroutine z_bc_to_br(bc,br,lbw,ubw)
    complex(kind=dp), dimension(:,:), intent(in) :: bc
    complex(kind=dp), dimension(:,:), intent(out) :: br
    integer(kind=int32), intent(in) :: lbw, ubw

    integer(kind=int32) :: d, j, n
    n=size(bc,2)
    do d=1,lbw+1
       do j=1,n-lbw+d-1
          br(lbw-d+j+1,d)=bc(ubw+lbw+2-d,j)
       end do
    end do
    do d=lbw+2,lbw+ubw+1
       do j=1,n-d+lbw+1
          br(j,d)=bc(ubw+lbw+2-d,j+d-lbw-1)
       end do
    end do
  end subroutine z_bc_to_br

  subroutine d_br_to_bc(br,bc,lbw,ubw)
    real(kind=dp), dimension(:,:), intent(out) :: bc
    real(kind=dp), dimension(:,:), intent(in) :: br
    integer(kind=int32), intent(in) :: lbw, ubw

    integer(kind=int32) :: d, j, n

    n=size(bc,2)
    do d=1,lbw+1
       do j=1,n-lbw+d-1
          bc(ubw+lbw+2-d,j)=br(lbw-d+j+1,d)
       end do
    end do
    do d=lbw+2,lbw+ubw+1
       do j=1,n-d+lbw+1
          bc(ubw+lbw+2-d,j+d-lbw-1)=br(j,d)
       end do
    end do
  end subroutine d_br_to_bc

  subroutine z_br_to_bc(br,bc,lbw,ubw)
    complex(kind=dp), dimension(:,:), intent(out) :: bc
    complex(kind=dp), dimension(:,:), intent(in) :: br
    integer(kind=int32), intent(in) :: lbw, ubw

    integer(kind=int32) :: d, j, n
    n=size(bc,2)
    do d=1,lbw+1
       do j=1,n-lbw+d-1
          bc(ubw+lbw+2-d,j)=br(lbw-d+j+1,d)
       end do
    end do
    do d=lbw+2,lbw+ubw+1
       do j=1,n-d+lbw+1
          bc(ubw+lbw+2-d,j+d-lbw-1)=br(j,d)
       end do
    end do
  end subroutine z_br_to_bc

  subroutine d_extract_diagonals_bc(a, n, b, lbw, ubw, lbwmax, ubwmax)
    real(kind=dp), target, dimension(n,n), intent(in) :: a
    real(kind=dp), dimension(lbwmax+ubwmax+1,n), intent(out) :: b
    integer(kind=int32), intent(in) :: n, lbw, ubw, lbwmax, ubwmax
    !
    integer(kind=int32) :: k, d
    ! put diagonals in b
    b=0.0_dp
    do d=1,ubw+1
       do k=ubw-d+2,n
          b(d,k) = a(d+k-ubw-1,k)
       end do
    end do
    do d=ubw+2, ubw+lbw+1
       do k=1,n-d+ubw+1
          b(d,k) = a(d+k-ubw-1,k)
       end do
    end do
  end subroutine d_extract_diagonals_bc

  subroutine z_extract_diagonals_bc(a, n, b, lbw, ubw, lbwmax, ubwmax)
    complex(kind=dp), target, dimension(n,n), intent(in) :: a
    complex(kind=dp), dimension(lbwmax+ubwmax+1,n), intent(out) :: b
    integer(kind=int32), intent(in) :: n, lbw, ubw, lbwmax, ubwmax
    !
    integer(kind=int32) :: k, d
    ! put diagonals in b
    b=0.0_dp
    do d=1,ubw+1
       do k=ubw-d+2,n
          b(d,k) = a(d+k-ubw-1,k)
       end do
    end do
    do d=ubw+2, ubw+lbw+1
       do k=1,n-d+ubw+1
          b(d,k) = a(d+k-ubw-1,k)
       end do
    end do
  end subroutine z_extract_diagonals_bc

  subroutine d_extract_diagonals_br(a, n, b, lbw, ubw, lbwmax, ubwmax)
    real(kind=dp), target, dimension(n,n), intent(in) :: a
    real(kind=dp), dimension(n,lbwmax+ubwmax+1), intent(out) :: b
    integer(kind=int32), intent(in) :: n, lbw, ubw, lbwmax, ubwmax
    !
    integer(kind=int32) :: j, d
    ! put diagonals in b
    b=0.0_dp
    do d=1,lbw+1
       do j=lbw-d+2, n
          b(j,d)=a(j,d+j-lbw-1)
       end do
    end do
    do d=lbw+2,ubw+lbw+1
       do j=1, n-d+lbw+1
          b(j,d)=a(j,d+j-lbw-1)
       end do
    end do
  end subroutine d_extract_diagonals_br

  subroutine z_extract_diagonals_br(a, n, b, lbw, ubw, lbwmax, ubwmax)
    complex(kind=dp), target, dimension(n,n), intent(in) :: a
    complex(kind=dp), dimension(n,lbwmax+ubwmax+1), intent(out) :: b
    integer(kind=int32), intent(in) :: n, lbw, ubw, lbwmax, ubwmax
    !
    integer(kind=int32) :: j, d
    ! put diagonals in b
    b=0.0_dp
    do d=1,lbw+1
       do j=lbw-d+2, n
          b(j,d)=a(j,d+j-lbw-1)
       end do
    end do
    do d=lbw+2,ubw+lbw+1
       do j=1, n-d+lbw+1
          b(j,d)=a(j,d+j-lbw-1)
       end do
    end do
  end subroutine z_extract_diagonals_br

  ! Impose a particular upper and lower variable bandwidth on a band
  ! matrix.  Used for generating random band matrices for tests.
  subroutine d_truncate_profile_br(br,n,lbw,lbwmax,ubwmax,lower,upper)
    real(kind=dp), dimension(n,lbwmax+ubwmax+1), intent(inout) :: br
    integer(kind=int32), intent(in) :: n, lbw, lbwmax, ubwmax
    integer(kind=int32), dimension(n), intent(in), optional :: lower, upper

    integer(kind=int32) :: j

    if (present(upper)) then
       do j=1,n
          br(j,lbw+2+upper(j):lbwmax+ubwmax+1) = 0.0_dp
       end do
    end if

    if (present(lower)) then
       do j=1,n
          br(j,1:lbw-lower(j)) = 0.0_dp
       end do
    end if

  end subroutine d_truncate_profile_br

  subroutine z_truncate_profile_br(br,n,lbw,lbwmax,ubwmax,lower,upper)
    complex(kind=dp), dimension(n,lbwmax+ubwmax+1), intent(inout) :: br
    integer(kind=int32), intent(in) :: n, lbw, lbwmax, ubwmax
    integer(kind=int32), dimension(n), intent(in), optional :: lower, upper

    integer(kind=int32) :: j

    if (present(upper)) then
       do j=1,n
          br(j,lbw+2+upper(j):lbwmax+ubwmax+1) = (0.0_dp,0.0_dp)
       end do
    end if

    if (present(lower)) then
       do j=1,n
          br(j,1:lbw-lower(j)) = (0.0_dp,0.0_dp)
       end do
    end if
  end subroutine z_truncate_profile_br

  subroutine d_truncate_profile_bc(bc,n,ubw,lbwmax,ubwmax,lower,upper)
    real(kind=dp), dimension(lbwmax+ubwmax+1,n), intent(inout) :: bc
    integer(kind=int32), intent(in) :: n, ubw, lbwmax, ubwmax
    integer(kind=int32), dimension(n), intent(in), optional :: lower, upper

    integer(kind=int32) :: j

    if (present(upper)) then
       do j=1,n
          bc(1:ubw-upper(j),j)=0.0_dp
       end do
    end if

    if (present(lower)) then
       do j=1,n
          bc(ubw+2+lower(j):lbwmax+ubwmax+1,j)=0.0_dp
       end do
    end if

  end subroutine d_truncate_profile_bc

  subroutine z_truncate_profile_bc(bc,n,ubw,lbwmax,ubwmax,lower,upper)
    complex(kind=dp), dimension(lbwmax+ubwmax+1,n), intent(inout) :: bc
    integer(kind=int32), intent(in) :: n, ubw, lbwmax, ubwmax
    integer(kind=int32), dimension(n), intent(in), optional :: lower, upper

    integer(kind=int32) :: j

    if (present(upper)) then
       do j=1,n
          bc(1:ubw-upper(j),j)=(0.0_dp,0.0_dp)
       end do
    end if

    if (present(lower)) then
       do j=1,n
          bc(ubw+2+lower(j):lbwmax+ubwmax+1,j)=(0.0_dp,0.0_dp)
       end do
    end if

  end subroutine z_truncate_profile_bc

  real(kind=dp) function d_band_norm_bc(bc,n,lbw,ubw,lbwmax, &
       ubwmax,d0,d1) result(x)
    real(kind=dp), dimension(lbwmax+ubwmax+1,n), intent(in) :: bc
    integer(kind=int32), intent(in) :: n, lbw, ubw, lbwmax, ubwmax, d0, d1

    integer(kind=int32) :: k,l,l0,l1
    real(kind=dp) :: y, tmp

    l0 = ubw+1-d1
    l1 = ubw+1-d0

    y=0.0_dp
    ! maxabs for the specified band.
    do l=max(1,l0),min(ubw+1,l1)
       do k=ubw-l+2,n
          tmp=abs(bc(l,k))
          if (tmp > y) y=tmp
       end do
    end do
    do l=max(l0,ubw+2), min(ubw+lbw+1,l1)
       do k=1,n-l+ubw+1
          tmp=abs(bc(l,k))
          if (tmp > y) y=tmp
       end do
    end do
    x=0.0_dp;
    if (y==0.0_dp) return
    y=(2.0_dp)**(floor(log(y)/log(2.0_dp)))
    y=1.0_dp/y
    if (y > 1e150_dp .or. y < 1e-150_dp) then
       do l=max(1,l0),min(ubw+1,l1)
          do k=ubw-l+2,n
             tmp=bc(l,k)*y
             x=x+tmp*tmp
          end do
       end do
       do l=max(l0,ubw+2), min(ubw+lbw+1,l1)
          do k=1,n-l+ubw+1
             tmp=bc(l,k)*y
             x=x+tmp*tmp
          end do
       end do
       x=sqrt(x)/y
    else
       do l=max(1,l0),min(ubw+1,l1)
          do k=ubw-l+2,n
             tmp=bc(l,k)
             x=x+tmp*tmp
          end do
       end do
       do l=max(l0,ubw+2), min(ubw+lbw+1,l1)
          do k=1,n-l+ubw+1
             tmp=bc(l,k)
             x=x+tmp*tmp
          end do
       end do
       x=sqrt(x)
    end if
  end function d_band_norm_bc

  real(kind=dp) function z_band_norm_bc(bc,n,lbw,ubw,lbwmax, &
       ubwmax,d0,d1) result(x)
    complex(kind=dp), dimension(lbwmax+ubwmax+1,n), intent(in) :: bc
    integer(kind=int32), intent(in) :: n, lbw, ubw, lbwmax, ubwmax, d0, d1

    integer(kind=int32) :: k,l,l0,l1
    real(kind=dp) :: y, tmp

    l0 = ubw+1-d1
    l1 = ubw+1-d0

    y=0.0_dp
    ! maxabs for the specified band.
    do l=max(1,l0),min(ubw+1,l1)
       do k=ubw-l+2,n
          tmp=abs(bc(l,k))
          if (tmp > y) y=tmp
       end do
    end do
    do l=max(l0,ubw+2), min(ubw+lbw+1,l1)
       do k=1,n-l+ubw+1
          tmp=abs(bc(l,k))
          if (tmp > y) y=tmp
       end do
    end do

    x=0.0_dp;
    if (y==0.0_dp) return
    y=(2.0_dp)**(floor(log(y)/log(2.0_dp)))
    y=1.0_dp/y
    if (y > 1e150_dp .or. y < 1e-150_dp) then
       do l=max(1,l0),min(ubw+1,l1)
          do k=ubw-l+2,n
             tmp=abs(bc(l,k)*y)
             x=x+tmp*tmp
          end do
       end do
       do l=max(l0,ubw+2), min(ubw+lbw+1,l1)
          do k=1,n-l+ubw+1
             tmp=abs(bc(l,k)*y)
             x=x+tmp*tmp
          end do
       end do
       x=sqrt(x)/y
    else
       do l=max(1,l0),min(ubw+1,l1)
          do k=ubw-l+2,n
             tmp=abs(bc(l,k))
             x=x+tmp*tmp
          end do
       end do
       do l=max(l0,ubw+2), min(ubw+lbw+1,l1)
          do k=1,n-l+ubw+1
             tmp=abs(bc(l,k))
             x=x+tmp*tmp
          end do
       end do
       x=sqrt(x)
    end if
  end function z_band_norm_bc

  real(kind=dp) function d_band_norm_br(br,n,lbw,ubw,lbwmax, &
       ubwmax,d0,d1) result(x)
    real(kind=dp), dimension(n,lbwmax+ubwmax+1), intent(in) :: br
    integer(kind=int32), intent(in) :: n, lbw, ubw, lbwmax, ubwmax, d0, d1

    integer(kind=int32) :: j,l,l0,l1
    real(kind=dp) :: y, tmp

    l0 = d0+lbw+1
    l1 = d1+lbw+1

    y=0.0_dp
    ! maxabs for the specified band.

    do l=max(1,l0),min(lbw+1,l1)
       do j=lbw-l+2,n
          tmp=abs(br(j,l))
          if (tmp > y) y=tmp
       end do
    end do
    do l=max(l0,lbw+2), min(ubw+lbw+1,l1)
       do j=1,n-l+lbw+1
          tmp=abs(br(j,l))
          if (tmp > y) y=tmp
       end do
    end do

    x=0.0_dp;
    if (y==0.0_dp) return
    y=(2.0_dp)**(floor(log(y)/log(2.0_dp)))
    y=1.0_dp/y
    if (y > 1e150_dp .or. y < 1e-150_dp) then
       do l=max(1,l0),min(lbw+1,l1)
          do j=lbw-l+2,n
             tmp=br(j,l)*y
             x=x+tmp*tmp
          end do
       end do
       do l=max(l0,lbw+2), min(ubw+lbw+1,l1)
          do j=1,n-l+lbw+1
             tmp=br(j,l)*y
             x=x+tmp*tmp
          end do
       end do
       x=sqrt(x)/y
    else
       do l=max(1,l0),min(lbw+1,l1)
          do j=lbw-l+2,n
             tmp=br(j,l)
             x=x+tmp*tmp
          end do
       end do
       do l=max(l0,lbw+2), min(ubw+lbw+1,l1)
          do j=1,n-l+lbw+1
             tmp=br(j,l)
             x=x+tmp*tmp
          end do
       end do
       x=sqrt(x)
    end if
  end function d_band_norm_br

  real(kind=dp) function z_band_norm_br(br,n,lbw,ubw,lbwmax, &
       ubwmax,d0,d1) result(x)
    complex(kind=dp), dimension(n,lbwmax+ubwmax+1), intent(in) :: br
    integer(kind=int32), intent(in) :: n, lbw, ubw, lbwmax, ubwmax, d0, d1

    integer(kind=int32) :: j,l,l0,l1
    real(kind=dp) :: y, tmp

    l0 = d0+lbw+1
    l1 = d1+lbw+1

    y=0.0_dp
    ! maxabs for the specified band.
    do l=max(1,l0),min(lbw+1,l1)
       do j=lbw-l+2,n
          tmp=abs(br(j,l))
          if (tmp > y) y=tmp
       end do
    end do
    do l=max(l0,lbw+2), min(ubw+lbw+1,l1)
       do j=1,n-l+lbw+1
          tmp=abs(br(j,l))
          if (tmp > y) y=tmp
       end do
    end do

    x=0.0_dp;
    if (y==0.0_dp) return
    y=(2.0_dp)**(floor(log(y)/log(2.0_dp)))
    y=1.0_dp/y
    if (y > 1e150_dp .or. y < 1e-150_dp) then
       do l=max(1,l0),min(lbw+1,l1)
          do j=lbw-l+2,n
             tmp=abs(br(j,l)*y)
             x=x+tmp*tmp
          end do
       end do
       do l=max(l0,lbw+2), min(ubw+lbw+1,l1)
          do j=1,n-l+lbw+1
             tmp=abs(br(j,l)*y)
             x=x+tmp*tmp
          end do
       end do
       x=sqrt(x)/y
    else
       do l=max(1,l0),min(lbw+1,l1)
          do j=lbw-l+2,n
             tmp=abs(br(j,l))
             x=x+tmp*tmp
          end do
       end do
       do l=max(l0,lbw+2), min(ubw+lbw+1,l1)
          do j=1,n-l+lbw+1
             tmp=abs(br(j,l))
             x=x+tmp*tmp
          end do
       end do
       x=sqrt(x)
    end if
  end function z_band_norm_br

end module mod_band_types

max_stars_repo_name: m-a-stewart/BandGivensWeight

max_issues_repo_name: m-a-stewart/BandGivensWeight

max_forks_repo_name: m-a-stewart/BandGivensWeight


text:
*   Subroutines to integrate the RGEs for the gauge and Yukawa couplings

      SUBROUTINE
     .       ODEINT(YSTART,NVAR,X1,X2,EPS,DERIVS,RKQS,IFAIL)

*   Driver subroutine to integrate Ordinary Differential Equations
*   using the Runge-Kutta 5th order method with adaptative step size
*   (from Numerical Recipes)
*   IFAIL=1 stepsize smaller than minimum in ODEINT
*   IFAIL=2 too many steps in ODEINT
*   IFAIL=3 stepsize underflow in RKQS
*   IFAIL=4 max(|dy(i)/dx|*(1+x)/(1+|y(i)|)) > 1/eps^2

      IMPLICIT NONE

      INTEGER NVAR,NMAX,IFAIL,MAXSTP,NSTP,I
      PARAMETER(MAXSTP=10000,NMAX=500)

      DOUBLE PRECISION YSTART(NVAR),Y(NMAX),YSCAL(NMAX)
      DOUBLE PRECISION DYDX(NMAX),X,X1,X2,EPS,TINY
      DOUBLE PRECISION H1,HMIN,H,HDID,HNEXT

      EXTERNAL DERIVS,RKQS

      IFAIL=0
      TINY=EPS**4
      H1=DSQRT(EPS)*DABS(X2-X1)
      HMIN=EPS**2*DABS(X2-X1)

      X=X1
      H=SIGN(H1,X2-X1)
      DO I=1,NVAR
       Y(I)=YSTART(I)
      ENDDO

      DO NSTP=1,MAXSTP

       CALL DERIVS(NVAR,X,Y,DYDX)

       DO I=1,NVAR
        IF(DABS(DYDX(I))*(1d0+X)/(1d0+DABS(Y(I))).GT.EPS**(-2))IFAIL=4
       ENDDO
       IF(IFAIL.GT.0)RETURN

       DO I=1,NVAR
        YSCAL(I)=DABS(Y(I))+DABS(H*DYDX(I))+TINY
       ENDDO

       CALL RKQS(Y,DYDX,NVAR,X,X2,H,EPS,YSCAL,HDID,HNEXT,DERIVS,IFAIL)
       IF(IFAIL.GT.0)RETURN

       IF(DABS(Y(2)-5d0/3d0*Y(1)).LT.EPS)THEN
        X2=X
        DO I=1,NVAR
         YSTART(I)=Y(I)
        ENDDO
        RETURN
       ENDIF

       IF(ABS(HNEXT).LT.HMIN)THEN
        IFAIL=1
        RETURN
       ENDIF

       H=HNEXT

      ENDDO

      IFAIL=2

      RETURN
      END

      SUBROUTINE RKQS(Y,DYDX,N,X,X2,HTRY,EPS,YSCAL,HDID,
     .       HNEXT,DERIVS,IFAIL)

*   Stepper subroutine for ODEINT

      IMPLICIT NONE

      INTEGER IFAIL,I,N,NMAX
      PARAMETER(NMAX=500)

      DOUBLE PRECISION EPS,HDID,HNEXT,HTRY,X,X2,DYDX(N),Y(N),YSCAL(N)
      DOUBLE PRECISION ERRMAX,H,HTEMP,XNEW,YERR(NMAX),YTEMP(NMAX)
      DOUBLE PRECISION SAFETY,PGROW,PSHRNK,ERRCON,G,GTEMP

      EXTERNAL DERIVS

      SAFETY=.9d0
      PGROW=-.2d0
      PSHRNK=-.25d0
      ERRCON=(5d0/SAFETY)**(1d0/PGROW)
      H=HTRY
      G=Y(2)-5d0/3d0*Y(1)

1     CALL RKCK(Y,DYDX,N,X,H,YTEMP,YERR,DERIVS)

      ERRMAX=0d0
      DO I=1,N
       ERRMAX=MAX(ERRMAX,DABS(YERR(I)/YSCAL(I)))
      ENDDO
      ERRMAX=ERRMAX/EPS

      IF(ERRMAX.GT.1d0)THEN
       HTEMP=SAFETY*H*(ERRMAX**PSHRNK)
       H=SIGN(MAX(DABS(HTEMP),.1d0*DABS(H)),H)
       XNEW=X+H
       IF(XNEW.EQ.X)THEN
        IFAIL=3
        RETURN
       ENDIF
       GOTO 1
      ENDIF

      GTEMP=YTEMP(2)-5d0/3d0*YTEMP(1)
      IF(GTEMP.LT.-EPS)THEN
       IFAIL=-1
       X2=X+H
       H=G/(G-GTEMP)*H
       GOTO 1
      ENDIF

      HDID=H
      X=X+H
      DO I=1,N
       Y(I)=YTEMP(I)
      ENDDO

      IF(ERRMAX.GT.ERRCON)THEN
       HNEXT=SAFETY*H*(ERRMAX**PGROW)
      ELSE
        HNEXT=5d0*H
      ENDIF

      IF(X+HNEXT.GT.X2 .AND. IFAIL.EQ.-1) HNEXT=X2-X

      RETURN

      END

      SUBROUTINE RKCK(Y,DYDX,N,X,H,YOUT,YERR,DERIVS)

*   Algorithm subroutine for ODEINT

      IMPLICIT NONE

      INTEGER I,N,NMAX
      PARAMETER(NMAX=500)

      DOUBLE PRECISION H,X,DYDX(N),Y(N),YERR(N),YOUT(N)
      DOUBLE PRECISION AK2(NMAX),AK3(NMAX),AK4(NMAX)
      DOUBLE PRECISION AK5(NMAX),AK6(NMAX),YTEMP(NMAX)
      DOUBLE PRECISION A2,A3,A4,A5,A6,C1,C3,C4,C6
      DOUBLE PRECISION B21,B31,B32,B41,B42,B43,B51,B52,B53,B54
      DOUBLE PRECISION B61,B62,B63,B64,B65,DC1,DC3,DC4,DC5,DC6

      PARAMETER(A2=.2d0, A3=.3d0, A4=.6d0, A5=1d0, A6=.875d0,
     .       B21=.2d0, B31=3d0/40d0, B32=9d0/40d0, B41=.3d0, B42=-.9d0,
     .       B43=1.2d0, B51=-11d0/54d0, B52=2.5d0, B53=-70d0/27d0,
     .       B54=35d0/27d0, B61=1631d0/55296d0, B62=175d0/512d0,
     .       B63=575d0/13824d0, B64=44275d0/110592d0,
     .       B65=253d0/4096d0, C1=37d0/378d0, C3=250d0/621d0,
     .       C4=125d0/594d0, C6=512d0/1771d0, DC1=C1-2825d0/27648d0,
     .       DC3=C3-18575d0/48384d0, DC4=C4-13525d0/55296d0,
     .       DC5=-277d0/14336d0, DC6=C6-.25d0)

      EXTERNAL DERIVS

      DO I=1,N
       YTEMP(I)=Y(I)+B21*H*DYDX(I)
      ENDDO
      CALL DERIVS(N,X+A2*H,YTEMP,AK2)
      DO I=1,N
       YTEMP(I)=Y(I)+H*(B31*DYDX(I)+B32*AK2(I))
      ENDDO
      CALL DERIVS(N,X+A3*H,YTEMP,AK3)
      DO I=1,N
       YTEMP(I)=Y(I)+H*(B41*DYDX(I)+B42*AK2(I)+B43*AK3(I))
      ENDDO
      CALL DERIVS(N,X+A4*H,YTEMP,AK4)
      DO I=1,N
       YTEMP(I)=Y(I)+H*(B51*DYDX(I)+B52*AK2(I)+B53*AK3(I)+B54*AK4(I))
      ENDDO
      CALL DERIVS(N,X+A5*H,YTEMP,AK5)
      DO I=1,N
       YTEMP(I)=Y(I)+H*(B61*DYDX(I)+B62*AK2(I)+B63*AK3(I)+B64*AK4(I)+
     .  B65*AK5(I))
      ENDDO
      CALL DERIVS(N,X+A6*H,YTEMP,AK6)
      DO I=1,N
       YOUT(I)=Y(I)+H*(C1*DYDX(I)+C3*AK3(I)+C4*AK4(I)+C6*AK6(I))
      ENDDO
      DO I=1,N
       YERR(I)=H*(DC1*DYDX(I)+DC3*AK3(I)+DC4*AK4(I)+DC5*AK5(I)+DC6*
     .  AK6(I))
      ENDDO

      RETURN
      END

      SUBROUTINE DERIVS(N,X,Y,F)

*   2-loop Renormalization group equations for g1, g2, g3,
*   lambda, kappa, htop, hbot, htau to be integrated by ODEINT

      IMPLICIT NONE

      INTEGER N
      DOUBLE PRECISION X,Y(N),F(N),PI,c2
      DOUBLE PRECISION G1,G2,G3,L2,K2,HT2,HB2,HTAU2

      PI=4d0*DATAN(1d0)
      c2=1d0/(16d0*PI**2)

      G1=Y(1)
      G2=Y(2)
      G3=Y(3)
      L2=Y(4)
      K2=Y(5)
      HT2=Y(6)
      HB2=Y(7)
      HTAU2=Y(8)

      F(1)= 11d0*g1**2
     .    + c2*g1**2*(199d0/9d0*g1 + 9d0*g2 + 88d0/3d0*g3
     .    - 2d0*L2 - 26d0/3d0*HT2 - 14d0/3d0*HB2 - 6d0*HTAU2)

      F(2)= g2**2
     .    + c2*g2**2*(3d0*g1 + 25d0*g2 + 24d0*g3
     .    - 2d0*L2 - 6d0*HT2 - 6d0*HB2 - 2d0*HTAU2)

      F(3)= -3d0*g3**2
     .    + c2*g3**2*(11d0/3d0*g1 + 9d0*g2 + 14d0*g3
     .    - 4d0*HT2 - 4d0*HB2)

      F(4)= L2*(-g1 - 3d0*g2
     .    + 4d0*L2 + 2d0*K2 + 3d0*HT2 + 3d0*HB2 + HTAU2)
     .    + c2*L2*(L2*(2d0*g1 + 6d0*g2
     .    - 10d0*L2 - 12d0*K2 - 9d0*HT2 - 9d0*HB2 - 3d0*HTAU2)
     .    + 23d0/2d0*g1**2 + 15d0/2d0*g2**2 + 3d0*g1*g2
     .    + 4d0/3d0*g1*HT2 - 2d0/3d0*g1*HB2 + 2d0*g1*HTAU2
     .    + 16d0*g3*HT2 + 16d0*g3*HB2 - 3d0*HTAU2**2
     .    - 9d0*HT2**2 - 9d0*HB2**2 - 6d0*HT2*HB2 - 8d0*K2**2)

      F(5)= 6d0*K2*(L2 + K2)
     .    + c2*6d0*K2*(-4d0*K2**2 + L2*(g1 + 3d0*g2
     .    - 2d0*L2 - 4d0*K2 - 3d0*HT2 - 3d0*HB2 - HTAU2))

      F(6)= HT2*(-13d0/9d0*g1 - 3d0*g2 - 16d0/3d0*g3
     .    + L2 + 6d0*HT2 + HB2)
     .    + c2*HT2*(HT2*(2d0*g1 + 6d0*g2 + 16d0*g3
     .    - 3d0*L2 - 22d0*HT2 - 5d0*HB2)
     .    + 2d0/3d0*g1*HB2 - 3d0*L2**2 - 2d0*L2*K2
     .    - 4d0*L2*HB2 - 5d0*HB2**2 - HB2*HTAU2 - L2*HTAU2
     .    + 2743d0/162d0*g1**2 + 15d0/2d0*g2**2- 16d0/9d0*g3**2
     .    + 5d0/3d0*g1*g2 + 136d0/27d0*g1*g3 + 8d0*g2*g3)

      F(7)= HB2*(-7d0/9d0*g1 - 3d0*g2 - 16d0/3d0*g3
     .    + L2 + HT2 + 6d0*HB2 + HTAU2)
     .    + c2*HB2*(HB2*(2d0/3d0*g1 + 6d0*g2 + 16d0*g3
     .    - 3d0*L2 - 5d0*HT2 - 22d0*HB2 - 3d0*HTAU2)
     .    + 4d0/3d0*g1*HT2 - 3d0*L2**2 - 2d0*L2*K2
     .    - 4d0*L2*HT2 - 5d0*HT2**2 - 3d0*HTAU2**2 + 2d0*g1*HTAU2
     .    + 1435d0/162d0*g1**2 + 15d0/2d0*g2**2 - 16d0/9d0*g3**2
     .    + 5d0/3d0*g1*g2 + 40d0/27d0*g1*g3 + 8d0*g2*g3)

      F(8)= HTAU2*(-3d0*g1 - 3d0*g2
     .    + L2 + 3d0*HB2 + 4d0*HTAU2)
     .    + c2*HTAU2*(-10d0*HTAU2**2 - 9d0*HTAU2*HB2 - 9d0*HB2**2
     .    - 3d0*HB2*HT2 + HTAU2*(2d0*G1 + 6d0*G2)
     .    + HB2*(-2d0/3d0*G1 + 16*G3)
     .    - L2*(3d0*HTAU2 + 3d0*HT2 + 3d0*L2 + 2d0*K2)
     .    + 75d0/2d0*G1**2 + 3d0*G1*G2 + 15d0/2d0*G2**2)

      RETURN
      END