w3sbs1md.F90

Go to the documentation of this file.

 
 
#include "w3macros.h"
 
!/ ------------------------------------------------------------------- /
MODULE w3sbs1md
  !/
  !/                  +-----------------------------------+
  !/                  | WAVEWATCH III                SHOM |
  !/                  |            F. Ardhuin             |
  !/                  |                        FORTRAN 90 |
  !/                  | Last update :         14-Nov-2010 |
  !/                  +-----------------------------------+
  !/
  !/    15-Jul-2005 : Origination.                        ( version 3.07 )
  !/    23-Jun-2006 : Formatted for submitting code for   ( version 3.09 )
  !/                  inclusion in WAVEWATCH III.
  !/    10-May-2007 : adapt from version 2.22.SHOM        ( version 3.10.SHOM )
  !/    14-Nov-2010 : include scaling factor and clean up ( version 3.14 )
  !/
  !  1. Purpose :
  !
  !     This module computes a scattering term
  !     based on the theory by Ardhuin and Magne (JFM 2007)
  !
  !  2. Variables and types :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !     ----------------------------------------------------------------
  !
  !  3. Subroutines and functions :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !      W3SBS1    Subr. Public   bottom scattering
  !      INSBS1    Subr. Public   Corresponding initialization routine.
  !     ----------------------------------------------------------------
  !
  !  4. Subroutines and functions used :
  !
  !      Name      Type  Module   Description
  !     ----------------------------------------------------------------
  !      STRACE    Subr. W3SERVMD Subroutine tracing.
  !     ----------------------------------------------------------------
  !
  !  5. Remarks :
  !
  !
  !  6. Switches :
  !
  !     !/S  Enable subroutine tracing.
  !
  !  7. Source code :
  !/
  !/ ------------------------------------------------------------------- /
  !/
  !
  PUBLIC
  !/
  !/ Public variables
  !/
  REAL, DIMENSION(:,:), ALLOCATABLE   :: botspec
  INTEGER, PARAMETER :: nkscat = 30                 !number of wavenumbers
  DOUBLE PRECISION  ,DIMENSION(:,:,:) ,  ALLOCATABLE :: scatmatv !scattering matrices
  DOUBLE PRECISION  ,DIMENSION(:,:,:) ,  ALLOCATABLE :: scatmata !original matrix
  DOUBLE PRECISION  ,DIMENSION(:,:)   ,  ALLOCATABLE :: scatmatd
  CHARACTER(len=10)                   :: botspec_indicator
  INTEGER            :: nkbx, nkby
  REAL               :: dkbx, dkby, kwmin, kwmax
  REAL, PARAMETER    :: scattcutoff=0.
  REAL               :: curtx, curty
  !/
