twcal2.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE twcal2
!                   *****************
!
!***********************************************************************
! ARTEMIS   V7P4                                     Nov 2017
!***********************************************************************
!
!brief    DISCRETISES AN ENERGY SPECTRUM IN NDALE x NPALE CELLS
!+                OF EQUAL ENERGY. THE RESULT IS A LIST OF
!+                DIRECTIONS CORRESPONDING TO EACH BAND AND
!+                FOR EACH LIQUID BOUNDARY NODE.
!
!history  N.DURAND (HRW)
!+        August 2017
!+        V7P3
!+   New. Computes BDALEs: spatially varying directional components
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| BDALE          |<--| ISO ENERGY DIRECTIONS
!|                |   |     FOR EACH LIQUID BOUNDARY NODE
!| PMIN           |-->| MINIMUM FREQUENCY FOR SPECTRUM
!| PMAX           |-->| MAXIMUM FREQUENCY FOR SPECTRUM
!| TETMIN         |-->| MAXIMUM VALUE FOR THE PROPAGATION ANGLE
!| TETMAX         |-->| MAXIMUM VALUE FOR THE PROPAGATION ANGLE
!| NPALE          |-->| NUMBER OF DISCRETISATION BAND FOR PERIODS
!| NDALE          |-->| NUMBER OF DISCRETISATION BAND FOR DIRECTIONS
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief
      USE declarations_artemis
      USE declarations_telemac, ONLY : kinc,ksort
      USE interface_artemis, ONLY : stwc1
      USE interface_parallel, ONLY : p_max
!
      USE declarations_special
      IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
      INTEGER                       :: IERR
      INTEGER                       :: IPTFR
      INTEGER                       :: JTMP,IMIN,IMAX
      INTEGER                       :: NPASF,NPASD
      INTEGER                       :: IDALE,IDD,IFF,I
      INTEGER                       :: ISPEC
!
      DOUBLE PRECISION              :: FMIN,FMAX
      DOUBLE PRECISION              :: DTETA,DF,DTETA2
      DOUBLE PRECISION              :: SUMB,POIDS,SUMD,VAR,SUMICI
!
      INTEGER,ALLOCATABLE           :: IDTWC(:)
!
      DOUBLE PRECISION,ALLOCATABLE  :: SDIRTWC(:,:)
      DOUBLE PRECISION,ALLOCATABLE  :: SDIR(:)
!
      INTEGER                       :: IDDMINART
      DOUBLE PRECISION,ALLOCATABLE  :: DIR_ART(:)
      DOUBLE PRECISION,ALLOCATABLE  :: SDIR_ART(:)
!
      INTEGER                       :: PRINTMSG
!
!=======================================================================
!
!     MIN/MAX FREQUENCY
      fmin = 1.d0 / pmax
      fmax = 1.d0 / pmin
!
!     NUMBER OF INTEGRATION INTERVALS FOR THE TRAPEZOIDS METHOD : DIRECTIONS
      npasd = 500*ndale
!
!     NUMBER OF INTEGRATION INTERVALS FOR THE TRAPEZOIDS METHOD : FREQUENCIES
      npasf = 500*npale
!
!     WIDTH OF AN INTEGRATION INTERVAL : DIRECTIONS
      dteta = (tetmax-tetmin)/float(npasd)
!
!     WIDTH OF AN INTEGRATION INTERVAL : FREQUENCIES
      df = (fmax-fmin)/float(npasf)
!
!=======================================================================
!
      ALLOCATE(idtwc(2*ndale+2),stat=ierr)
      CALL check_allocate(ierr,'TWCAL2:IDTWC')
      ALLOCATE(sdirtwc(nspec,ndir+1),stat=ierr)
      CALL check_allocate(ierr,'TWCAL2:SDIRTWC')
      ALLOCATE(sdir(npasd+1),stat=ierr)
      CALL check_allocate(ierr,'TWCAL2:SDIR')
      ALLOCATE(dir_art(npasd+1),stat=ierr)
      CALL check_allocate(ierr,'TWCAL2:DIR_ART')
      ALLOCATE(sdir_art(npasd+1),stat=ierr)
      CALL check_allocate(ierr,'TWCAL2:SDIR_ART')
!
      printmsg = 0
!
!-----------------------------------------------------------------------
!
!     INTEGRAL OF THE TOMAWAC SPECTRUM (TRAPEZOIDS METHOD)
!     FOR EACH DIRECTION COMPONENT
!
      DO ispec = 1,nspec
