This module computes a scattering term based on the theory by Ardhuin and Magne (JFM 2007). More...

Functions/Subroutines
subroutine	w3sbs1 (A, CG, WN, DEPTH, CX1, CY1, TAUSCX, TAUSCY, S, D)
	Bottom scattering source term. More...

subroutine	insbs1 (inistep)
	Initialization for bottom scattering source term routine. More...

Variables
real, dimension(:,:), allocatable	botspec

integer, parameter	nkscat = 30

double precision, dimension(:,:,:), allocatable	scatmatv

double precision, dimension(:,:,:), allocatable	scatmata

double precision, dimension(:,:), allocatable	scatmatd

character(len=10)	botspec_indicator

integer	nkbx

integer	nkby

real	dkbx

real	dkby

real	kwmin

real	kwmax

real, parameter	scattcutoff =0.

real	curtx

real	curty

Detailed Description

This module computes a scattering term based on the theory by Ardhuin and Magne (JFM 2007).

Author: F. Ardhuin

Date: 14-Nov-2010

Copyright: Copyright 2009-2022 National Weather Service (NWS), National Oceanic and Atmospheric Administration. All rights reserved. WAVEWATCH III is a trademark of the NWS. No unauthorized use without permission.

Function/Subroutine Documentation

◆ insbs1()

subroutine w3sbs1md::insbs1 ( integer, intent(in) inistep )

Initialization for bottom scattering source term routine.

Parameters

[in] inistep

Author: F. Ardhuin

Date: 23-Jun-2006

Definition at line 474 of file w3sbs1md.F90.

     !/
     !/                  +-----------------------------------+
     !/                  | 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 ----------------------------------------------------- /
     !/

References botspec, w3gdatmd::dden, w3servmd::diagonalize(), dkbx, dkby, w3gdatmd::dth, w3gdatmd::ecos, w3gdatmd::esin, file(), kwmax, kwmin, w3gdatmd::nk, nkbx, nkby, nkscat, w3gdatmd::nspec, w3gdatmd::nth, scatmata, scatmatd, scatmatv, scattcutoff, w3gdatmd::sig, and w3servmd::strace().

Referenced by w3sbs1().

◆ w3sbs1()

subroutine w3sbs1md::w3sbs1	(	real, dimension(nth,nk), intent(in)	A,
		real, dimension(nk), intent(in)	CG,
		real, dimension(nk), intent(in)	WN,
		real, intent(in)	DEPTH,
		real, intent(in)	CX1,
		real, intent(in)	CY1,
		real, intent(out)	TAUSCX,
		real, intent(out)	TAUSCY,
		real, dimension(nspec), intent(out)	S,
		real, dimension(nspec), intent(out)	D
	)

Bottom scattering source term.

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.

Parameters

[in]	A	Action density spectrum (1-D)
[in]	CG	Group velocities
[in]	WN	Wavenumbers
[in]	DEPTH	Mean water depth
[in]	CX1	Current components at ISEA
[in]	CY1	Current components at ISEA
[out]	TAUSCX	Change of wave momentum due to scattering
[out]	TAUSCY	Change of wave momentum due to scattering
[out]	S	Source term (1-D version)
[out]	D	Diagonal term of derivative (1-D version)

Author: F. Ardhuin

Date: 23-Jun-2006

Definition at line 114 of file w3sbs1md.F90.

     !/
     !/                  +-----------------------------------+
     !/                  | 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 ----------------------------------------------------- /
     !/

References botspec, curtx, curty, w3gdatmd::dden, dkbx, dkby, w3gdatmd::dsii, w3gdatmd::dth, w3gdatmd::ec2, w3gdatmd::ecos, w3gdatmd::esin, constants::grav, insbs1(), kwmax, kwmin, w3gdatmd::mapth, w3gdatmd::mapwn, w3gdatmd::nk, nkbx, nkby, nkscat, w3gdatmd::nspec, w3gdatmd::nth, constants::pi, scatmatd, scatmatv, scattcutoff, w3gdatmd::sig, w3gdatmd::sig2, and w3servmd::strace().

