The TELEMAC-MASCARET system  trunk
twcal2.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE twcal2
3 ! *****************
4 !
5 !***********************************************************************
6 ! ARTEMIS V7P4 Nov 2017
7 !***********************************************************************
8 !
9 !brief DISCRETISES AN ENERGY SPECTRUM IN NDALE x NPALE CELLS
10 !+ OF EQUAL ENERGY. THE RESULT IS A LIST OF
11 !+ DIRECTIONS CORRESPONDING TO EACH BAND AND
12 !+ FOR EACH LIQUID BOUNDARY NODE.
13 !
14 !history N.DURAND (HRW)
15 !+ August 2017
16 !+ V7P3
17 !+ New. Computes BDALEs: spatially varying directional components
18 !
19 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
20 !| BDALE |<--| ISO ENERGY DIRECTIONS
21 !| | | FOR EACH LIQUID BOUNDARY NODE
22 !| PMIN |-->| MINIMUM FREQUENCY FOR SPECTRUM
23 !| PMAX |-->| MAXIMUM FREQUENCY FOR SPECTRUM
24 !| TETMIN |-->| MAXIMUM VALUE FOR THE PROPAGATION ANGLE
25 !| TETMAX |-->| MAXIMUM VALUE FOR THE PROPAGATION ANGLE
26 !| NPALE |-->| NUMBER OF DISCRETISATION BAND FOR PERIODS
27 !| NDALE |-->| NUMBER OF DISCRETISATION BAND FOR DIRECTIONS
28 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
29 !
30  USE bief
32  USE declarations_telemac, ONLY : kinc,ksort
33  USE interface_artemis, ONLY : stwc1
34  USE interface_parallel, ONLY : p_max
35 !
37  IMPLICIT NONE
38 !
39 !-----------------------------------------------------------------------
40 !
41  INTEGER :: IERR
42  INTEGER :: IPTFR
43  INTEGER :: JTMP,IMIN,IMAX
44  INTEGER :: NPASF,NPASD
45  INTEGER :: IDALE,IDD,IFF,I
46  INTEGER :: ISPEC
47 !
48  DOUBLE PRECISION :: FMIN,FMAX
49  DOUBLE PRECISION :: DTETA,DF,DTETA2
50  DOUBLE PRECISION :: SUMB,POIDS,SUMD,VAR,SUMICI
51 !
52  INTEGER,ALLOCATABLE :: IDTWC(:)
53 !
54  DOUBLE PRECISION,ALLOCATABLE :: SDIRTWC(:,:)
55  DOUBLE PRECISION,ALLOCATABLE :: SDIR(:)
56 !
57  INTEGER :: IDDMINART
58  DOUBLE PRECISION,ALLOCATABLE :: DIR_ART(:)
59  DOUBLE PRECISION,ALLOCATABLE :: SDIR_ART(:)
60 !
61  INTEGER :: PRINTMSG
62 !
63 !=======================================================================
64 !
65 ! MIN/MAX FREQUENCY
66  fmin = 1.d0 / pmax
67  fmax = 1.d0 / pmin
68 !
69 ! NUMBER OF INTEGRATION INTERVALS FOR THE TRAPEZOIDS METHOD : DIRECTIONS
70  npasd = 500*ndale
71 !
72 ! NUMBER OF INTEGRATION INTERVALS FOR THE TRAPEZOIDS METHOD : FREQUENCIES
73  npasf = 500*npale
74 !
75 ! WIDTH OF AN INTEGRATION INTERVAL : DIRECTIONS
76  dteta = (tetmax-tetmin)/float(npasd)
77 !
78 ! WIDTH OF AN INTEGRATION INTERVAL : FREQUENCIES
79  df = (fmax-fmin)/float(npasf)
80 !
81 !=======================================================================
82 !
83  ALLOCATE(idtwc(2*ndale+2),stat=ierr)
84  CALL check_allocate(ierr,'TWCAL2:IDTWC')
85  ALLOCATE(sdirtwc(nspec,ndir+1),stat=ierr)
86  CALL check_allocate(ierr,'TWCAL2:SDIRTWC')
87  ALLOCATE(sdir(npasd+1),stat=ierr)
88  CALL check_allocate(ierr,'TWCAL2:SDIR')
89  ALLOCATE(dir_art(npasd+1),stat=ierr)
90  CALL check_allocate(ierr,'TWCAL2:DIR_ART')
91  ALLOCATE(sdir_art(npasd+1),stat=ierr)
92  CALL check_allocate(ierr,'TWCAL2:SDIR_ART')
93 !
94  printmsg = 0
95 !
96 !-----------------------------------------------------------------------
97 !
98 ! INTEGRAL OF THE TOMAWAC SPECTRUM (TRAPEZOIDS METHOD)
99 ! FOR EACH DIRECTION COMPONENT
100 !
101  DO ispec = 1,nspec
102 !
103  DO idd = 1,ndir+1
104  sumd=0.d0
105  DO iff = 1,npasf+1
106 ! IF FREQ AND/OR DIR ON THE BOUNDARY OF THE INTEGRATION DOMAIN : CONTRIBUTION/2
107  poids=1.d0
108  IF(iff.EQ.1 .OR. iff.EQ.npasf+1) poids=0.5d0*poids
109  IF(idd.EQ.1 .OR. idd.EQ.ndir+1) poids=0.5d0*poids
110  var=stwc1(fmin+float(iff-1)*df,s_tom%DIR(idd),s_tom,ispec)
111  sumd = sumd + poids*var*df
112  ENDDO
113  sdirtwc(ispec,idd) = sumd
114 ! IF(DEBUG.GT.0) WRITE(LU,*) ISPEC,S_TOM%DIR(IDD),
115 ! & SDIRTWC(ISPEC,IDD)
116  ENDDO
117 !
118  ENDDO ! ISPEC = 1,NSPEC
119  IF(debug.GT.0)
120  & WRITE(lu,*) '< TWCAL2: SDIRTWC COMPUTED FOR 1 TO ',nspec
121 !
122 !-----------------------------------------------------------------------
123 !
124 ! FOR THE DIRECTION RANGE SPECIFIED IN ARTEMIS
125 !
126 ! INDICES COVERING THE RANGE
127  IF(tetmax-tetmin.EQ.360.d0) THEN
128  imin = 1
129  ELSEIF(tetmin.LE.0) THEN
130  DO jtmp = 1,ndir+1
131  IF(s_tom%DIR(jtmp).LE.tetmin+360.d0) imin = jtmp
132  ENDDO
133  ELSE
134  DO jtmp = 1,ndir+1
135  IF(s_tom%DIR(jtmp).LE.tetmin) imin = jtmp
136  ENDDO
137  ENDIF
138 !
139  IF(tetmax-tetmin.EQ.360.d0) THEN
140  imax = ndir+1
141  ELSEIF(tetmax.LT.0) THEN
142  DO jtmp = ndir+1,1,-1
143  IF(s_tom%DIR(jtmp).GE.tetmax+360.d0) imax = jtmp
144  ENDDO
145  ELSE
146  DO jtmp = ndir+1,1,-1
147  IF(s_tom%DIR(jtmp).GE.tetmax) imax = jtmp
148  ENDDO
149  ENDIF
150 !
151  IF(debug.GT.0) WRITE(lu,*)'< TWCAL2: TETMIN,IMIN,S_TOM%DIR(IMIN)'
152  IF(debug.GT.0) WRITE(lu,*) tetmin,imin,s_tom%DIR(imin)
153  IF(debug.GT.0) WRITE(lu,*)'< TWCAL2: TETMAX,IMAX,S_TOM%DIR(IMAX)'
154  IF(debug.GT.0) WRITE(lu,*) tetmax,imax,s_tom%DIR(imax)
155 !
156 !=======================================================================
157 !
158 ! FOR EACH ARTEMIS BOUNDARY NODE
159 !
160  DO iptfr = 1 , nptfr
161 !
162  IF (lihbor%I(iptfr).EQ.kinc .OR. lihbor%I(iptfr).EQ.ksort) THEN
163 !
164  printmsg = printmsg + 1
165 !
166 !-----------------------------------------------------------------------
167 !
168  IF(debug.GT.0 .AND. printmsg.EQ.1) THEN
169  WRITE(lu,*)
170  WRITE(lu,*) '> TWCAL2: ',
171  & 'STARTS INTERPOLATION OF DIR. DISTRIBUTION TO ARTEMIS NODES'
172  ENDIF
173 !
174 ! INTERPOLATES THE SPECTRUM FOR INCIDENT BOUNDARY POINTS, ONLY
175 !
176  DO idd = 1,npasd+1
177  sdir(idd) = 0.d0
178  ENDDO
179 !
180  CALL fasp_sp( s_tom%XOUTER,s_tom%YOUTER,sdirtwc,nspec,
181  & x(mesh%NBOR%I(iptfr)),y(mesh%NBOR%I(iptfr)),sdir,iptfr )
182 !
183  CALL stwc2( imin,imax,npasd+1,dir_art,sdir )
184  dteta2 = dir_art(2)-dir_art(1)
185 !
186  IF(debug.GT.0 .AND. printmsg.EQ.1)
187  & WRITE(lu,*) '< TWCAL2: INTERPOLATION TO ARTEMIS NODES ENDED'
188 !
189 !-----------------------------------------------------------------------
190 !
191 ! REORGANISE THE ARTEMIS SPECTRUM PRIOR TO SLICING
192 !
193  iddminart = 1
194  DO idd = 2,npasd
195  IF(sdir(idd).LT.sdir(iddminart)) iddminart = idd
196  ENDDO
197 !
198 ! REORGANISED SPECTRUM IN SDIR_ART(NPASD+1)
199  DO idd = 1,npasd
200  i=idd-iddminart+1
201  IF (idd.LT.iddminart) i=i+npasd
202  sdir_art(i)=sdir(idd)
203  ENDDO
204  sdir_art(npasd+1)=sdir_art(1)
205 !
206 ! INTEGRAL OF THE SPECTRUM (TRAPEZOIDS METHOD)
207 !
208  sumb = 0.d0
209  DO idd = 1,npasd+1
210  poids=1.d0
211  IF(idd.EQ.1 .OR. idd.EQ.npasd+1) poids=0.5d0*poids
212  sumb = sumb + poids*sdir_art(idd)*dteta2*degrad
213  ENDDO
214 !
215 ! SIGNIFICANT WAVE HEIGHT CORRESPONDING TO TOTAL ENERGY STORAGE
216 !
217  hscal=4.d0*sqrt(sumb)
218 !
219  IF(debug.GT.0 .AND. printmsg.EQ.1) THEN
220  WRITE(lu,*) '> TWCAL2: STARTS SPECTRUM DISCRETISATION'
221  WRITE(lu,*) ' DIRECTIONS'
222  ENDIF
223 !
224 ! =======================================================
225 ! DIRECTIONS
226 ! =======================================================
227 ! DIVIDES THE SPECTRUM INTO 2*NDALE BANDS OF EQUAL ENERGY
228  sumb = sumb/float(2*ndale)
229 !
230 ! FIRST TERM OF DIRECTION DISCRETIZATION
231  i=1
232  idtwc(i) = 1
233 !
234 ! IDENTIFIES THE ANGLES EVERY (I)*SUMB (I=1,NDALE)
235  sumici = 0.d0
236 !
237  DO idd = 1,npasd+1
238  sumici = sumici + sdir_art(idd)*dteta2*degrad
239 !
240 ! CHECKS IF SUMB = OVERALL SUMB/(2 NDALE) IS REACHED AND SAVES TETA
241 !
242  IF (sumici.GE.sumb*float(i) .OR. idd.EQ.npasd+1) THEN
243  i=i+1
244  idtwc(i) = idd
245  ENDIF
246  ENDDO
247 !
248 ! IDTWC HOLDS INDICES FOR THE DIRECTIONS TO BE CONSIDERED IN ARTEMIS
249 !
250 ! THESE ARE BASED ON CONTINUOUS SPECTRUM AND MAY NEED ADJUSTING,
251 ! IF PRINTED OUT FOR CHECKS, AS DONE FOR BDALE BELOW
252 !
253  DO i=1,2*ndale-1,2
254  idale = (i+1)/2
255  bdale%ADR(idale)%P%R(iptfr)=
256  & float(idtwc(i+1)-1)*dteta2+dir_art(iddminart)
257  IF(bdale%ADR(idale)%P%R(iptfr).GT.dir_art(npasd+1))
258  & bdale%ADR(idale)%P%R(iptfr)=
259  & bdale%ADR(idale)%P%R(iptfr)-dir_art(npasd+1)+dir_art(1)
260 !
261  IF(bdale%ADR(idale)%P%R(iptfr).GT.180d0)
262  & bdale%ADR(idale)%P%R(iptfr)=
263  & bdale%ADR(idale)%P%R(iptfr)-360.d0
264  ENDDO
265 !
266  ENDIF ! LIHBOR%I(IPTFR).EQ.KINC
267 !
268 !-----------------------------------------------------------------------
269 !
270  IF(debug.GT.0 .AND. printmsg.EQ.1)
271  & WRITE(lu,*) '< TWCAL2: SPECTRUM DISCRETISATION ENDED'
272 !
273  ENDDO ! IPTFR = 1 , NPTFR
274 !
275  DEALLOCATE(idtwc,sdirtwc,sdir)
276  DEALLOCATE(dir_art,sdir_art)
277 !
278 !-----------------------------------------------------------------------
279 !
280 ! FINDS BOUNDARY NODE CLOSEST TO REFERENCE POINT
281 ! (USED FOR PRINTOUTS ONLY IN ARTEMIS.F)
282 !
283  iptfr_ref = 0
284  IF (nptfr.GT.0) CALL twcclosest
285  IF (ncsize.GT.1) iptfr_ref = p_max(iptfr_ref)
286  IF(debug.GT.0)
287  & WRITE(lu,*) 'CLOSEST: REFERENCE BOUNDARY POINT:',iptfr_ref
288 !
289 !-----------------------------------------------------------------------
290 !
291  RETURN
292  END
type(bief_obj), target lihbor
double precision function stwc1(F, DIR, SPEC, I)
Definition: stwc1.f:7
subroutine twcal2
Definition: twcal2.f:4
double precision, dimension(:), pointer y
subroutine stwc2(IMIN, IMAX, N, DIR2, SDIR)
Definition: stwc2.f:7
type(bief_obj), target bdale
type(bief_mesh), target mesh
integer, parameter kinc
double precision, dimension(:), pointer x
subroutine fasp_sp(XRELV, YRELV, ZRELV, NP, X, Y, Z, I)
Definition: fasp_sp.f:7
integer, parameter ksort
subroutine twcclosest
Definition: twcclosest.f:4
Definition: bief.f:3