!
        DO idd = 1,ndir+1
          sumd=0.d0
          DO iff = 1,npasf+1
!           IF FREQ AND/OR DIR ON THE BOUNDARY OF THE INTEGRATION DOMAIN :  CONTRIBUTION/2
            poids=1.d0
            IF(iff.EQ.1 .OR. iff.EQ.npasf+1) poids=0.5d0*poids
            IF(idd.EQ.1 .OR. idd.EQ.ndir+1)  poids=0.5d0*poids
            var=stwc1(fmin+float(iff-1)*df,s_tom%DIR(idd),s_tom,ispec)
            sumd = sumd + poids*var*df
          ENDDO
          sdirtwc(ispec,idd) = sumd
!         IF(DEBUG.GT.0) WRITE(LU,*) ISPEC,S_TOM%DIR(IDD),
!    &                               SDIRTWC(ISPEC,IDD)
        ENDDO
!
      ENDDO ! ISPEC = 1,NSPEC
      IF(debug.GT.0)
     &    WRITE(lu,*) '< TWCAL2: SDIRTWC COMPUTED FOR 1 TO ',nspec
!
!-----------------------------------------------------------------------
!
!     FOR THE DIRECTION RANGE SPECIFIED IN ARTEMIS
!
!     INDICES COVERING THE RANGE
      IF(tetmax-tetmin.EQ.360.d0) THEN
        imin = 1
      ELSEIF(tetmin.LE.0) THEN
        DO jtmp = 1,ndir+1
          IF(s_tom%DIR(jtmp).LE.tetmin+360.d0) imin = jtmp
        ENDDO
      ELSE
        DO jtmp = 1,ndir+1
          IF(s_tom%DIR(jtmp).LE.tetmin) imin = jtmp
        ENDDO
      ENDIF
!
      IF(tetmax-tetmin.EQ.360.d0) THEN
        imax = ndir+1
      ELSEIF(tetmax.LT.0) THEN
        DO jtmp = ndir+1,1,-1
          IF(s_tom%DIR(jtmp).GE.tetmax+360.d0) imax = jtmp
        ENDDO
      ELSE
        DO jtmp = ndir+1,1,-1
          IF(s_tom%DIR(jtmp).GE.tetmax) imax = jtmp
        ENDDO
      ENDIF
!
      IF(debug.GT.0) WRITE(lu,*)'< TWCAL2: TETMIN,IMIN,S_TOM%DIR(IMIN)'
      IF(debug.GT.0) WRITE(lu,*) tetmin,imin,s_tom%DIR(imin)
      IF(debug.GT.0) WRITE(lu,*)'< TWCAL2: TETMAX,IMAX,S_TOM%DIR(IMAX)'
      IF(debug.GT.0) WRITE(lu,*) tetmax,imax,s_tom%DIR(imax)
!
!=======================================================================
!
!     FOR EACH ARTEMIS BOUNDARY NODE
!
      DO iptfr = 1 , nptfr
!
        IF (lihbor%I(iptfr).EQ.kinc .OR. lihbor%I(iptfr).EQ.ksort) THEN
!
          printmsg = printmsg + 1
!
!-----------------------------------------------------------------------
!
          IF(debug.GT.0 .AND. printmsg.EQ.1) THEN
            WRITE(lu,*)
            WRITE(lu,*) '> TWCAL2: ',
     &      'STARTS INTERPOLATION OF DIR. DISTRIBUTION TO ARTEMIS NODES'
          ENDIF
!
!     INTERPOLATES THE SPECTRUM FOR INCIDENT BOUNDARY POINTS, ONLY
!
          DO idd = 1,npasd+1
            sdir(idd) = 0.d0
          ENDDO
!
          CALL fasp_sp( s_tom%XOUTER,s_tom%YOUTER,sdirtwc,nspec,
     &         x(mesh%NBOR%I(iptfr)),y(mesh%NBOR%I(iptfr)),sdir,iptfr )
!
          CALL stwc2( imin,imax,npasd+1,dir_art,sdir )
          dteta2 = dir_art(2)-dir_art(1)
!
          IF(debug.GT.0 .AND. printmsg.EQ.1)
     &     WRITE(lu,*) '< TWCAL2: INTERPOLATION TO ARTEMIS NODES ENDED'