CONTAINS
 
  SUBROUTINE w3sbs1(A, CG, WN, DEPTH, CX1, CY1,      &
       TAUSCX, TAUSCY, S, D)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |            F. Ardhuin             |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         23-Jun-2006 |
    !/                  +-----------------------------------+
    !/
    !/    15-Jul-2005 : Origination.                        ( version 3.07 )
    !/    23-Jun-2006 : Formatted for submitting code for   ( version 3.09 )
    !/                  inclusion in WAVEWATCH III.
    !/
    !  1. Purpose :
    !
    !     Bottom scattering source term
    !
    !  2. Method :
    !
    !     Without current, goes through a diagonalization of the matrix
    !     problem  S(f,:) = M(f,:,:)**E(f,:)
    !     With current, integrates the source term along the resonant locus
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !  A         R.A.  I   Action density spectrum (1-D)
    !       CG        R.A.  I   Group velocities.
    !       WN        R.A.  I   Wavenumbers.
    !       DEPTH     Real  I   Mean water depth.
    !       S         R.A.  O   Source term (1-D version).
    !       D         R.A.  O   Diagonal term of derivative (1-D version).
    !       CX1-Y1    R.A.  I   Current components at ISEA.
    !       TAUSCX-Y  R.A.  I   Change of wave momentum due to scattering
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      STRACE    Subr. W3SERVMD Subroutine tracing.
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3SRCE    Subr. W3SRCEMD Source term integration.
    !      W3EXPO    Subr.   N/A    Point output post-processor.
    !      GXEXPO    Subr.   N/A    GrADS point output post-processor.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !       None.
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/S  Enable subroutine tracing.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants
    USE w3gdatmd, ONLY: nk, nth, nspec, sig, dth, dden, &
         ecos, esin, ec2, mapth, mapwn, &
         sig2, dsii
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    !/
    !
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    REAL, INTENT(IN)        :: CG(NK), WN(NK), DEPTH
    REAL, INTENT(IN)        :: A(NTH,NK)
    REAL, INTENT(IN)        :: CX1, CY1
    REAL, INTENT(OUT)       :: TAUSCX, TAUSCY
    REAL, INTENT(OUT)       :: S(NSPEC), D(NSPEC)
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER         :: ISPEC, IK, NSCUT, ITH, ITH2, i, j,iajust,iajust2
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    LOGICAL, SAVE   :: FIRST = .true.
    INTEGER         :: MATRICES = 0
 
    REAL            :: R1, R2, R3
 
    REAL            :: WN2(NSPEC, NTH), Ka(NSPEC),        &
         Kb(NSPEC, NTH), WNBOT(NSPEC, NTH), &
         B(NSPEC, NTH)
    REAL            :: kbotxi, kbotyi, xbk,   &
         ybk,integral, kbotx, kboty, count,count2
 
    INTEGER         :: ibk, jbk, ik2
    REAL            ::  SIGP,KU, KPU, CGK, CGPK, WN2i, xk2, Ap, kcutoff, ECC2, &
         variance , integral1,integral1b,integral2, SB(NK,NTH), integral3,&
         ajust,absajust,aa,bb,LNORM,UdotL,KdotKP,MBANDC
    REAL            :: KD, Kfactor, kscaled, kmod , CHECKSUM, ETOT
    REAL            :: SMATRIX(NTH,NTH),SMATRIX2(NTH,NTH)
    DOUBLE PRECISION :: AVECT(NTH)
 
    curtx=cx1
    curty=cy1
 
    count=0
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3SBS1')
#endif
    !
    ! 0.  Initializations ------------------------------------------------ *
    !
    !     **********************************************************
    !     ***    The initialization routine should include all   ***
    !     *** initialization, including reading data from files. ***
    !     **********************************************************
    !
    IF ( first ) THEN
      CALL insbs1( 1 )
      first  = .false.
    END IF
    IF (( (abs(cx1)+abs(cy1)).EQ.0.).AND.(matrices.EQ.0) ) THEN
      kwmin=max(max(dkbx,dkby),sig(1)**2/grav)
      kwmax=min(nkbx*dkbx,nkby*dkby)*0.25
      WRITE(*,*) 'k range:',kwmin,kwmax,sig(1)**2/grav
      CALL insbs1( 2 )
      matrices  = 1
    END IF
    !
    ! 1.  Sets scattering term to zero
    !
    d = 0.
    s = 0.
    tauscx=0.
    tauscy=0.
    !
    ! 3.  Bottom scattering ================================================== *
    !
    IF ( depth*wn(1) .LE. 6 ) THEN
      !
      ! 3.a Ardhuin and Herbers JFM 2000: no current
      !
      IF ((abs(cx1)+abs(cy1).EQ.0.).AND.(matrices.EQ.1)) THEN
        DO ik=1,nk
          kd=wn(ik)*depth
          IF ( kd .LE. 6 .AND.wn(ik).LT.kwmax ) THEN
            ! Test on kwmax means that scattering is not computed if interaction goes beyond the shortest resolved
            ! bottom component. This should probably be replaced by a warning...
            kfactor=(wn(ik)**4)*sig(ik)*pi*4.             &
                 /(sinh(2*kd)*(2*kd+sinh(2*kd)))
            kscaled=(nkscat-2)*(wn(ik)-kwmin)/(kwmax-kwmin)
            avect=dble(a(:,ik))
            IF (kscaled.LT.0) THEN
              ibk=0
              kmod=0.
            ELSE
              ibk=int(kscaled)
              kmod=mod(kscaled,1.0)
            END IF
            s((ik-1)*nth+1:ik*nth)                              &
                 =real(matmul(scatmatv(ibk,:,:),kfactor*scatmatd(ibk,:) &
                 *matmul(transpose(scatmatv(ibk,:,:)),avect))*(1.-kmod))
            s((ik-1)*nth+1:ik*nth)                              &
                 =s((ik-1)*nth+1:ik*nth)                           &
                 +real(matmul(scatmatv(ibk+1,:,:),kfactor*scatmatd(ibk+1,:) &
                 *matmul(transpose(scatmatv(ibk+1,:,:)),avect))*kmod)
            checksum=abs(sum(s((ik-1)*nth+1:ik*nth) ))
            etot=sum(a(:,ik))
            IF (checksum.GT.0.01*etot) WRITE(*,*)         &
                 'Energy not conserved:',ik,depth,checksum,etot
          ELSE
            s((ik-1)*nth+1:ik*nth)=0.
          END IF
        END DO
      ELSE
        ! 3.b
        !  Case with current (Ardhuin and Magne JFM 2007)
        ! Compute k' (WN2) from k (WN) and U (CX1, CY1)
        ! using : k'=(Cg+k.U/k)/(Cg+k'.U/k')
        !
        DO ith2=1, nth
 
          DO ispec=1, nspec
 
            ku=cx1 * ecos(mapth(ispec))+cy1 * esin(mapth(ispec))
            kpu=cx1 * ecos(ith2)+  cy1 * esin(ith2)
            cgk=cg(mapwn(ispec))
            IF ((cgk+kpu).LT.0.1*cgk) kpu=-0.9*cgk
            IF ((cgk+ku).LT.0.1*cgk)  ku=-0.9*cg(mapwn(ispec))
            wn2(ispec,ith2)= wn(mapwn(ispec))*(cgk+ku)/(cgk+kpu)
          END DO
 
        END DO
        !
        ! 3.c Compute the coupling coefficient as a product of two terms
        !
        !   K=0.5*pi k'^2 * M(k,k')^2 / [sig*sig' *(k'*Cg'+k'.U)]
        !                                      (Magne and Ardhuin JFM 2007)
        !
        !   K=Ka(k)*Kb(k,k',theta')
        !
        !   Ka = ...
        !   here Mc is neglected
        !
        DO ispec=1, nspec
          ka(ispec)= 4*pi*sig2(ispec) * wn(mapwn(ispec))  /    &
               sinh(min(2*wn(mapwn(ispec))*depth,20.))
 
          DO ith2=1, nth
            ku=cx1 * ecos(mapth(ispec))+cy1 * esin(mapth(ispec))
            kpu=cx1 * ecos(ith2)+  cy1 * esin(ith2)
            sigp=sqrt(grav*wn2(ispec,ith2)*tanh(wn2(ispec,ith2)*depth))
            cgpk=sigp*(0.5+wn2(ispec,ith2)*depth &
                 /sinh(min(2*wn2(ispec,ith2)*depth,20.)))/wn2(ispec, ith2)
 
            kb(ispec, ith2)= wn2(ispec, ith2)**3      &
                 *ec2(1+abs(mapth(ispec)-ith2)) /        &
                 (                                      &
                 2*wn2(ispec, ith2)*depth +             &
                 sinh(min(2*wn2(ispec,ith2)*depth,20.)) &
                 *(1+wn2(ispec,ith2)*kpu*2/sigp)  &
                 )
 
            !
            !  Other option for computing also Mc
            !
            !  UdotL=WN(MAPWN(ISPEC))*KU-KPU*WN2(ISPEC,ITH2)
            !  KdotKP=EC(1+ABS(MAPTH(ISPEC)-ITH2))*WN2(ISPEC,ITH2)*WN(MAPWN(ISPEC))
            !  LNORM=sqrt(WN(MAPWN(ISPEC))**2+WN2(ISPEC, ITH2)**2-2*KdotKP)
            !  MBANDC=grav*KdotKP &
            !        /(COSH(MIN(WN2(ISPEC,ITH2)*DEPTH,20.))*COSH(MIN(WN(MAPWN(ISPEC))*DEPTH,20.)
            !        +(UdotL*(SIGP*(WN(MAPWN(ISPEC))**2-KdotKP)+SIG2(ISPEC)*(KdotKP-WN2(ISPEC, ITH2)**2)) &
            !           - UdotL**2*(KdotKP-SIGP*SIG2(ISPEC)*(SIGP*SIG2(ISPEC)+UdotL**2)/GRAV**2)) &
            !           /(LNORM*(UdotL**2/(GRAV*LNORM)-TANH(MIN(LNORM*DEPTH,20.)))*COSH(MIN(LNORM*DEPTH,20.)))
            !  Kb(ISPEC,ITH2)= WN2(ISPEC, ITH2)**2
            !                   /((SIG2(ISPEC)*SIGP*WN2(ISPEC, ITH2)*(CGPK+KPU)) &
            !                  *MBANDC**2
            !
          END DO
        END DO
        !
        ! 3.a Bilinear interpolation of the bottom spectrum BOTSPEC
        !     along the locus -> B(ISPEC,ITH2)
        !
        b(:,:)=0
        DO ispec=1, nspec
          kcutoff=scattcutoff*wn(mapwn(ispec))
          DO ith2=1,nth
            kbotx=wn(mapwn(ispec))*ecos(mapth(ispec)) - &
                 wn2(ispec, ith2) * ecos(ith2)
            kboty=wn(mapwn(ispec))*esin(mapth(ispec)) - &
                 wn2(ispec, ith2) * esin(ith2)
            !
            ! 3.a.1 test if the bottom wavenumber is larger than the cutoff
            !        otherwise the interaction is set to zero
 
            IF ((kbotx**2+kboty**2)>(kcutoff**2)) THEN
 
              kbotxi=real(nkbx-mod(nkbx,2))/2.+1.+kbotx/dkbx   ! The MOD(nkbx,2) is either 1 or 0
              kbotyi=real(nkby-mod(nkby,2))/2.+1.+kboty/dkby   ! k=0 is at ik=(nkbx-1)/2+1 if kkbx is odd
 
              ibk=max(min(int(kbotxi),nkbx-1),1)
              xbk=mod(kbotxi,1.0)
              jbk=max(min(int(kbotyi),nkby-1),1)
              ybk=mod(kbotyi,1.0)
 
              b(ispec,ith2)=(                &
                   (botspec(ibk,jbk)*(1-ybk)+      &
                   botspec(ibk,jbk+1)*ybk)*(1-xbk) &
                   +               &
                   (botspec(ibk+1,jbk)*(1-ybk)+    &
                   botspec(ibk+1,jbk+1)*ybk)*xbk   &
                   )
            END IF
          END DO
        END DO
        !
        ! 4. compute Sbscat
        ! 4.a linear interpolation of A(k', theta') -> Ap
        !
 
        !  4.b computation of the source term
        integral2=0.
        integral3=0.
        smatrix(:,:)=0.
        DO ispec=1, nspec
          integral=0
          DO ith2=1, nth
            iajust=1
            DO i=2,nk
              if(wn2(ispec,ith2).GE.wn(i)) iajust=i
            END DO
            iajust=max(iajust,1)
            iajust2=min(iajust+1,nk)
            IF (iajust.EQ.iajust2) THEN
              ap=a(ith2,iajust)
            ELSE
              bb=(wn2(ispec,ith2)-wn(iajust))/(wn(iajust2)-wn(iajust))
              aa=(wn(iajust2)-wn2(ispec,ith2))/(wn(iajust2)-wn(iajust))
              ap=(a(ith2,iajust)*aa+a(ith2,iajust2)*bb)
            END IF
 
            integral=integral + ka(ispec)*kb(ispec, ith2)*b(ispec,ith2)*               &
                 ( ap*wn(mapwn(ispec))/wn2(ispec,ith2)- a(mapth(ispec),mapwn(ispec))) *dth
            ! the factor WN/WN2 accounts for the fact that N(K) and N(K')
            ! have different Jacobian transforms from kx,ky to k,theta
 
            integral1=integral1+kb(ispec, ith2)*b(ispec,ith2)*ap*wn(mapwn(ispec))/wn2(ispec,ith2)*dth &
                 *dth*dsii(mapwn(ispec))/cg(mapwn(ispec))
            integral1b=integral1b+kb(ispec, ith2)*b(ispec,ith2)*a(mapth(ispec),mapwn(ispec))*dth  &
                 *dth*dsii(mapwn(ispec))/cg(mapwn(ispec))
          END DO
          s(ispec)=s(ispec)+integral
 
          integral2=integral2+s(ispec)*dth*dsii(mapwn(ispec))/cg(mapwn(ispec))
          integral3=integral3+abs(s(ispec))*dth*dsii(mapwn(ispec))/cg(mapwn(ispec))
        END DO
      END IF
#ifdef W3_T
      print*,'BOTTOM SCAT CHECKSUM:',integral2,integral3,integral1,integral1b
#endif
 
#ifdef W3_T
      DO ith=1,120
        WRITE(6,'(120G15.7)') smatrix(ith,:)
      END DO
#endif
    END IF
 
    !/
    !/ End of W3SBS1 ----------------------------------------------------- /
    !/
  END SUBROUTINE w3sbs1
  !/ ------------------------------------------------------------------- /
 
  SUBROUTINE insbs1( inistep )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |            F. Ardhuin             |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         23-Jun-2006 |
    !/                  +-----------------------------------+
    !/
    !/    23-Jun-2006 : Origination.                        ( version 3.09 )
    !/
    !  1. Purpose :
    !
    !     Initialization for bottom scattering source term routine.
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      STRACE    Subr. W3SERVMD Subroutine tracing.
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3SBS1    Subr. W3SBS1MD Corresponding source term.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !       None.
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/S  Enable subroutine tracing.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE w3gdatmd, ONLY: nk, nth, nspec, sig, dth, dden, ecos, esin
    USE w3servmd, ONLY: diagonalize
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN)        :: inistep
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    INTEGER         :: I, J, K1, K2, IK, JK, NROT
    REAL            :: kbotx, kboty, kcurr, kcutoff, variance
    REAL            :: kbotxi, kbotyi, xk, yk
    DOUBLE PRECISION, ALLOCATABLE,DIMENSION(:,:) :: AMAT, V
    DOUBLE PRECISION, ALLOCATABLE,DIMENSION(:)   :: D
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'INSBS1')
#endif
    !
    IF (inistep.EQ.1) THEN
      !
      ! 1.  Reads bottom spectrum
      !
      OPEN(183,file= 'bottomspectrum.inp', status='old')
      READ(183,*) nkbx, nkby
      READ(183,*) dkbx, dkby
      WRITE(*,*) 'Bottom spec. dim.:', nkbx, nkby, dkbx, dkby
      ALLOCATE(botspec(nkbx, nkby))
      DO i=1, nkbx
        READ(183,*) botspec(i,:)
      END DO
      CLOSE(183)
      variance=0
      DO i=1,nkbx
        DO j=1,nkby
          variance=variance+botspec(i,j)*dkbx*dkby
        END DO
      END DO
      WRITE(*,*) 'Bottom variance:', variance
      !
    ELSE
      !
      ! 2.  Precomputed the scatering matrices for zero current
      !
      ! The Scattering source term is expressed as a matrix problem for
      ! a list of wavenumbers k0
      ! in the range of wavenumbers used in the model.
      ! i.e. S(k0,theta)=Kfactor*SCATMATA ** TRANSPOSE (E(k0,theta))
      !
      ! in which
      !
      ! Kfactor is a scalar computed in CALCSOURCE as
      !        Kfactor=tailfactor*(Kp(I,J)**4)*2.*pi*FREQ(J)*pi*4./(SINH(HND)*(HND+SINH(HND)))
      !
      ! SCATMATA is a square matrix of size NTH*NTH
      !
      ! S(k0,theta) and E(k0,theta) are the vectors giving the directional source term
      !   and spectrum at a fixed wavenumber
      !
      ALLOCATE(scatmata(0:nkscat-1,1:nth,1:nth))
      ALLOCATE(amat(nth,nth))
      DO i=0,nkscat-1
        ! kcurr is the current surface wavenumber for which
        ! the scattering matrices are evaluated
        kcurr=kwmin+i*(kwmax-kwmin)/(nkscat-2)
        kcutoff=scattcutoff*kcurr
        DO k1=1,nth
          DO k2=1,nth
            kbotx=-kcurr*(ecos(k2)-ecos(k1))
            kboty=-kcurr*(esin(k2)-esin(k1))
            amat(k1,k2)=0.
            ! Tests if the bottom wavenumber is larger than the cutoff
            ! Otherwise the interaction is set to zero
            IF ((kbotx**2+kboty**2) > (kcutoff**2)) THEN
              !WARNING : THERE MAY BE A BUG : spectrum not symmetric when
              ! nkbx is odd !!
 
              kbotxi=real(nkbx)/2.+1.+kbotx/dkbx
              kbotyi=real(nkby)/2.+1.+kboty/dkby
              !WRITE(6,*) 'Bottom wavenumber i:',kbotxi,kbotyi
              ik=int(kbotxi)
              xk=mod(kbotxi,1.0)
              jk=int(kbotyi)
              yk=mod(kbotyi,1.0)
              IF (ik.GE.nkbx) ik=nkbx-1
              IF (jk.GE.nkby) jk=nkby-1
              IF (ik.LT.1) ik=1
              IF (jk.LT.1) jk=1
              ! Bilinear interpolation of the bottom spectrum
              amat(k1,k2)=((botspec(ik,jk  )  *(1-yk)           &
                   +botspec(ik,jk+1)  *yk    )*(1-xk)   &
                   +(botspec(ik+1,jk)  *(1-yk)           &
                   +botspec(ik+1,jk+1)*yk)    *xk)      &
                   *(ecos(k1)*ecos(k2)+esin(k1)*esin(k2))**2
            END IF
          END DO
          amat(k1,k1)=amat(k1,k1)-sum(amat(k1,:))
        END DO
        amat(:,:)=dth*(amat(:,:)+transpose(amat(:,:)))*0.5
        !makes sure the matrix is exactly symmetric
        !which should already be the case if the bottom
        ! spectrum is really symmetric
        scatmata(i,:,:)=amat(:,:)
      END DO
      ALLOCATE(scatmatd(0:nkscat-1,nth))
      ALLOCATE(scatmatv(0:nkscat-1,nth,nth))
      ALLOCATE(v(nth,nth))
      ALLOCATE(d(nth))
      DO i=0,nkscat-1
        amat(:,:)=scatmata(i,:,:)
        !
        !diagonalizes the matrix A
        !D is a vector with the eigenvalues, V is the matrix made of the
        !eigenvectors so that VD2Vt=A  with D2(i,j)=delta(i,j)D(i)
        !and VVt=Id, so that exp(A)=Vexp(D2)Vt
        !
        CALL diagonalize(amat,d,v,nrot)
        scatmatd(i,:)=d(:)    !eigen values
        scatmatv(i,:,:)=v(:,:) !eigen vectors
        kcurr=kwmin+i*(kwmax-kwmin)/(nkscat-2)
        WRITE(*,*) 'Scattering matrix diagonalized for k=  ',kcurr,',',i+1,'out of ',nkscat
      END DO
    END IF
 
    !/
    !/ End of INSBS1 ----------------------------------------------------- /
    !/
  END SUBROUTINE insbs1
  !/
  !/ End of module INSBS1MD -------------------------------------------- /
  !/
END MODULE w3sbs1md