w3gig1md.F90

Go to the documentation of this file.

#include "w3macros.h"
!/ ------------------------------------------------------------------- /
MODULE w3gig1md
  !/
  !/                  +-----------------------------------+
  !/                  | WAVEWATCH III                     |
  !/                  |     A. Rawat and  F. Ardhuin      |
  !/                  |                        FORTRAN 90 |
  !/                  | Last update :         05-Jul-2012 |
  !/                  +-----------------------------------+
  !/
  !/    31-Mar-2010 : Origination.                        ( version 4.07 )
  !/
  !  1. Purpose :
  !
  !     This module computes :
  !        - the second order spectrum, in particular for infragravity waves
  !
  !  2. Variables and types :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !     ----------------------------------------------------------------
  !
  !  3. Subroutines and functions :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !      W3SREF    Subr. Public   Reflection of waves (shorline, islands...)
  !     ----------------------------------------------------------------
  !
  !  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
  !/
  !
  !/
CONTAINS
  !/ ------------------------------------------------------------------- /
 
  FUNCTION df1f2theta(s1,s2,WN1,WN2,theta,DEPTH)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         29-Nov-1999 |
    !/                  +-----------------------------------+
    !/                                Based on INCYMD of the GLA GCM.
    !/
    !/    18-Oct-1998 : Final FORTRAN 77                    ( version 1.18 )
    !/    29-Nov-1999 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/
    !  1. Purpose :
    !
    !      Computes the coupling coefficient between waves of frequencies f1 and f2
    !      and an angle theta.
    !       This is for the surface elevation variance
    !      See Okihiro et al. 1992
    !       Code adapted from Matlab by Arshad Rawat, 2012.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       NYMD    Int.   I   Old date in YYMMDD format.
    !       M       Int.   I   +/- 1 (Day adjustment)
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      STRACE    Subr. W3SERVMD Subroutine tracing.
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !     Any subroutine.
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/S  Enable subroutine tracing using STRACE.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    USE constants
 
    IMPLICIT NONE
 
    REAL, INTENT(IN)        :: s1,s2,theta,depth
    REAL                    :: df1f2theta,wn1,wn2
    REAL                    :: k1,k2,co,cok1,cok2,k3,c1,c2,c3,c4
    REAL                    :: c1b,s3,sk2,g2,g
 
    k1=wn1
    k2=wn2
    co=cos(theta)
    g2=grav**2
    s3=s1+s2
    k3=sqrt(k1**2+k2**2+2*k1*k2*co)
    g=grav
    sk2=g*k3*tanh(k3*depth)
 
    c1=-(k1*k2*co)/(s1*s2)
    c1b=(s3**2-s1*s2)/g2
    c2=s3
    c3=(s3**2-sk2)*s1*s2
 
    ! C4 is Hasselmann's D times i
 
    c4=s3*(k1*k2*co-((s1*s2)**2)/g2)+0.5*(s1*k2**2+s2*k1**2-s1*s2*(s2**3+s1**3)/g2)
 
    df1f2theta=g*(0.5*(c1+c1b)+(c2*c4/c3));
 
    RETURN
  END FUNCTION df1f2theta
 
 
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3addig(E,DEPTH,WN,CG,IACTION)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III                     |
    !/                  |     A. Rawat and  F. Ardhuin      |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         05-Jul-2012 |
    !/                  +-----------------------------------+
    !/
    !/    31-Mar-2010 : Origination.                        ( version 4.07 )
    !/
    !  1. Purpose :
    !
    !     This subroutine computes :
    !        - the second order spectrum, in particular for infragravity waves
    !  2. Method :
    !     Uses 2nd order coupling coefficient (Biesel 1952, Hasselmann 1962)
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       E         R.A. I/O   Energy density spectrum (1-D), f-theta
    !       DEPTH     Real I     Water depth
    !       WN        R.A.       wavenumbers
    !       CG        R.A.       group velocities
    !       IACTION   Int  I     Switch to specify if the input spectrum
    !                            is E(f,theta) or A(k,theta)
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      STRACE    Subr. W3SERVMD Subroutine tracing.
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3SREF    Subr. W3REF1MD Shoreline reflection source term
    !      W3EXPO    Subr.   N/A    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 w3dispmd
    USE w3gdatmd, ONLY: nk, nth, nspec, sig, th, dth, dden,  &
         ecos, esin, ec2, mapth, mapwn, &
         dsip, iobpd, gtype, ungtype, igpars
 
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    !/
    !
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    REAL, INTENT(INOUT)     :: E(NSPEC)
    REAL, INTENT(IN)        :: DEPTH
    REAL, INTENT(IN)        :: WN(NK)
    REAL, INTENT(IN)        :: CG(NK)
    INTEGER, INTENT(IN)     :: IACTION
 
    !*****************************************************************************
    ! Computes the "second order spectrum" (only difference interaction, not sum)
    !*****************************************************************************
    ! Reads in the wave frequency-directional spectrum
    !
 
 
    INTEGER     :: NKIG,iloc,NSPECIG
    INTEGER     :: i,iIG,IFR,IK,ith,ith1,ith2,itime,I2, ISP1, ISP2, ISP3
    INTEGER , DIMENSION(:,:),      ALLOCATABLE :: ifr2c
 
    REAL :: d,deltaf,dfIG,CG2
    REAL :: WN1,K1,K2,Dkx,Dky,Eadd,thetaIG,memo
 
    REAL   , DIMENSION(:),        ALLOCATABLE :: df,fIG,II,Efmall
    REAL   , DIMENSION(:,:),      ALLOCATABLE :: wfr1,Efth
    REAL   , DIMENSION(:),        ALLOCATABLE :: EfthIG
    REAL   , DIMENSION(:,:,:,:),  ALLOCATABLE :: DD
    REAL   , DIMENSION(NSPEC)                 :: ESPEC
    CHARACTER(120) ::path,filename,filename2
 
 
    ! Defines the spectral domain for the IG computation
    nkig=igpars(5)
    nspecig=nkig*nth
 
    ALLOCATE(dd(nkig,nk,nth,nth))
    ALLOCATE(wfr1(nkig,nk))
    ALLOCATE(ifr2c(nkig,nk))
    ALLOCATE(efthig(nspecig))
    efthig(:)=0.
 
    !  WRITE(*,*) 'Computing coupling coefficient for SURFACE ELEVATION'
 
    IF (iaction.EQ.0) THEN
      espec=e
    ELSE
      DO ik = 1,nk
        DO ith = 1, nth
          isp1=ith+(ik-1)*nth
          espec(isp1)=e(isp1)*sig(ik)*tpi / cg(ik)
        END DO
      END DO
    END IF
    !
    DO iig=1,nkig
      DO ifr=1,nk
        CALL wavnu1 (sig(ifr)+sig(iig),depth,wn1,cg2)
        DO ith1=1,nth
          DO ith2=1,nth
            !
            ! This is the coupling coefficient for the SURFACE ELEVATION. See .e.g. forristall (JPO 2000)
            !
            dd(iig,ifr,ith1,ith2)=(df1f2theta(sig(ifr)+sig(iig),-sig(ifr), wn1,wn(ifr), &
                 (abs(th(ith1)-th(ith2))+pi),depth))**2
 
          END DO
        END DO
        !
        ! weights
        !
        wfr1(iig,ifr)=dble(dsip(ifr))*dth
        !
        ! Computes indices for a proper integration over the spectral domain using Rectangle's rule
        ! since we integrate E(f)*E(f+fIG)*df  for a fixed fIG
 
        iloc=1
 
        if (sig(iig) < 0.5*dsip(ifr))THEN
          ifr2c(iig,ifr)=ifr
        else
          iloc=minloc(abs((sig(1:nk)-dsip(1:nk))-(sig(iig)+sig(ifr))), 1)
          !find(f-df< (fIG(iIG)+f(ifr)))
          if (iloc /= 0) THEN
            ifr2c(iig,ifr)=iloc  ! index of frequency f+fIG
          else
            ifr2c(iig,ifr)=nk
          end if
 
          !wfr1(iIG,ifr)=0.0
        end if
      end do
    end do
 
 
    DO iig=1,nkig
      DO ifr = 1,nk-1
 
        ! AR calculating k1 and k2 before loops on th1 and th2
 
        k1=wn(ifr)
        k2=wn(ifr2c(iig,ifr))
 
        DO ith1 = 1,nth
          DO ith2 = 1,nth
 
            ! Adds the effect of interaction of frequency f(ifr), theta(ith1) with f(ifr)+fIG(:), theta(ith2)
            isp1 = ith1 + (ifr2c(iig,ifr)-1)*nth
            isp2 = ith2 + (ifr-1)*nth
 
            eadd=dd(iig,ifr,ith1,ith2)*wfr1(iig,ifr) &
                 *espec(isp1)*espec(isp2) ! Rectangle rule by AR
            dkx=k2*cos(dble(dth*ith2))- k1*cos(dble(dth*ith1))
            dky=k2*sin(dble(dth*ith2))- k1*sin(dble(dth*ith1))
 
            thetaig=atan2(dky,dkx)
 
            if (thetaig.LT.0) thetaig=2*pi+thetaig
            ! Finding corresponding index of theta IG in theta array
            !I=INT((thetaIG/(2*pi))*nth)
 
            i=minloc(abs(thetaig-th), 1)-1
            if (i==0) i=nth
            isp3 = i + (iig-1)*nth
            !            memo=EfthIG(ISP3)
            efthig(isp3)= efthig(isp3)+eadd;
            !              IF (EfthIG(ISP3).NE.EfthIG(ISP3).AND.Eadd.NE.0) WRITE(6,*) 'EfthIG:',IIG, IFR, ITH1,ITH2,ISP3, &
            !                                                                        EfthIG(ISP3),Eadd,memo
          END DO
        END DO
      end do
    end do
 
    !   ESPEC(1:NSPECIG)=ESPEC(1:NSPECIG)+EfthIG(:)
    espec(1:nspecig)=efthig(:)
 
    IF (iaction.EQ.0) THEN
      DO isp1=1,nspecig
        e(isp1)=espec(isp1)
      END DO
    ELSE
      DO ik = 1,nkig
        DO ith = 1, nth
          isp1=ith+(ik-1)*nth
          e(isp1)=espec(isp1)*cg(ik)/(sig(ik)*tpi)
        END DO
      END DO
    END IF
 
    ! OPEN(5555,FILE='testos.dat',status='unknown')
    ! WRITE(5555,*) E,EfthIG !f,fIG,tet!ifr2c !Efth, !!, Efth,
 
    !/
    !/ End of W3ADDIG ----------------------------------------------------- /
    !/
  END SUBROUTINE w3addig
  !/ ------------------------------------------------------------------- /
 
  !/
  !/ End of module W3GIG1MD -------------------------------------------- /
  !/
END MODULE w3gig1md