!
!-----------------------------------------------------------------------
!
!     REORGANISE THE ARTEMIS SPECTRUM PRIOR TO SLICING
!
          iddminart = 1
          DO idd = 2,npasd
            IF(sdir(idd).LT.sdir(iddminart)) iddminart = idd
          ENDDO
!
!         REORGANISED SPECTRUM IN SDIR_ART(NPASD+1)
          DO idd = 1,npasd
            i=idd-iddminart+1
            IF (idd.LT.iddminart) i=i+npasd
            sdir_art(i)=sdir(idd)
          ENDDO
          sdir_art(npasd+1)=sdir_art(1)
!
!     INTEGRAL OF THE SPECTRUM (TRAPEZOIDS METHOD)
!
          sumb = 0.d0
          DO idd = 1,npasd+1
            poids=1.d0
            IF(idd.EQ.1 .OR. idd.EQ.npasd+1)  poids=0.5d0*poids
            sumb = sumb + poids*sdir_art(idd)*dteta2*degrad
          ENDDO
!
!     SIGNIFICANT WAVE HEIGHT CORRESPONDING TO TOTAL ENERGY STORAGE
!
          hscal=4.d0*sqrt(sumb)
!
          IF(debug.GT.0 .AND. printmsg.EQ.1) THEN
            WRITE(lu,*) '> TWCAL2: STARTS SPECTRUM DISCRETISATION'
            WRITE(lu,*) '           DIRECTIONS'
          ENDIF
!
!         =======================================================
!                              DIRECTIONS
!         =======================================================
!         DIVIDES THE SPECTRUM INTO 2*NDALE BANDS OF EQUAL ENERGY
          sumb = sumb/float(2*ndale)
!
!         FIRST TERM OF DIRECTION DISCRETIZATION
          i=1
          idtwc(i) = 1
!
!         IDENTIFIES THE ANGLES EVERY (I)*SUMB (I=1,NDALE)
          sumici = 0.d0
!
          DO idd = 1,npasd+1
            sumici = sumici + sdir_art(idd)*dteta2*degrad
!
!           CHECKS IF SUMB = OVERALL SUMB/(2 NDALE) IS REACHED AND SAVES TETA
!
            IF (sumici.GE.sumb*float(i) .OR. idd.EQ.npasd+1) THEN
              i=i+1
              idtwc(i) = idd
            ENDIF
          ENDDO
!
!         IDTWC HOLDS INDICES FOR THE DIRECTIONS TO BE CONSIDERED IN ARTEMIS
!
!         THESE ARE BASED ON CONTINUOUS SPECTRUM AND MAY NEED ADJUSTING,
!         IF PRINTED OUT FOR CHECKS, AS DONE FOR BDALE BELOW
!
          DO i=1,2*ndale-1,2
            idale = (i+1)/2
            bdale%ADR(idale)%P%R(iptfr)=
     &               float(idtwc(i+1)-1)*dteta2+dir_art(iddminart)
            IF(bdale%ADR(idale)%P%R(iptfr).GT.dir_art(npasd+1))
     &        bdale%ADR(idale)%P%R(iptfr)=
     &        bdale%ADR(idale)%P%R(iptfr)-dir_art(npasd+1)+dir_art(1)
!
            IF(bdale%ADR(idale)%P%R(iptfr).GT.180d0)
     &        bdale%ADR(idale)%P%R(iptfr)=
     &        bdale%ADR(idale)%P%R(iptfr)-360.d0
          ENDDO
!
        ENDIF ! LIHBOR%I(IPTFR).EQ.KINC
!
!-----------------------------------------------------------------------
!
      IF(debug.GT.0 .AND. printmsg.EQ.1)
     &  WRITE(lu,*) '< TWCAL2: SPECTRUM DISCRETISATION ENDED'
!
      ENDDO ! IPTFR = 1 , NPTFR
!
      DEALLOCATE(idtwc,sdirtwc,sdir)
      DEALLOCATE(dir_art,sdir_art)
!
!-----------------------------------------------------------------------
!
!     FINDS BOUNDARY NODE CLOSEST TO REFERENCE POINT
!    (USED FOR PRINTOUTS ONLY IN ARTEMIS.F)
!
      iptfr_ref = 0
      IF (nptfr.GT.0) CALL twcclosest
      IF (ncsize.GT.1) iptfr_ref = p_max(iptfr_ref)
      IF(debug.GT.0)
     &      WRITE(lu,*) 'CLOSEST: REFERENCE BOUNDARY POINT:',iptfr_ref
!
!-----------------------------------------------------------------------
!
      RETURN
      END