forked from eriove/pdstrip
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathpdstrip.f90
2224 lines (2146 loc) · 118 KB
/
pdstrip.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
!Last modification December 2014
module pdstripsubprograms
implicit none
logical:: output=.false.
integer, parameter:: nmumax=100 !max. number of wave angles, regular waves
integer, parameter:: nofmax=240 !max. number of offset points on a section and the free surface
integer, parameter:: nfremax=52 !max. number of frequencies for hydrodynamic section calculations
integer, parameter:: nsemax=100 !max. number of offset sections
integer, parameter:: nprmax=100 !max. number of pressure points per section
integer, parameter:: nfs=55 !panel number on free suface on each side of the body
integer, parameter:: nomma=200 !max. no. of frequencies in transfer function calculation
integer, parameter:: nvmax=50 !max. no. of ship speeds
integer, parameter:: nbmax=200 !max. no. of points where motions should be determined
logical:: sym !true if all offset sections are symmetric
logical:: subm !true if a section is fully submerged
logical:: lsect !true if hydrodynamic section calculations required
logical:: lsign !t to compute significant amplitudes in natural seaways; f otherw.
logical:: ls !true if intersection forces and moments to be determined
integer:: ise,nse !index and number of sections
integer:: npres,npres1 !number of pressure points per section; npres1=max(npres,1)
logical:: lpres=.true. !true if pressure transfer functions to be written in pdstrip.out
integer:: nof(nsemax) !number of offset points on a section
integer:: ngap !number of `gaps' of a section
integer:: gap(10) !offset point index where a gap begins
integer:: nmu !number of wave angles in hydrodynamic section calculations
real, parameter:: pi=3.14159
real:: wangl(nmumax) !wave angles in hydrodynamic section calculations
real:: g !gravity acceleration
real:: rho !water density
real:: zwl !z coordinate of the still waterline
real:: zbot !z coordinate of the water bottom
real:: tg !reference draft (same for all sections)
real:: waven !wave number
real:: x(0:nsemax+1) !x coordinates of sections
real:: yof(nsemax,nofmax) !y coordinates of offset points (section index, point intex)
real:: zof(nsemax,nofmax) !z coordinates of offset points (section index, point intex)
real:: bcross(nsemax) !maximum breadth of sections used for cross-flow resistance
real:: tcross(nsemax) !maximum draft of sections used for cross-flow resistance
real:: area(0:nsemax) !submerged area of sections
real:: ys(0:nsemax+1) !y coordinate of center of gravity of section area
real:: zs(0:nsemax+1) !z coordinate of center of gravity of section area
real:: zdrift !z coordinate for water particle drift velocity
character(80):: text !text describing the ship and its loading case
character(80):: offsetfile !name of the file containing the offset point data
interface operator(.mprod.); module procedure prodmatr; end interface
interface operator(.vprod.); module procedure vectprod; module procedure cvectprod;
end interface !vector product
contains
subroutine sectiondata !input section data; call of sectionhydrodynamics
! for deep and shallow water; section motion and wave excitation; with and without symmetry
integer:: i,ii,j !indices
integer:: im1 !i-1 except for i=1
integer:: nfre=52 !number of frequencies
real:: fp(nfremax) !wave frequency parameter omega^2*t/g
real:: om(nfremax) !circular wave frequency
real:: da !section area increment
real:: yof1(nofmax),zof1(nofmax) !y, z coordinates of offset points of one section only
complex:: addedm(3,3,nfremax) !complex added mass matrix of a section = a+d/(i*omega_e)
complex:: diff(3,nmumax,nfremax) !diffraction force amplitude vector
complex:: frkr(3,nmumax,nfremax) !Froude-Krilow force amplitude vector
complex:: pr(nprmax,3+nmumax,nfremax) !pressure amplitude at pressure points on section contour
fp=(/0.01,0.015, 0.020,.025,.031,.04,.05,.063,.08,.10,.12,.15,.19,.24,.31,.39, &
.47,.55,.62,.70,.80,.90,1.0,1.1,1.25,1.4,1.55,1.7,1.9,2.2,2.4,2.7,3.0,3.3,3.6,4.0, &
4.5,5.2,6.,7.,8.,9.,10.5,12.0,13.5,15.,17.,20.,25.,30.,40.,50./) !standard frequency parameters
!nfre=3; fp(1:3)=(/0.1,1.,10./) !only for rapid testing
read(5,'(a)')text
write(6,'(/1x,1a80)')text
write(6,'(a)')' In the following input data x,y,z are directed forward, to port side, up'
read(5,*)g,rho,zwl,zbot,zdrift
write(6,'(/'' Gravity acceleration '',f 8.3)')g
write(6,'( '' Water density '',f 8.3)')rho
write(6,'( '' z of still waterline '',f 8.3)')zwl
write(6,'( '' z of water bottom '',g12.3)')zbot
write(6,'( '' z for water drift velocity '',f 8.3)')zdrift
read(5,*)nmu,wangl(1:nmu)
write(6,'( '' Wave encounter angles '',(10f8.3))')wangl(1:nmu)
read(5,*)offsetfile
write(6,'(a,a)')' Offset file name: ',trim(offsetfile)
if (nmu>nmumax) call stop1('Too many wave angles')
if (zbot>0) write(6,*)'*** positive bottom z coordinate appears unreasonable'
zwl=-zwl; zbot=-zbot; zdrift=-zdrift !change from input to output coord. system
if (lsect) then
write(20,'(a)')text; write(20,*)nfre !write initial data on file sectionsresults
endif
open(9,file=offsetfile,status='old')
read(9,*)nse,sym,tg !number of sections, symmetry?, reference draft
write(6,'(a,i8 )') ' Number of sections ',nse
write(6,'(a,a )') ' Symmetry? ',merge(' yes',' no',sym)
write(6,'(a,f8.3)') ' Reference draft ',tg
if (nse>nsemax) call stop1('Too many sections')
if (tg<=0) call stop1('Reference draft should be positive')
if (sym.and.npres/=0) call stop1('Pressure calculation requires sym=.false')
if (npres==1.or.npres==2) &
call stop1('<=2 pressure points per section not admitted. 0 or >=15 recommended')
om(1:nfre)=sqrt(fp(1:nfre)*g/tg) !circular fr. resulting from frequency parameters
Sections: do ise=1,nse
read(9,*)x(ise),nof(ise),ngap,gap(1:ngap) !section x coord., no. of offsets, gap data
read(9,*)yof(ise,1:nof(ise)),zof(ise,1:nof(ise)) !offsets read from file 9
write(6,'(" Section no.",i3," at x=",f10.3)') ise,x(ise)
if (ngap>0) write(6,'(a,10i5)') ' Gaps following point no.',gap(1:ngap)
write(6,'(" y:",8f10.3)') yof(ise,1:nof(ise))
write(6,'(" z:",8f10.3)') zof(ise,1:nof(ise))
if (sym.and.yof(ise,1)/=0) call stop1('First y must be 0 for symmetrical body')
yof(ise,1:nof(ise))=-yof(ise,1:nof(ise))
zof(ise,1:nof(ise))=-zof(ise,1:nof(ise)) !change from input to output coord. system
area(ise)=0; ys(ise)=0; zs(ise)=0 !compute area etc. gaps covered 2 times if sym=.false.
Offsets: do i=1,nof(ise)
im1=merge(nof(ise),i-1,i==1)
if (ngap>0.and.any(gap(1:ngap)==im1)) cycle
da=(yof(ise,im1)*zof(ise,i)-zof(ise,im1)*yof(ise,i))/2
area(ise)=area(ise)+da
ys(ise)=ys(ise)+da*(yof(ise,im1)+yof(ise,i))/3
zs(ise)=zs(ise)+da*(zof(ise,im1)+zof(ise,i))/3
enddo Offsets
bcross(ise)=maxval(yof(ise,1:nof(ise)))-minval(yof(ise,1:nof(ise)))
tcross(ise)=maxval(zof(ise,1:nof(ise)))-minval(zof(ise,1:nof(ise)))
if (ngap>0) tcross(ise)=tcross(ise)-merge(1.,0.5,sym)* &
sum(abs(zof(ise,gap(1:ngap))-zof(ise,gap(1:ngap)+1)))
if (sym) then
da=yof(ise,nof(ise))*(zof(ise,nof(ise))-zof(ise,1))
area(ise)=2*area(ise)+da
ys(ise)=0
zs(ise)=2*zs(ise)+da*(2*zof(ise,nof(ise))+zof(ise,1))/3
bcross(ise)=2*bcross(ise)
endif
ys(ise)=ys(ise)/area(ise)
zs(ise)=zs(ise)/area(ise)
if ((.not.sym.and.(yof(ise,1)-yof(ise,nof(ise)))**2+(zof(ise,1)-zof(ise,nof(ise)))**2<1e-4) &
.or.(sym.and.abs(yof(ise,nof(ise)))<1e-3)) then
subm=.true.; write(6,*) 'fully submerged'
else
subm=.false.
if (abs(zof(ise,nof(ise))-zwl)>0.001) call stop1('Last z must be equal to zwl')
if (.not.sym.and.abs(zof(ise,1)-zwl)>0.001) call stop1('First z must be equal to zwl')
endif
yof1(1:nof(ise))=yof(ise,1:nof(ise)); zof1(1:nof(ise))=zof(ise,1:nof(ise))
if (lsect) then
output=ise==18
call sectionhydrodynamics(nof(ise),yof1,zof1,nfre,om,addedm,diff,frkr,pr)
Frequencies: do i=1,nfre
write(20,*)om(i),nmu,wangl(1:nmu)*pi/180 !store on file sectionresults: wave data
write(20,*)(om(i)**2*addedm(j,1:3,i),j=1,3) !radiation force
write(20,*)(diff(j,1:nmu,i),j=1,3) !diffraction force
write(20,*)(frkr(j,1:nmu,i),j=1,3) !Froude-Krilow force
if (npres>0) write(20,*)(ii,pr(ii,1:3+nmu,i),ii=1,npres) !pressures
!if(i==16)print *,'pr4',om(i),pr(1,2,16)
enddo Frequencies
endif
enddo Sections
if (lsect) &
write(20,*)npres+sum(wangl(1:nmu))+sum(x(1:nse))+sum((/(yof(ise,1:nof(ise)),ise=1,nse)/)) &
+sum((/(zof(ise,1:nof(ise)),ise=1,nse)/))+g+rho+zwl+1/zbot !test number for unchanged input data
end subroutine sectiondata
subroutine sectionhydrodynamics(nof1,yof1,zof1,nfre,om,addedm,diff,frkr,pr)
integer, intent(in):: nof1 !number of offset points on one section
integer:: nfre !number of wave frequencies
integer:: i,ii !indices
integer:: iom !index of frequency
integer:: ileft(nprmax) !index of offset point following a pressure point
real:: yof1(nofmax),zof1(nofmax) !y and z coordinates of offset points of one section
real:: om(nfremax) !circular frequency
real:: girth(nofmax) !girth length of section up to offset points
real:: u !girth length up to current pressure point
real:: aint(nprmax) !interpolation factor for pressure point between offsets
complex:: addedm(3,3,nfremax) !complex added mass matrices of one section, all frequencies
complex:: diff(3,nmumax,nfremax) !diffraction force amplitudes of one section, all frequencies
complex:: frkr(3,nmumax,nfremax) !Froude-Krilow force amplitudes, one section, all frequ.
complex:: addedmprel(2,2) !preliminary added mass matrix (1 or 2 degrees of freedom)
complex:: diffprel(2,nmumax) !preliminary diffraction force
complex:: frkrprel(2,nmumax) !preliminary Froude-Krilow force
complex:: pr(nprmax,3+nmumax,nfremax) !pressure amplitudes
if (nfre.gt.nfremax) call stop1('Too many frequencies')
call testsuitability(nof1,yof1,zof1) !discretization of section suitable?
girth(1)=0.
do i=2,nof1
if (ngap>0.and.any(gap(1:ngap)==i-1)) then
girth(i)=girth(i-1); cycle
endif
girth(i)=girth(i-1)+sqrt((yof1(i)-yof1(i-1))**2+(zof1(i)-zof1(i-1))**2)
enddo
PressurePoints: do ii=1,npres !ii index of pressure points
u=(ii-1.0)/(npres-1)*girth(nof1) !u= girth length to pressure point ii
i=which(nof1-1,u<=girth(2:nof1))+1 !pressure point between ileft-1 and ileft
ileft(ii)=i
aint(ii)=(u-girth(i-1))/(girth(i)-girth(i-1)) !intp=val(ileft-1)+aint*(val(ileft)-val(ileft-1))
write(24,*)yof1(i-1)+aint(ii)*(yof1(i)-yof1(i-1)),zof1(i-1)+aint(ii)*(zof1(i)-zof1(i-1))
enddo PressurePoints
Frequencies: do iom=1,nfre
Symmetry: if (sym) then
addedm(1:3,1:3,iom)=0
call addedmassexcitations(nof1,yof1,zof1,om(iom),1,addedm(2,2,iom),diff(2,1:nmu,iom), &
frkr(2,1:nmu,iom))
call addedmassexcitations(nof1,yof1,zof1,om(iom),2,addedmprel,diffprel,frkrprel)
addedm(1:3:2,1:3:2,iom)=addedmprel
diff(1:3:2,1:nmu,iom)=diffprel(1:2,1:nmu)
frkr(1:3:2,1:nmu,iom)=frkrprel(1:2,1:nmu)
else symmetry
call addedmassexcitationa(nof1,yof1,zof1,om(iom),addedm(1,1,iom),diff(1,1,iom), &
frkr(1,1,iom),pr(1,1,iom),ileft,aint)
endif Symmetry
enddo Frequencies
end subroutine sectionhydrodynamics
function pq(dy,dz) !pq=source (flux 2 pi) potential. dy,dz vector from source point to actual point
real:: pq,dy,dz
pq=0.5*log(dy**2+dz**2)
end function pq
function pqn(dy1,dz1,dy2,dz2) !pqn=flux (area/time) between p1, p2 due to source of strength 2 pi.
real:: pqn,dy1,dz1,dy2,dz2 !pqn positive for flow to the left when looking from p1 to p2
pqn=-atan2(dy2*dz1-dz2*dy1,dy1*dy2+dz1*dz2)
end function pqn
function pqb(yof,zof,yq,zq) !pq possibly with bottom mirror source
real:: pqb,yof,zof,yq,zq
pqb=pq(yof-yq,zof-zq)
if (zbot<1e6) pqb=pqb+pq(yof-yq,zof-(2*zbot-zq))
end function pqb
function pqnb(yof1,zof1,yof2,zof2,yq,zq) !pqn possibly with bottom mirror source
real:: pqnb,yof1,zof1,yof2,zof2,yq,zq
pqnb=pqn(yof1-yq,zof1-zq,yof2-yq,zof2-zq)
if (zbot<1e6) pqnb=pqnb+pqn(yof1-yq,zof1-(2*zbot-zq),yof2-yq,zof2-(2*zbot-zq))
end function pqnb
function pqs(yof,zof,yq,zq,mp) !=pq with y=0 mirror source
integer:: mp !and possibly with bottom mirror source
real:: pqs,yof,zof,yq,zq
pqs=pq(yof-yq,zof-zq)+mp*pq(yof+yq,zof-zq)
if (zbot<1000.) pqs=pqs+pq(yof-yq,zof-(2*zbot-zq))+mp*pq(yof+yq,zof-(2*zbot-zq))
end function pqs
function pqsn(yof1,zof1,yof2,zof2,yq,zq,mp) !=pqn with y=0 mirror source
integer:: mp
real:: pqsn,yof1,zof1,yof2,zof2,yq,zq !and possibly with bottom mirror sources
pqsn=pqn(yof1-yq,zof1-zq,yof2-yq,zof2-zq)+mp*pqn(yof1+yq,zof1-zq,yof2+yq,zof2-zq)
if (zbot<1000.) pqsn=pqsn+pqn(yof1-yq,zof1-(2*zbot-zq),yof2-yq,zof2-(2*zbot-zq)) &
+mp*pqn(yof1+yq,zof1-(2*zbot-zq),yof2+yq,zof2-(2*zbot-zq))
end function pqsn
subroutine addedmassexcitationa(nof1,yof1,zof1,om,addedm1,diff1,frkr1,pr,ileft,aint)
! added mass and wave excitation, deep and shallow water, asymmetrical section allowed
! bottom at zbot, surface at zwl. determines also pressure at interpolated points.
! for deep water use zbot=1.e6. water density 1 assumed.
! computes complex added mass = mass-i/omegae*damping for sway,heave,roll (3*3 matrix)
! and diffraction and Froude-Krilow forces diff,frkr for nmu wave angles wangl (degree)
! section offset points from starboard waterline point over basis to port wl point
! +y to starboard; +z downwards. first and last point should be on waterline.
integer, intent(in):: nof1 !no. of offset points on a section
integer:: ileft(nprmax) !index of offset point following a pressure point
integer:: k,k1 !offset point indices
integer:: nf !number of free-surface panels (undamped) on either side
integer:: iw !index of wave angle
integer:: is !points to row if equation solver detects singularity
integer:: i,ii,j,l,m !indices
real:: aint(nprmax) !interpolation factor for pressure point between offsets
real:: yof1(nofmax),zof1(nofmax) !y and z coordinate of offset points on one section
real:: t !water depth
real:: fact1,fact2,factor,fact !factors
real:: dy,dz !range of contour segment in y and z direction
real:: sinw !sin(wave angle)
real:: y1 !y coordinate of a section offset point
real:: om !circular frequency
real:: h !point distance on free surface near to body
real:: hend !point distance on free surface far from body
real:: yq(0:nofmax),zq(0:nofmax) !y and z coordinates of sources
complex:: a(0:nofmax,0:nofmax+3+nmumax)!coefficient matrix to determine source strengths
complex:: ciom !i*omega
complex:: detl !log(determinant) of the coefficient matrix a
complex:: f !force amplitude
complex:: addedm1(3,3) !complex added mass matrix for one frequency, one section
complex:: pr(nprmax,3+nmumax) !pressure amplitudes at pressure points
complex:: phiatoffsets(nofmax) !pressure amplitudes at offset points
complex:: diff1(3,nmu) !diffraction force amplitude, one frequency, one section
complex:: frkr1(3,nmu) !Froude-Krilow force amplitude, one frequency, one section
complex:: zw1,zw2,zw1a,zw2a !complex intermediate numbers
if (nof1.gt.100) call stop1('*** Too many (>100) points on section.')
if (nof1+nfs>nofmax) call stop1('*** Too many body + waterline panels')
if (nmu.gt.nmumax) call stop1('*** Too many angles')
if (.not.subm.and.(abs(zof1(1)-zwl)>0.001.or.abs(zof1(nof1)-zwl) >0.001)) &
call stop1('*** First or last section point not on waterline')
addedm1=0.; diff1=0; pr=0
TwoFreeSurfaceDiscretisations: do nf=25,28,3
ciom=(0.,1.)*om
t=zbot-zwl !t=water depth
waven=max(om**2/g,om/sqrt(g*t)) !wave number
if (waven*t<6) then
do i=1,10 !for shallow water: iterative determination of wave number
waven=waven-(waven*tanh(waven*t)-om**2/g)/(tanh(waven*t)+waven/cosh(waven*t)**2)/2
enddo
fact1=g*exp( waven*t)/2/cosh(waven*t); fact2=g*exp(-waven*t)/2/cosh(waven*t)
else
fact1=g; fact2=0.
endif
h=sqrt((yof1(nof1)-yof1(nof1-1))**2+(zof1(nof1)-zof1(nof1-1))**2) !near distance on free surface
hend=2*pi/(12*waven) !far-off distance of points on free surface=wavelength/12
i=0 ! compute free surface grid points yof1,zof1 and all source points yq,zq
BodySegments: do k=1,nof1-1
if (ngap>0.and.any(gap(1:ngap)==k)) cycle
i=i+1
yq(i)=(yof1(k)+yof1(k+1))/2-(zof1(k+1)-zof1(k))/20 !source points within section
zq(i)=(zof1(k)+zof1(k+1))/2+(yof1(k+1)-yof1(k))/20
enddo BodySegments
FreeSurfaceSegments: do k=nof1,nof1+nfs-1
h=min(h*1.5,hend)
yof1(k+1)=yof1(k)-h
yof1(nfs+k+1)=merge(yof1(1),yof1(nfs+k),k==nof1)+h
zof1(k+1)=zwl
zof1(nfs+k+1)=zwl
yq(k-ngap)=yof1(k)-0.5*h
zq(k-ngap)=zwl-1.0*h
yq(nfs+k-ngap)=merge(yof1(1),yof1(nfs+k),k==nof1)+0.5*h
zq(nfs+k-ngap)=zwl-1.0*h
enddo FreeSurfaceSegments
yq(0)=0 !zero source
zq(0)=zwl-abs(yof1(nof1+nfs))/2
a(0:nof1-ngap+2*nfs-1,nof1-ngap+2*nfs:nof1-ngap+2*nfs+2+nmu)=0. !inhomogeneous terms = 0
frkr1=0
!Body boundary condition (panel integrals) and Fr.-Krilow forces
j=0 !j=index of body segments excluding gaps
BodySegm: do k=1,nof1-1
if (ngap>0.and.any(gap(1:ngap)==k)) cycle
j=j+1
dy=yof1(k+1)-yof1(k)
dz=zof1(k+1)-zof1(k)
a(j,nof1-ngap+2*nfs )=-ciom*dz
a(j,nof1-ngap+2*nfs+1)=+ciom*dy
a(j,nof1-ngap+2*nfs+2)=+ciom*(yof1(k+1)**2-yof1(k)**2+zof1(k+1)**2-zof1(k)**2)/2
WaveAngles: do iw=1,nmu
sinw=sin(wangl(iw)*pi/180)
zw1=cmplx(-waven*dz,waven*sinw*dy)
if (abs(zw1)>51.) call stop1('Frequency too high')
if (abs(zw1).lt.0.001) zw1=0.001
zw1a=-conjg(zw1)
zw2=exp(waven*cmplx(-zof1(k)+zwl,yof1(k)*sinw))* &
merge((exp(zw1)-1.)/zw1,1.+0.5*zw1,abs(zw1 )>0.01)
zw2a=exp(waven*cmplx( zof1(k)-zwl,yof1(k)*sinw))* &
merge((exp(zw1a)-1.)/zw1a,1.+0.5*zw1a,abs(zw1a)>0.01)
a(j,nof1-ngap+2*nfs+2+iw)=-fact1*waven/om*cmplx(dy,sinw*dz)*zw2*(0.,1.)
frkr1(1,iw)=frkr1(1,iw)+fact1*zw2*dz*rho
frkr1(2,iw)=frkr1(2,iw)-fact1*zw2*dy*rho
frkr1(3,iw)=frkr1(3,iw)-fact1*zw2*(yof1(k+1)**2-yof1(k)**2+zof1(k+1)**2-zof1(k)**2)/2*rho
Shallow: if (fact2/=0) then
a(j,nof1-ngap+2*nfs+2+iw)=a(j,nof1-ngap+2*nfs+2+iw) &
+fact2*waven/om*cmplx(dy,-sinw*dz)*zw2a*(0.,1.)
frkr1(1,iw)=frkr1(1,iw)+fact2*zw2a*dz*rho
frkr1(2,iw)=frkr1(2,iw)-fact2*zw2a*dy*rho
frkr1(3,iw)=frkr1(3,iw)-fact2*zw2a*(yof1(k+1)**2-yof1(k)**2+zof1(k+1)**2-zof1(k)**2)/2*rho
endif Shallow
enddo WaveAngles
Columns: do i=0,nof1-ngap+2*nfs-1
a(j,i)=pqnb(yof1(k),zof1(k),yof1(k+1),zof1(k+1),yq(i),zq(i))
enddo Columns
enddo BodySegm
!Free surface condition near to body
FreeSSegmentsNear: do k=nof1,nof1+2*nf+1
k1=merge(k-nf+nfs-1,k,k>nof1+nf) !k1=k on y>0-side if k<=nof1+nf
y1=merge(yof1(1),yof1(k1),k==nof1+nf+1)
dy=yof1(k1+1)-y1
AllSources: do i=0,nof1-ngap+2*nfs-1 !2-point integration (.316,.684) correct in case of k1=i
a(k1-ngap,i)=pqnb(y1,zwl,yof1(k1+1),zwl,yq(i),zq(i))*merge(1,-1,k==k1) &
+om**2/g*abs(dy)*(pqb(y1+0.316*dy,zwl,yq(i),zq(i))+pqb(y1+0.684*dy,zwl,yq(i),zq(i)))/2
enddo AllSources
enddo FreeSSegmentsNear
!Free surface condition far from body (radiation condition)
FarSegments: do k=nof1+nf+1,nof1+2*nfs-nf-2
k1=merge(k+nf+1,k,k>=nof1+nfs)
dy=yof1(k1+1)-yof1(k1)
factor=(real(k1-(nof1+nf))/(nfs-nf))**2
if (k/=k1) factor=(real(k1-(nof1+nfs+nf))/(nfs-nf))**2
AllColumns: do i=0,nof1-ngap+2*nfs-1
a(k1-ngap,i)=waven*abs(dy)* &
(pqb(yof1(k1)+0.316*dy,zwl,yq(i),zq(i))+pqb(yof1(k1)+0.684*dy,zwl,yq(i),zq(i)))/2 &
-(0.,1.)*(pqb(yof1(k1+1)-factor*dy,zwl,yq(i),zq(i))-pqb(yof1(k1)-factor*dy,zwl,yq(i),zq(i)))
enddo AllColumns
enddo FarSegments
a(0,0:nof1-ngap+2*nfs-1)=1. ! sum of source strengths = 0
call simqcd(a(0,0),nof1-ngap+2*nfs,3+nmu,nofmax+1,is,1.e-5,detl) !solve equation system
if (is.ne.0) call stop1('*** Singular equation system. Coinciding offset points?')
!if(output.and.nf==28.and.abs(om-0.4288)<1e-3)print *,'Q',a(0:nof1-ngap+2*nfs-1,nof1-ngap+2*nfs+1)
!Bei L"ucke keine Unbekannte angesetzt!
!Added mass and diffraction calc. by pressure integration over section contour
Forces: do m=1,3 !for transverse force, vertical force, roll moment
MotionsAndWaves: do l=1,3+nmu !for sway, heave, roll
f=0 !force
AllBodySegments: do k=1,nof1-1
if (ngap>0.and.any(gap(1:ngap)==k)) cycle
dy=yof1(k+1)-yof1(k)
dz=zof1(k+1)-zof1(k)
if (m==1) then; fact=-dz !horizontal force
elseif (m==2) then; fact=+dy !vertical force
else; fact=+0.5*(yof1(k+1)**2-yof1(k)**2+zof1(k+1)**2-zof1(k)**2) !m==3; moment
endif
ForAllSources: do i=0,nof1-ngap+2*nfs-1
f=f+a(i,nof1-ngap+2*nfs-1+l)*fact*(pqb(yof1(k)+0.316*dy,zof1(k)+0.316*dz,yq(i),zq(i)) &
+pqb(yof1(k)+0.684*dy,zof1(k)+0.684*dz,yq(i),zq(i)))/2 !/2 for two-point integr.
enddo ForAllSources
enddo AllBodySegments
RadiationDiffraction: if (l<=3) then !added mass = -force /acceleration; pr.=-rho i om phi
addedm1(m,l)=addedm1(m,l)+0.5*ciom*f/(-om**2)*rho !*0.5 for average nf=25,28
else radiationdiffraction
diff1(m,l-3)=diff1(m,l-3)-0.5*ciom*f*rho
endif RadiationDiffraction
enddo MotionsAndWaves
enddo Forces
!Calculation of pressures
WithPressure: if (npres>0) then
Motions: do l=1,3+nmu
Offsets: do k=1,nof1
phiatoffsets(k)=0.
Sources: do i=0,nof1-ngap+2*nfs-1
phiatoffsets(k)=phiatoffsets(k)+a(i,nof1-ngap+2*nfs-1+l)*pqb(yof1(k),zof1(k),yq(i),zq(i))
enddo Sources
enddo Offsets
!print *,'pato',l,phiatoffsets(1:nof1)
PressurePoints: do ii=1,npres
pr(ii,l)=pr(ii,l)-ciom*rho*(phiatoffsets(ileft(ii)-1)+ & !pressure at pressure point
aint(ii)*(phiatoffsets(ileft(ii))-phiatoffsets(ileft(ii)-1)))/2 !sum for 2 nf values
enddo PressurePoints
enddo Motions
endif WithPressure
enddo TwoFreeSurfaceDiscretisations
end subroutine addedmassexcitationa
subroutine addedmassexcitations(nof1,yof1,zof1,om,ndof,addedm1,diff1,frkr1)
! added mass and wave excitation for deep and shallow water, symmetrical section
! bottom at zbot, surface at zwl. for deep water use zbot=1.e6. water density 1 assumed.
! complex added mass = mass-i/omegae*damping for vertical (ndof=1; 1x1 matrix) or
! horizontal/roll motion (ndof=2; 2x2 matrix) as well as vertical (ndof=1) diffraction
! and Froude-Krilow forces or horizontal force and roll moment (ndof=2)
! section offset points from basis (y=0)to port wl point
! +y to starboard; +z downwards. last point should be on waterline
integer, intent(in):: nof1 !number of offsets on one half section
integer:: ndof !number of degrees of freedom: 1 heave, 2 sway/roll
integer:: mp !+1 if ndof=1; -1 if ndof=2
integer:: i,j !indices
integer:: iw !index of wave angle
integer:: is !points to row if equation solver detects singularity
integer:: l !motion index
integer:: m !force or degree-of-freedom index
integer:: k !index of offset points on section
integer:: nf !number of free-surface panels (undamped)
real:: om !circular frequency
real:: yof1(nofmax),zof1(nofmax) !y and z coordinates of offset points of one section
real:: yq(0:nofmax),zq(0:nofmax) !y and z coordinates of sources
real:: t !water depth
real:: fact1,fact2,factor,fact !factors
real:: h !point distance on free surface near to body
real:: hend !point distance on free surface far from body
real:: dy,dz !range of contour segment in y and z direction
real:: sinw !sin(wave angle)
complex:: a(0:nofmax,0:nofmax+3+nmumax)!coefficient matrix for source strengths determination
complex:: ciom !i*om
complex:: detl !log(determinant) of a
complex:: f(2,2+nmumax) !force amplitude
complex:: addedm1(ndof,ndof) !complex added mass matrix
complex:: diff1(ndof,nmu) !diffraction force amplitude, one section
complex:: frkr1(ndof,nmu) !Froude-Krilow force amplitude, one section
complex:: zw1,zw2,zw1a,zw2a !intermediate values
if (nof1.gt.50) call stop1('*** Too many points on section. should be <=50')
if (ndof.ne.1.and.ndof.ne.2) call stop1('*** Incorrect ndof. should be 1 or 2')
if (nmu.gt.nmumax) call stop1('*** Too many angles')
if (nof1+nfs>nofmax) call stop1('*** Too many body + waterline panels')
addedm1=0.
diff1=0.
twofreesurfacediscretisations: do nf=25,28,3 !uses two values nf; results averaged
mp=3-2*ndof
ciom=(0.,1.)*om
t=zbot-zwl !t=water depth
waven=max(om**2/g,om/sqrt(g*t))
if (waven*t<6) then
do i=1,10 !for shallow water: iterative determination of wave number
waven=waven-(waven*tanh(waven*t)-om**2/g)/(tanh(waven*t)+waven/cosh(waven*t)**2)/2
enddo
fact1=g*exp( waven*t)/2/cosh(waven*t); fact2=g*exp(-waven*t)/2/cosh(waven*t)
else
fact1=g; fact2=0.
endif
h=sqrt((yof1(nof1)-yof1(nof1-1))**2+(zof1(nof1)-zof1(nof1-1))**2) !near distance on free surf.
hend=2*pi/(12*waven) !far-off distance of points on free surface=wavelength/12
!Compute free surface grid points yof1,zof1 and all source points yq,zq
i=0
BodySegments: do k=1,nof1-1
if (ngap>0.and.any(gap(1:ngap)==k)) cycle
i=i+1
yq(i)=(yof1(k)+yof1(k+1))/2-(zof1(k+1)-zof1(k))/20 !source points within section
zq(i)=(zof1(k)+zof1(k+1))/2+(yof1(k+1)-yof1(k))/20
enddo BodySegments
FreeSurfaceSegments: do k=nof1,nof1+nfs-1
h=min(h*1.5,hend)
yof1(k+1)=yof1(k)-h
zof1(k+1)=zwl
yq(k-ngap)=yof1(k)-0.5*h
zq(k-ngap)=zwl-1.0*h
enddo FreeSurfaceSegments
yq(0)=0
zq(0)=zwl-abs(yof1(nof1+nfs))/2
a(0:nof1-ngap+nfs-1,nof1-ngap+nfs:nof1-ngap+nfs-1+ndof+nmu)=0. !inhomogeneous terms = 0
frkr1(1:ndof,1:nmu)=0
!Body boundary condition (panel integrals) and Fr.-Krilow forces
j=0
BodySegm:do k=1,nof1-1
if (ngap>0.and.any(gap(1:ngap)==k)) cycle
j=j+1
dy=yof1(k+1)-yof1(k)
dz=zof1(k+1)-zof1(k)
if (mp.eq.1) then
a(j,nof1-ngap+nfs)=+ciom*dy
else
a(j,nof1-ngap+nfs )=-ciom*dz !inhomogeneous terms for body motion
a(j,nof1-ngap+nfs+1)=+ciom*(yof1(k+1)**2-yof1(k)**2+zof1(k+1)**2-zof1(k)**2)/2
endif
WaveAngles: do iw=1,nmu
sinw=sin(wangl(iw)*pi/180)
zw1=cmplx(-waven*dz,waven*sinw*dy)
!print *,'om',om,waven,dz,zw1
if (abs(zw1)>51) call stop1('Frequency too high')
if (abs(zw1).lt.0.001) zw1=0.001 ! stop 'frequency too low' seems unnecessary
zw1a=-conjg(zw1)
zw2 =exp(waven*cmplx(-zof1(k)+zwl,yof1(k)*sinw))* &
merge((exp(zw1 )-1.)/zw1,1.+0.5*zw1,abs(zw1)>0.01)
zw2a=exp(waven*cmplx( zof1(k)-zwl,yof1(k)*sinw))* &
merge((exp(zw1a)-1.)/zw1a,1.+0.5*zw1a,abs(zw1a)>0.01)
if (mp.eq.1) then !heave motion
a(j,nof1-ngap+nfs-1+ndof+iw)=-fact1*waven/om*real(cmplx(dy,sinw*dz)*zw2)*(0.,1.)
frkr1(1,iw)=frkr1(1,iw)-fact1*2*real(zw2 )*dy*rho
Shallow: if (fact2/=0) then
a(j,nof1-ngap+nfs-1+ndof+iw)=a(j,nof1-ngap+nfs-1+ndof+iw) &
+fact2*waven/om*real(cmplx(dy,-sinw*dz)*zw2a)*(0.,1.)
frkr1(1,iw)=frkr1(1,iw)-fact2*2*real(zw2a)*dy*rho
endif Shallow
else !sway/roll motion
a(j,nof1-ngap+nfs-1+ndof+iw)=-fact1*waven/om*real(cmplx(-sinw*dz,dy)*zw2)
frkr1(1,iw)=frkr1(1,iw)+fact1*(0.,2.)*aimag(zw2)*dz*rho
frkr1(2,iw)=frkr1(2,iw)-fact1*(0.,1.)*aimag(zw2)* &
(yof1(k+1)**2-yof1(k)**2+zof1(k+1)**2-zof1(k)**2)*rho
ShallowWater: if (fact2/=0) then
a(j,nof1-ngap+nfs-1+ndof+iw)=a(j,nof1-ngap+nfs-1+ndof+iw) &
+fact2*waven/om*real(cmplx( sinw*dz,dy)*zw2a)
frkr1(1,iw)=frkr1(1,iw)+fact2*(0.,2.)*aimag(zw2a)*dz*rho
frkr1(2,iw)=frkr1(2,iw)-fact2*(0.,1.)*aimag(zw2a)* &
(yof1(k+1)**2-yof1(k)**2+zof1(k+1)**2-zof1(k)**2)*rho
endif ShallowWater
endif !heave or sway/roll
enddo WaveAngles
Columns: do i=ndof-1,nof1-ngap+nfs-1
a(j,i)=pqsn(yof1(k),zof1(k),yof1(k+1),zof1(k+1),yq(i),zq(i),mp)
enddo Columns
enddo BodySegm
!Free-surface condition near to body
FreeSSegmentsNear: do k=nof1,nof1+nf
dy=yof1(k+1)-yof1(k)
AllSources:do i=ndof-1,nof1-ngap+nfs-1 !2-point integration for correct integral at k=i
a(k-ngap,i)=pqsn(yof1(k),zwl,yof1(k+1),zwl,yq(i),zq(i),mp) &
+om**2/g*abs(dy)*(pqs(yof1(k)+0.316*dy,zwl,yq(i),zq(i),mp) &
+pqs(yof1(k)+0.684*dy,zwl,yq(i),zq(i),mp))/2
enddo AllSources
enddo FreeSSegmentsNear
!Free-surface condition far from body (radiation condition)
FarSegments: do k=nof1+nf+1,nof1+nfs-1
dy=yof1(k+1)-yof1(k)
factor=(real(k-nof1-nf)/(nfs-nf-1))**2
AllColumns: do i=ndof-1,nof1-ngap+nfs-1
a(k-ngap,i)=waven*abs(dy)*(pqs(yof1(k)+0.316*dy,zwl,yq(i),zq(i),mp) &
+pqs(yof1(k)+0.684*dy,zwl,yq(i),zq(i),mp))/2 &
-(0.,1.)*(pqs(yof1(k+1)-factor*dy,zwl,yq(i),zq(i),mp) &
-pqs(yof1(k )-factor*dy,zwl,yq(i),zq(i),mp))
enddo AllColumns
enddo FarSegments
if (ndof==1) a(0,0:nof1-ngap+nfs-1)=1. !sum of sources = 0 only in case ndof=1
call simqcd(a(ndof-1,ndof-1),nof1-ngap+nfs+1-ndof,ndof+nmu,nofmax+1,is,1.e-5,detl) !solve equations
if (is.ne.0) call stop1('*** Singular equation system. Coinciding offset points?')
Forces: do m=1,ndof !ndof=1 heave, ndof=2 sway and roll
MotionsAndWaves: do l=1,ndof+nmu !m force index, l motion index
f(m,l)=0
AllBodySegments: do k=1,nof1-1 !k section point index
if (ngap>0.and.any(gap(1:ngap)==k)) cycle
dy=yof1(k+1)-yof1(k)
dz=zof1(k+1)-zof1(k)
if (mp.eq.1) then; fact=+dy !heave motion
else
if (m.eq.1) fact=-dz
if (m.eq.2) fact=+0.5*(yof1(k+1)**2-yof1(k)**2+zof1(k+1)**2-zof1(k)**2)
endif
ForAllSources: do i=ndof-1,nof1-ngap+nfs-1
f(m,l)=f(m,l)+a(i,nof1-ngap+nfs-1+l)*fact* &
(pqs(yof1(k)+0.316*dy,zof1(k)+0.316*dz,yq(i),zq(i),mp) &
+pqs(yof1(k)+0.684*dy,zof1(k)+0.684*dz,yq(i),zq(i),mp))/2
enddo ForAllSources
enddo AllBodySegments
if (l<=ndof) then
addedm1(m,l)=addedm1(m,l)+ciom*f(m,l)/(-om**2)*rho
else
diff1(m,l-ndof)=diff1(m,l-ndof)-ciom*f(m,l)*rho
endif !/2 because of average for nf=25,28; *2 because of symmetry
enddo MotionsAndWaves
enddo Forces
enddo TwoFreeSurfaceDiscretisations
end subroutine addedmassexcitations
subroutine testsuitability(nof1,yof1,zof1)
! Tests whether body sources (locally 1/20 panel length inside section contour) are globally inside
! of contour and outside of inner contour being locally 1/9 panel length inside contour
integer, intent(in):: nof1 !no. of offset points on one section
!integer:: ngap,gap(10) !no. and position of gaps in section contour
integer:: i,j,k !indices
integer:: gap1(10) !no. and point indices of gaps including mirror image
real:: yof1(:),zof1(:) !y, z of offset points, including mirror image
real:: yof(nof1+1),zof(nof1+1) !y, z of points inside section contour
real:: ysave,zsave !offsets of previous contour point
real:: dy,dz !coordinate differences of offset points
real:: yq,zq !y, z of sources
real:: yof2(nofmax),zof2(nofmax) !y, z of original offset points (possibly of half section)
real:: angle !angle under which section contour is seen from a source
real:: ay,az,by,bz !y, z coordinate difference between source and offset point
if(abs(zof1(1)-zof1(nof1))>1e-4)call stop1('First and last z of a section differ')
yof(1:nof1)=yof1(1:nof1); zof(1:nof1)=zof1(1:nof1)
yof(nof1+1)=0.; zof(nof1+1)=-1e5 !additional point high above waterline
TwoTimes: do j=1,2 !Test whether sources inside contour; then interior parallel whether sources outside
Sources: do k=1,nof1-1 !for all source points
!print *,'?',k,ngap,gap(1:ngap),any(gap(1:ngap)==k)
if(ngap>0.and.any(gap(1:ngap)==k)) cycle !not at gaps
yq=(yof1(k)+yof1(k+1))/2-(zof1(k+1)-zof1(k))/20 !generate prelim. source points from original offs.
zq=(zof1(k)+zof1(k+1))/2+(yof1(k+1)-yof1(k))/20
!print *,'yq,zq',k,yq,zq
angle=0
Offsets: do i=1,nof1+1 !compute angle under which the closed section contour is seen from source
if(i==1) then
ay=yq-yof(nof1+1)
az=zq-zof(nof1+1)
else
ay=by
az=bz
endif
by=yq-yof(i)
bz=zq-zof(i)
angle=angle+atan2(ay*bz-az*by,ay*by+az*bz)
!print *,k,i,angle
enddo Offsets
if(abs(angle-merge(6.283,0.,j==1))>1.)call stop1('Unsuitable section discretisation')
enddo Sources
if(j==2)return
do i=1,nof1 !Generate interior contour. Additional high point unchanged.
if(i==1.or.(ngap>0.and.any(gap(1:ngap)==i-1))) then
dy=yof(i+1)-yof(i)
dz=zof(i+1)-zof(i)
elseif(i==nof1.or.(ngap>0.and.any(gap(1:ngap)==i))) then
dy=yof(i)-ysave
dz=zof(i)-zsave
else
dy=(yof(i+1)-ysave)/2
dz=(zof(i+1)-zsave)/2
endif
ysave=yof(i); zsave=zof(i)
yof(i)=yof(i)-dz/9 !interior parallel to contour
zof(i)=zof(i)+dy/9
enddo
!print *,'yof',yof(1:nof1+1); print *,'zof',zof(1:nof1+1)
enddo TwoTimes
end subroutine testsuitability
function which(n,l) !first i for which l(i)=.true., i=1..n. 0 if l(1:n)=.false
integer:: n,which; logical, dimension(n):: l
do which=1,n; if (l(which)) return; enddo
which=0
end function which
function betr2(c)
real:: betr2
complex:: c
betr2=c*conjg(c)
end function betr2
subroutine simqcd(a,n,nr,im,i,s,detl)
! Gauss algorithm for a complex linear equation system. Full matrix, several inhomogeneous vectors.
! a(row,column)=complex coefficient matrix + NEGATIVE right-hand sides as additional columns. When
! finished right-hand sides are replaced by the solutions. Coefficient matrix is destroyed.
! n=no. of equations (rows), nr no. of right-hand sides, im range of first index of a.
! i=0 after normal solutions, /=0 if a pivot is <s (singular or near-singular matrix a).
! detl ist ln(determinant of coefficient matrix). With column pivoting.
real:: s
integer:: n,nr,im,i,j,k,nnr,mmr,k1,k2,imax
complex:: biga,save,detl,a(*)
nnr=(n+nr)*im
mmr=n*im
k1=-im
detl=0.
do i=1,n
k1=k1+im+1
biga=a(k1)
imax=k1
do j=k1+1,k1+n-i
if (abs(biga).lt.abs(a(j))) then
biga=a(j)
imax=j
endif
enddo
if (abs(biga).lt.s) return
detl=detl+log(biga)
a(imax)=a(k1)
do k=k1+im,nnr,im
imax=imax+im
save=-a(imax)/biga
a(imax)=a(k)
a(k)=save
k2=k1
do j=k+1,k+n-i
k2=k2+1
a(j)=a(j)+a(k2)*save
enddo
enddo
enddo
do i=mmr,nnr-1,im
k1=mmr+1
do j=n,2,-1
k1=k1-im
k2=i
save=a(i+j)
do k=k1,k1+j-2
k2=k2+1
a(k2)=a(k2)+a(k)*save
enddo
enddo
enddo
i=0
end subroutine simqcd
function fillcomplmatr(nr,nc,elements) !generates complex matrix from complex elements
integer:: nr,nc !,ir,ic,i
complex:: fillcomplmatr(nr,nc),elements(nr*nc) !elements: 1st row, 2nd row ...
fillcomplmatr=reshape(elements,(/nr,nc/),(/(0.,0.)/),(/2,1/))
!i=0
!do ir=1,nr; do ic=1,nc
! i=i+1
! fillcomplmatr(ir,ic)=elements(i)
!enddo; enddo
end function fillcomplmatr
function fillrealmatr(nr,nc,elements) !generates complex matrix from real elements
integer:: nr,nc !,ir,ic,i
real:: elements(nr*nc) !elements: 1st row, 2nd row ...
complex:: fillrealmatr(nr,nc)
fillrealmatr=reshape(elements,(/nr,nc/),(/0./),(/2,1/))
!i=0
!do ir=1,nr; do ic=1,nc
! i=i+1
! fillrealmatr(ir,ic)=elements(i)
!enddo; enddo
end function fillrealmatr
function prodmatr(ma,mb) !Product of COMPLEX matrices
integer:: nra,nca,nrb,ncb,i,j
complex, intent(in),dimension(:,:):: ma,mb
complex, dimension(size(ma,1),size(mb,2)):: prodmatr
nra=size(ma,1); nca=size(ma,2)
nrb=size(mb,1); ncb=size(mb,2)
if (nca/=nrb) call stop1('Incompatible matrices in prodmatr')
do i=1,nra; do j=1,ncb
prodmatr(i,j)=sum(ma(i,1:nca)*mb(1:nrb,j))
enddo; enddo
end function prodmatr
function vectprod(a,b)
real,intent(in):: a(3),b(3)
real:: vectprod(3)
vectprod(1)=a(2)*b(3)-a(3)*b(2)
vectprod(2)=a(3)*b(1)-a(1)*b(3)
vectprod(3)=a(1)*b(2)-a(2)*b(1)
end function vectprod
function cvectprod(a,b)
complex,intent(in):: a(3),b(3)
complex:: cvectprod(3)
cvectprod(1)=a(2)*b(3)-a(3)*b(2)
cvectprod(2)=a(3)*b(1)-a(1)*b(3)
cvectprod(3)=a(1)*b(2)-a(2)*b(1)
end function cvectprod
subroutine writematr(name,matr) !writes matr below text 'name'
integer:: ir,nr
complex, dimension(:,:):: matr
character(*):: name
nr=size(matr,1)
write(*,*)name
do ir=1,nr
write(*,'(6(2g10.3,1x))')matr(ir,:)
enddo
end subroutine writematr
function logmatr(ma,mb) !log of each element of ma so that |ma(i,j)-mb(i,j)| becomes small
integer:: nra,nca,nrb,ncb,i,j,k !principal value of log if mb(1,1)=999.
real:: ar
complex, intent(in), dimension(:,:):: ma,mb
complex:: logmatr(size(ma,1),size(ma,2))
nra=size(ma,1); nca=size(ma,2)
nrb=size(mb,1); ncb=size(mb,2)
if (mb(1,1)/=(999.,0.).and.(nra/=nrb.or.nca/=ncb)) call stop1('Incompatible matrices in logmatr')
AllElements: do i=1,nra; do j=1,nca
logmatr(i,j)=merge(log(ma(i,j)),(-60.,0.),ma(i,j)/=(0.,0.))
if (mb(1,1)/=(999.,0.)) then !adapt phase angle for small difference to mb
ar=aimag(logmatr(i,j)-mb(i,j))
k=0
do
if (abs(ar)<=3.1416) exit !intentionally a little more than pi
k=k-sign(1.1,ar) !k = no. of required 2pi phase angle corrections
ar=ar-sign(2*pi,ar)
enddo
if (k/=0) logmatr(i,j)=logmatr(i,j)+cmplx(0.,k*2.*pi)
endif
enddo; enddo AllElements
end function logmatr
subroutine signampl !Computes significant amplitudes in natural seaways
logical:: lfirst=.true. !true for first seaway, then false
integer:: i !index
integer:: ifr, nfr !index and number of offsets of frequency valuation curve
integer:: nb !number of points where motions are determined
integer:: iom, nom !index and no. of frequencies in transfer function calc.
integer:: iv,nv !index and number of speeds in transfer function calc.
integer:: imu !index of wave angles in transfer function calc.
integer:: iart,nart !index and no. of transfer functions excl. pressures & drift
integer:: i1 !index of acceleration quantities
integer:: is2 !index of intersections (for sectional forces and moments)
integer:: ib !index of motion points
integer:: ipart,npart !index and no. of sub-intervals within one wave angle interval
integer:: iommax(6*(nsemax+(nbmax*3)/2),nvmax) !frequency index of maximum of transfer function
real:: rla(nomma) !wave length
real:: vf(nvmax) !ship speeds
real:: ybetr2(6*(nsemax+(nbmax*3)/2),nvmax,2*nmumax) !square of response amplitude operators
real:: om(0:nomma+1) !wave circular frequencies
real:: ampken(6*(nsemax+(nbmax*3)/2),nvmax) !significant amplitudes
real:: bewf(30) !acceleration valuation factors
real:: muf(-nmumax:2*nmumax) !wave angles
real:: muin(-nmumax:2*nmumax) !integral of wave energy spreading function over wave angle
real:: fr(10) !frequencies for defining acceleration valuation function
real:: maxdampken(6*(nsemax+(nbmax*3)/2),nvmax) !max. contrib. of one frequency interval to ampken
real:: fxi(nvmax,nmumax) !long. drift force in regular wave per wave amplitude^2
real:: feta(nvmax,nmumax) !transverse drift force in regular wave per wave amplitude^2
real:: barfxi(nvmax) !long. drift force in a seaway
real:: barfeta(nvmax) !transverse drift force in a seaway
real:: dotxdr(nvmax,nmumax) !reduced long. drift velocity in a regular wave/wave ampl.^2
real:: dotydr(nvmax,nmumax) !reduced transv. drift velocity in a regular wave/wave ampl.^2
real:: bardotxdr !reduced longitudinal water drift velocity in a seaway
real:: bardotydr !reduced transverse water drift velocity in a seaway
real:: ypr(nprmax,nsemax),zpr(nprmax,nsemax) !y and z of pressure points
real:: pres2(nprmax,nsemax,nvmax,nmumax) !pressure response amplitude squared
real:: h !significant wave height
real:: t !significant period corr. to c.o.gravity of wave spectrum
real:: mainmu !main direction of the seaway
real:: expon !e in angular spreading of wave energy acc. to cos^e(mu-mainmu)
real:: spuh !peak enhancement factor in Jonswap spectrum
real:: ompeak !circular frequency at peak of wave spectrum
real:: domsp !interval of circular frequency for spectrum integration
real:: sumrm !sum of mu integrals; used to normalize the mu integrals
real:: rmgl,rmbe !left end and range of wave angle interval
real:: omgl,ombe !left end and range of wave circular frequency interval
real:: mu !wave angle
real:: omzw !intermediate circular frequency within one frequency interval
real:: bvh !non-dim. breadth of peak enhancement in Jonswap spectrum
real:: omv !non-dimensional frequency in Jonswap spectrum
real:: sjonsw !Jonswap spectral ordinate
real:: rml !mass of suspended load
real:: skip !used to skip values when reading a file
real:: ome !encounter circular frequency
real:: fq !encounter frequency, ome/(2pi)
real:: fint !quantity for interpolation of acceleration valuation function
real:: bw !valuation factor for accelerations
real:: dampken !contribution to significant amplitude
real:: spi !integral of wave spectrum over frequency
real:: signpres(nprmax,nsemax,nvmax) !significant pressure amplitude
complex:: cpres(nprmax,nsemax) !complex pressure amplitude
complex:: cskip !complex quantity to skip values when reading a file
character(80):: text1 !text describing the case
if (.not.lsign) stop
rewind 20; rewind 21; rewind 24
write(6,'(1x,131("-")/)')
write(6,'(/1x,1a80)')text
write(6,*) 'Significant amplitudes in natural seaways'
open(22,status='scratch',form='unformatted')
write(6,*)'Frequency acceleration valuation factor'
read(5,*)nfr,(fr(ifr),bewf(ifr),ifr=1,nfr)
fr(1)=fr(1)+1.e-8
if (nfr.gt.30) call stop1('Too many frequencies')
write(6,'(f9.3,f20.3)')(fr(ifr),bewf(ifr),ifr=1,nfr)
if (npres>0) read(24,*)((ypr(i,ise),zpr(i,ise),i=1,npres),ise=1,nse)
Seaways: do
read(21,'(a)')text1 !general data read from file 21
if (text1.ne.text) then
write(6,*)text1; write(6,*)text
call stop1('File 21 does not fit to the input data')
endif
read(21,*)nb,ls,g,rho,rml !reads data for nmu,muf complemented to cover [-90 : 270] degrees
read(21,*)nom,(rla(iom),iom=1,nom),nv,(vf(iv),iv=1,nv), &
nmu,(muf(imu),imu=1,nmu),nse,(x(ise),ise=1,nse) !rotations not divided by k as in pdstrip.out
nart=6*(1+nb+merge(1,0,ls)*(nse-1))+2 !+1 for suspended weight rel., + 1 abs. motions
read(5,*)h,t,mainmu,expon,spuh !reads seaway datda
if (h.eq.0) stop
write(6,'(74("-"))')
write(6,'(/ &
&'' Significant wave height '',f10.3/ &
&'' Period T1 corr. to c.o.g. of spectrum '',f10.3/ &
&'' Main wave dir.(0 from stern, 90 from stb.) '',f10.3/ &
&'' Exponent n in cos^n angular spreading funct. '',f10.3/ &
&'' Peak enhancement factor (1 for P.-M.spectrum)'',f10.3)') h,t,mainmu,expon,spuh
! Preparations for integration
ompeak=(4.65+0.182*spuh)/t
write(6,'(/'' Peak period '',f10.3)')6.283/ompeak
domsp=0.025*ompeak
do iom=1,nom
waven=2.*pi/rla(iom)
om(iom)=sqrt(waven*g*tanh(waven*(zbot-zwl)))
enddo
!print *,'om',om(1:nom)
ampken(1:nart+3*nb,1:nv)=0.
maxdampken(1:nart+3*nb,1:nv)=0
if (npres>0) then !if pressure calc.
signpres(1:npres,1:nse,1:nv)=0.
barfxi(1:nv)=0; barfeta(1:nv)=0
bardotxdr=0; bardotydr=0
endif
! Determine integral over mu
sumrm=0
WaveAngles: do imu=1,nmu
rmgl=(merge(muf(imu-1),muf(nmu)-2*pi,imu>1)+muf(imu))/2. !left end of interval
rmbe=(merge(muf(imu+1),muf(1)+2*pi,imu<nmu)+muf(imu))/2.-rmgl !range of interval
if (rmbe>0.78.and.abs(rmgl-mainmu*pi/180)<0.5*pi) call stop1( &
'The specified wave angles are not sufficient for this seaway')
npart=rmbe/0.02+1
muin(imu)=0
Parts: do ipart=1,npart
mu=rmgl+(ipart-0.5)*rmbe/npart
if ( modulo(mainmu*pi/180-mu,2*pi).lt.pi/2. &
.or.modulo(mainmu*pi/180-mu,2*pi).gt.pi*3/2.) &
muin(imu)=muin(imu)+(cos(mainmu*pi/180-mu))**expon*rmbe/npart
enddo Parts
sumrm=sumrm+muin(imu)
enddo WaveAngles
muin(1:nmu)=muin(1:nmu)/sumrm
!print *,'muin',muin(1:nmu)
! Integration over omega for all responses
om(0)=max(5.*ompeak-om(1),om(1))
om(nom+1)=min(1.2*ompeak-om(nom),om(nom))
Frequencies: do iom=1,nom
omgl=(om(iom+1)+om(iom))/2 !left end of interval
ombe=(om(iom-1)+om(iom))/2.-omgl !range of interval
npart=ombe/domsp+1 !no. of sub-intervals
spi=0 !spectrum part integral
SubIntervals: do ipart=1,npart !integration of wave spectrum from omgl to omgl+ombe
omzw=omgl+(ipart-0.5)*ombe/npart
bvh=merge(0.09,0.07,omzw>ompeak)
omv=omzw/ompeak
sjonsw=h**2*t*(177.5-6.52*spuh)/(t*omzw)**5*exp(-1.25/omv**4)
if (abs(omv-1.).lt.0.36) sjonsw=sjonsw*spuh**exp(-(omv-1)**2/(2.*bvh**2))
spi=spi+sjonsw*ombe/npart !=wave spectrum, integral over current omega interval
enddo SubIntervals
Speeds: do iv=1,nv
Angles: do imu=1,nmu
if (lfirst) then !for first seaway: read transfer functions for all angles
read(21,*)(ybetr2(iart,iv,imu),iart=1,nart)
write(22)(ybetr2(iart,iv,imu),iart=1,nart) !write into unformatted scratch file
if (npres>0) then
read(24,*)((skip,i=1,4),(cpres(i,ise),cskip,cskip,i=1,npres),ise=1,nse)
pres2(1:npres,1:nse,iv,imu)=abs(cpres(1:npres,1:nse))**2
read(24,*)fxi(iv,imu),feta(iv,imu),dotxdr(1,imu),dotydr(iv,imu)
write(22)pres2(1:npres,1:nse,iv,imu),fxi(iv,imu),feta(iv,imu),dotxdr(1,imu),dotydr(iv,imu)
endif
else !for further seaways: read transfer functions from unformatted scratch file (faster)