Referenced by gxexpo(), w3exnc(), w3expo(), and w3srcemd::w3srce().

Variable Documentation

◆ botspec

real, dimension(:,:), allocatable w3sbs1md::botspec

Definition at line 79 of file w3sbs1md.F90.

79 REAL, DIMENSION(:,:), ALLOCATABLE :: BOTSPEC

Referenced by insbs1(), and w3sbs1().

◆ botspec_indicator

character(len=10) w3sbs1md::botspec_indicator

Definition at line 84 of file w3sbs1md.F90.

84 CHARACTER(len=10) :: botspec_indicator

◆ curtx

real w3sbs1md::curtx

Definition at line 88 of file w3sbs1md.F90.

88 REAL :: CURTX, CURTY

Referenced by w3sbs1().

◆ curty

real w3sbs1md::curty

Definition at line 88 of file w3sbs1md.F90.

Referenced by w3sbs1().

◆ dkbx

real w3sbs1md::dkbx

Definition at line 86 of file w3sbs1md.F90.

86 REAL :: dkbx, dkby, kwmin, kwmax

Referenced by insbs1(), and w3sbs1().

◆ dkby

real w3sbs1md::dkby

Definition at line 86 of file w3sbs1md.F90.

Referenced by insbs1(), and w3sbs1().

◆ kwmax

real w3sbs1md::kwmax

Definition at line 86 of file w3sbs1md.F90.

Referenced by insbs1(), and w3sbs1().

◆ kwmin

real w3sbs1md::kwmin

Definition at line 86 of file w3sbs1md.F90.

Referenced by insbs1(), and w3sbs1().

◆ nkbx

integer w3sbs1md::nkbx

Definition at line 85 of file w3sbs1md.F90.

85 INTEGER :: nkbx, nkby

Referenced by insbs1(), and w3sbs1().

◆ nkby

integer w3sbs1md::nkby

Definition at line 85 of file w3sbs1md.F90.

Referenced by insbs1(), and w3sbs1().

◆ nkscat

integer, parameter w3sbs1md::nkscat = 30

Definition at line 80 of file w3sbs1md.F90.

80 INTEGER, PARAMETER :: NKSCAT = 30 !number of wavenumbers

Referenced by insbs1(), and w3sbs1().

◆ scatmata

double precision, dimension(:,:,:), allocatable w3sbs1md::scatmata

Definition at line 82 of file w3sbs1md.F90.

82 DOUBLE PRECISION ,DIMENSION(:,:,:) , ALLOCATABLE :: SCATMATA !original matrix

Referenced by insbs1().

◆ scatmatd

double precision, dimension(:,:), allocatable w3sbs1md::scatmatd

Definition at line 83 of file w3sbs1md.F90.

83 DOUBLE PRECISION ,DIMENSION(:,:) , ALLOCATABLE :: SCATMATD

Referenced by insbs1(), and w3sbs1().

◆ scatmatv

double precision, dimension(:,:,:), allocatable w3sbs1md::scatmatv

Definition at line 81 of file w3sbs1md.F90.

81 DOUBLE PRECISION ,DIMENSION(:,:,:) , ALLOCATABLE :: SCATMATV !scattering matrices

Referenced by insbs1(), and w3sbs1().

◆ scattcutoff

real, parameter w3sbs1md::scattcutoff =0.

Definition at line 87 of file w3sbs1md.F90.

87 REAL, PARAMETER :: scattcutoff=0.

Referenced by insbs1(), and w3sbs1().

Functions/Subroutines

Variables

Detailed Description

Function/Subroutine Documentation

◆ insbs1()

◆ w3sbs1()

Variable Documentation

◆ botspec

◆ botspec_indicator

◆ curtx

◆ curty

◆ dkbx

◆ dkby

◆ kwmax

◆ kwmin

◆ nkbx

◆ nkby

◆ nkscat

◆ scatmata

◆ scatmatd

◆ scatmatv

◆ scattcutoff