-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathwwm_gridcf.F90
675 lines (624 loc) · 25.9 KB
/
wwm_gridcf.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
#include "wwm_functions.h"
!#define DEBUG
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE SPHERICAL_COORDINATE_DISTANCE(LON1, LON2, LAT1, LAT2, DIST)
USE DATAPOOL
IMPLICIT NONE
REAL(rkind), intent(in) :: LON1, LON2, LAT1, LAT2
REAL(rkind), intent(out) :: DIST
REAL(rkind) :: lon1rad, lon2rad, lat1rad, lat2rad
REAL(rkind) :: COEFF
REAL(rkind) :: x1, y1, z1, x2, y2, z2
REAL(rkind) :: scalprod
COEFF=PI/180.0_rkind
lon1rad=LON1*COEFF
lon2rad=LON2*COEFF
lat1rad=LAT1*COEFF
lat2rad=LAT2*COEFF
x1=cos(lon1rad)*cos(lat1rad)
y1=sin(lon1rad)*cos(lat1rad)
z1=sin(lat1rad)
x2=cos(lon2rad)*cos(lat2rad)
y2=sin(lon2rad)*cos(lat2rad)
z2=sin(lat2rad)
scalprod=x1*x2+y1*y2+z1*z2;
IF (scalprod .ge. 1) THEN
DIST=0;
ELSE
DIST=acos(scalprod)
END IF
END SUBROUTINE
!********************************************************************************
!* See formula for area of spherical triangle in *
!* http://math.stackexchange.com/questions/9819/area-of-a-spherical-triangle *
!* *
!********************************************************************************
SUBROUTINE SPHERICAL_COORDINATE_AREA(LON1, LON2, LON3, LAT1, LAT2, LAT3, AREA)
USE DATAPOOL
IMPLICIT NONE
REAL(rkind), intent(in) :: LON1, LON2, LON3, LAT1, LAT2, LAT3
REAL(rkind), intent(out) :: AREA
REAL(rkind) :: DistA, DistB, DistC, DistS
REAL(rkind) :: eTan1, eTan2, eTan3, eTan4
REAL(rkind) :: eProd, sqrtProd
CALL SPHERICAL_COORDINATE_DISTANCE(LON1, LON2, LAT1, LAT2, DistA)
CALL SPHERICAL_COORDINATE_DISTANCE(LON1, LON3, LAT1, LAT3, DistB)
CALL SPHERICAL_COORDINATE_DISTANCE(LON2, LON3, LAT2, LAT3, DistC)
DistS=(DistA + DistB + DistC)/2.0_rkind
eTan1=tan(DistS/2.0_rkind)
eTan2=tan((DistS - DistA)/2.0_rkind)
eTan3=tan((DistS - DistB)/2.0_rkind)
eTan4=tan((DistS - DistC)/2.0_rkind)
eProd=eTan1*eTan2*eTan3*eTan4
sqrtProd=SQRT(eProd)
AREA=4.0_rkind*ATAN(sqrtProd)
END SUBROUTINE
SUBROUTINE CORRECT_SINGLE_DXP(DXP)
USE DATAPOOL
IMPLICIT NONE
real(rkind), intent(inout) :: DXP
IF (DXP .le. -180_rkind) THEN
DXP=DXP + 360
END IF
IF (DXP .ge. 180) THEN
DXP=DXP - 360
END IF
END SUBROUTINE
SUBROUTINE INIT_SPATIAL_GRID
USE DATAPOOL
IMPLICIT NONE
REAL(rkind) :: TL1, TL2, TL3
REAL(rkind) :: TMPTLMIN
REAL(rkind) :: TMPTLMAX
REAL(rkind) :: DBLTMP, DXP1, DXP2, DXP3, DYP1, DYP2, DYP3
REAL(rkind) :: P1_XLOC, P2_XLOC, P3_XLOC
REAL(rkind) :: PROV1, PROV2, PROV3
REAL(rkind) :: AREA, AREA_RAD
INTEGER :: I1, I2, I3, TMPINE, NI(3)
INTEGER :: J1, J2, J3
INTEGER :: IP, IE, IEWRONG, IEWRONGSUM
LOGICAL :: LWRONG
#ifdef MPI_PARALL_GRID
REAL(rkind) :: AVETA_GL, TLMIN_GL, TLMAX_GL
AVETA_GL = ZERO; TLMIN_GL = ZERO; TLMAX_GL = ZERO
#endif
TMPTLMIN = 10.E10_rkind
TMPTLMAX = ZERO
TLMIN = VERYLARGE
TLMAX = VERYSMALL
AVETL = ZERO
AVETA = ZERO
LWRONG = .FALSE.
SELECT CASE (DIMMODE)
!
! *** One dimension mode
!
CASE (1)
IF (LVAR1D) THEN
DX1(0) = XP(2)- XP(1)
DX1(1) = DX1(0)
DX1(MNP) = XP(MNP) - XP(MNP-1)
DX1(MNP+1) = DX1(MNP)
DX2(0) = DX1(0)
DX2(MNP+1) = DX1(MNP)
DO IP = 2, MNP-1 ! Bandwith at gridpoints
DX1(IP) = (XP(IP)-XP(IP-1))/2. + (XP(IP+1)-XP(IP))/2.
END DO
DO IP = 2, MNP ! Stepwidth between gridpoints K and K-1
DX2(IP) = XP(IP) - XP(IP-1)
END DO
DX2(1) = DX1(0)
END IF
TL1 = 0.
DO IE = 1, MNP-1
TL1 = TL1 + ABS(XP(IE+1) - XP(IE))
END DO
AVETL = TL1/MyREAL(MNP)
!
! *** Two dimension mode
!
CASE(2)
IEWRONGSUM = 0
DO IE = 1, MNE
I1 = INE(1,IE)
I2 = INE(2,IE)
I3 = INE(3,IE)
NI = INE(:,IE)
! Special treatment for elements overpassing the dateline
P1_XLOC = XP(I1);
P2_XLOC = XP(I2);
P3_XLOC = XP(I3);
IF (LSPHE) THEN
! Kevin Martins' fix to close the mesh at the
! dateline: understand whether an element has a single node
! isolated on the other side, and bring it
! back. This approach does nothing with the element
! that contains the pole, and apparently is stable
! there.
IF ( (P1_XLOC - P2_XLOC).GT.180 .AND. &
& (P3_XLOC - P2_XLOC).GT.180 ) THEN
! In this case, P2 is 'isolated' to the East of the dateline; we "bring it back" towards the West
P2_XLOC = 180.0_rkind + ABS(-180.0_rkind - P2_XLOC)
ELSE IF ( (P2_XLOC - P1_XLOC).GT.180 .AND. &
& (P3_XLOC - P1_XLOC).GT.180 ) THEN
! In this case, P1 is 'isolated' to the East of the dateline; we "bring it back" towards the West
P1_XLOC = 180.0_rkind + ABS(-180.0_rkind - P1_XLOC)
ELSE IF ( (P1_XLOC - P3_XLOC).GT.180 .AND. &
& (P2_XLOC - P3_XLOC).GT.180 ) THEN
! In this case, P3 is 'isolated' to the East of the dateline; we "bring it back" towards the West
P3_XLOC = 180.0_rkind + ABS(-180.0_rkind - P3_XLOC)
ELSE IF ( (P1_XLOC - P2_XLOC).GT.180 .AND. &
& (P1_XLOC - P3_XLOC).GT.180 ) THEN
! In this case, P1 is 'isolated' to the West of the dateline; we "bring it back" towards the East
P1_XLOC = -180.0_rkind - ABS(180.0_rkind - P1_XLOC)
ELSE IF ( (P2_XLOC - P1_XLOC).GT.180 .AND. &
& (P2_XLOC - P3_XLOC).GT.180 ) THEN
! In this case, P2 is 'isolated' to the West of the dateline; we "bring it back" towards the East
P2_XLOC = -180.0_rkind - ABS(180.0_rkind - P2_XLOC)
ELSE IF ( (P3_XLOC - P1_XLOC).GT.180 .AND. &
& (P3_XLOC - P2_XLOC).GT.180 ) THEN
! In this case, P3 is 'isolated' to the West of the dateline; we "bring it back" towards the East
P3_XLOC = -180.0_rkind - ABS(180.0_rkind - P3_XLOC)
END IF
END IF
IF (IGRIDTYPE.ne.2) THEN
DXP1 = P2_XLOC - P1_XLOC
DYP1 = YP(I2) - YP(I1)
DXP2 = P3_XLOC - P2_XLOC
DYP2 = YP(I3) - YP(I2)
DXP3 = P1_XLOC - P3_XLOC
DYP3 = YP(I1) - YP(I3)
IF (APPLY_DXP_CORR) THEN
! the option APPLY_DXP_CORR does something
! similar to Kevin Martin's approach to close the
! mesh at the dateline. If the longitudinal
! distance DXP is >180deg or <-180deg it is
! adjusted subtracting or adding 360deg. This
! approach also acts on the element containing
! the pole, correcting a single longitudinal
! disnance of the three, making this element
! unstable.
! This option could be a substitute of Kevin
! Martin's fix if a special treatment of the
! element containing the pole is introduced.
! This option could be removed if it is a
! duplication of Kevin Martin's fix.
CALL CORRECT_SINGLE_DXP(DXP1)
CALL CORRECT_SINGLE_DXP(DXP2)
CALL CORRECT_SINGLE_DXP(DXP3)
END IF
IEN(1,IE) = - DYP2
IEN(2,IE) = DXP2
IEN(3,IE) = - DYP3
IEN(4,IE) = DXP3
IEN(5,IE) = - DYP1
IEN(6,IE) = DXP1
DBLTMP = (DXP3*DYP1 - DYP3*DXP1)*ONEHALF
IF (LSPHE .and. USE_EXACT_FORMULA_SPHERICAL_AREA) THEN
CALL SPHERICAL_COORDINATE_AREA(XP(I1), XP(I2), XP(I3), YP(I1), YP(I2), YP(I3), AREA_RAD)
AREA=AREA_RAD*RADDEG*RADDEG
#ifdef DEBUG
WRITE(STAT%FHNDL,*) 'IE=', IE
WRITE(STAT%FHNDL,*) 'I123=', I1, I2, I3
WRITE(STAT%FHNDL,*) 'XP123=', XP(I1), XP(I2), XP(I3)
WRITE(STAT%FHNDL,*) 'YP123=', YP(I1), YP(I2), YP(I3)
WRITE(STAT%FHNDL,*) ' AREA=', AREA
WRITE(STAT%FHNDL,*) 'DBLTMP=', DBLTMP
FLUSH(STAT%FHNDL)
#endif
ELSE
AREA=DBLTMP
END IF
TRIA(IE) = AREA
END IF
IF (TRIA(IE) .LT. 0.0) THEN
TMPINE = INE(2,IE)
INE(2,IE) = INE(3,IE)
INE(3,IE) = TMPINE
IF (ANY(INE(:,IE) .GT. MNP)) CALL WWM_ABORT('WRITE ELEMENT CONNECTION TALBE HAS NODENUMBERS GT MNP')
I2 = INE(2,IE)
I3 = INE(3,IE)
TRIA(IE) = -TRIA(IE)
PROV1=IEN(6,IE) ! DXP1
PROV2=IEN(2,IE) ! DXP2
PROV3=IEN(4,IE) ! DXP3
IEN(6,IE)=-PROV3
IEN(2,IE)=-PROV2
IEN(4,IE)=-PROV1
PROV1= - IEN(5,IE)
PROV2= - IEN(1,IE)
PROV3= - IEN(3,IE)
IEN(1,IE) = PROV2
IEN(3,IE) = PROV1
IEN(5,IE) = PROV3
LWRONG = .TRUE.
IEWRONG = IE
IEWRONGSUM = IEWRONGSUM + 1
WRITE(DBG%FHNDL,*) 'WRONG ELEMENT', IE
J1=INE(1,IE)
J2=INE(2,IE)
J3=INE(3,IE)
WRITE(DBG%FHNDL,*) 'NODENUMBERS I1=', J1
WRITE(DBG%FHNDL,*) 'NODENUMBERS I2=', J2
WRITE(DBG%FHNDL,*) 'NODENUMBERS I3=', J3
WRITE(DBG%FHNDL,*) 'XP1, YP1=', XP(J1), YP(J1)
WRITE(DBG%FHNDL,*) 'XP2, YP2=', XP(J2), YP(J2)
WRITE(DBG%FHNDL,*) 'XP3, YP3=', XP(J3), YP(J3)
WRITE(DBG%FHNDL,*) 'DXP1, DYP1=', DXP1, DYP1
WRITE(DBG%FHNDL,*) 'DXP2, DYP2=', DXP2, DYP2
WRITE(DBG%FHNDL,*) 'DXP3, DYP3=', DXP3, DYP3
WRITE(DBG%FHNDL,'(A40,6F15.8)') 'EDGELENGTHS OF THE WRONG ELEMENT', IEN(:,IE)
ELSE IF (TRIA(IE) .LT. THR) THEN
write(DBG%FHNDL,*) 'IE=', IE, ' TRIA=', TRIA(IE)
write(DBG%FHNDL,*) 'DXP1=', DXP1, ' DXP3=', DXP3
write(DBG%FHNDL,*) 'DYP1=', DYP1, ' DYP3=', DYP3
write(DBG%FHNDL,*) 'I123=', I1, I2, I3
write(DBG%FHNDL,*) 'XP,YP(I1)=', XP(I1), YP(I1)
write(DBG%FHNDL,*) 'XP,YP(I2)=', XP(I2), YP(I2)
write(DBG%FHNDL,*) 'XP,YP(I3)=', XP(I3), YP(I3)
CALL WWM_ABORT('too small triangles')
END IF
TL1 = SQRT(IEN(5,IE)**2 + IEN(6,IE)**2)
TL2 = SQRT(IEN(3,IE)**2 + IEN(4,IE)**2)
TL3 = SQRT(IEN(1,IE)**2 + IEN(2,IE)**2)
TMPTLMIN = MIN(TL1, TL2, TL3)
TMPTLMAX = MAX(TL1, TL2, TL3)
IF (TLMIN > TMPTLMIN) THEN
TLMIN = TMPTLMIN
END IF
IF (TLMAX < TMPTLMAX) THEN
TLMAX = TMPTLMAX
END IF
AVETA = AVETA+TRIA(IE)
END DO
IF (LGSE) THEN
DO IE = 1, MNE
IEND(1,1,IE) = DOT_PRODUCT(IEN(1:2,IE),IEN(1:2,IE))
IEND(1,2,IE) = DOT_PRODUCT(IEN(1:2,IE),IEN(3:4,IE))
IEND(1,3,IE) = DOT_PRODUCT(IEN(1:2,IE),IEN(5:6,IE))
IEND(2,1,IE) = DOT_PRODUCT(IEN(3:4,IE),IEN(1:2,IE))
IEND(2,2,IE) = DOT_PRODUCT(IEN(3:4,IE),IEN(3:4,IE))
IEND(2,3,IE) = DOT_PRODUCT(IEN(3:4,IE),IEN(5:6,IE))
IEND(3,1,IE) = DOT_PRODUCT(IEN(5:6,IE),IEN(1:2,IE))
IEND(3,2,IE) = DOT_PRODUCT(IEN(5:6,IE),IEN(3:4,IE))
IEND(3,3,IE) = DOT_PRODUCT(IEN(5:6,IE),IEN(5:6,IE))
END DO
ENDIF
#ifdef MPI_PARALL_GRID
CALL MPI_ALLREDUCE(TLMIN,TLMIN_GL,1,rtype,MPI_MIN,comm,ierr)
CALL MPI_ALLREDUCE(TLMAX,TLMAX_GL,1,rtype,MPI_MAX,comm,ierr)
CALL MPI_ALLREDUCE(AVETA,AVETA_GL,1,rtype,MPI_SUM,comm,ierr)
AVETA = AVETA_GL / MyREAL(ne_global)
AVETL = (TLMIN_GL+TLMAX_GL)/TWO
#else
AVETA = AVETA/MyREAL(MNE)
AVETL = (TLMIN+TLMAX)/TWO
#endif
#ifdef DEBUG
WRITE(STAT%FHNDL,*) 'AVETL=', AVETL
WRITE(STAT%FHNDL,*) 'TLMIN=', TLMIN
WRITE(STAT%FHNDL,*) 'TLMAX=', TLMAX
FLUSH(STAT%FHNDL)
#endif
IF (LWRONG) THEN
#ifdef MPI_PARALL_GRID
IF (myrank == 0) THEN
#endif
GRDCOR%FNAME = 'sysrenum.dat'
OPEN(GRDCOR%FHNDL, FILE=GRDCOR%FNAME, STATUS='UNKNOWN')
WRITE(GRDCOR%FHNDL,'(I10)') 0
WRITE(GRDCOR%FHNDL,'(I10)') MNP
IF (LSPHE) THEN
DO IP = 1, MNP
WRITE(GRDCOR%FHNDL,'(I10,3F15.8)') IP-1, XP(IP), YP(IP), DEP(IP)
END DO
ELSE
DO IP = 1, MNP
WRITE(GRDCOR%FHNDL,'(I10,3F15.6)') IP-1, XP(IP), YP(IP), DEP(IP)
END DO
ENDIF
WRITE(GRDCOR%FHNDL,'(I10)') MNE
DO IE = 1, MNE
WRITE(GRDCOR%FHNDL,'(5I10)') INE(1,IE)-1, INE(2,IE)-1, INE(3,IE)-1, 0, IE-1
END DO
CLOSE(GRDCOR%FHNDL)
#ifdef MPI_PARALL_GRID
END IF
#endif
WRITE(DBG%FHNDL,*) 'The Elements in your mesh are not correctly numbered!'
WRITE(DBG%FHNDL,*) 'New mesh is written to', TRIM(GRDCOR%FNAME)
WRITE(DBG%FHNDL,*) 'There are totally', IEWRONGSUM, 'wrong Elements'
WRITE(DBG%FHNDL,*) 'The last wrong element has the number', IEWRONG
CALL WWM_ABORT('SPATIAL GRID - ELEMENTS NOT CORRECTLY NUMBERED')
END IF
CASE DEFAULT
CALL WWM_ABORT('SPATIAL GRID - WRONG CASE - MUST BE 1D or 2D')
END SELECT
IF (LSPHE) THEN
AVETL = AVETL*REARTH*PI/180.0
DO IP = 1, MNP
INVSPHTRANS(IP,1) = ONE/(DEGRAD*REARTH*COS(YP(IP)*DEGRAD))
INVSPHTRANS(IP,2) = ONE/(DEGRAD*REARTH)
IF ( ABS(COS(YP(IP))) .LT. THR ) THEN
CALL WWM_ABORT('SPATIAL GRID INVSPHETRANS IS SINGULAR')
END IF
END DO
ELSE
INVSPHTRANS = 1.
END IF
IF (LTEST) THEN
SELECT CASE (DIMMODE)
CASE (1)
WRITE(STAT%FHNDL,101) AVETL, TLMIN, TLMAX
CASE (2)
WRITE(STAT%FHNDL,102) AVETA, AVETL, TLMIN, TLMAX
IF (ITEST > 100) THEN
WRITE(STAT%FHNDL,*) ' The element area = '
DO IE = 1, MNE
WRITE(STAT%FHNDL,*) 'IE = ', IE, TRIA(IE)
END DO
END IF
CASE DEFAULT
END SELECT
END IF
101 FORMAT (1X,'The averege element length = ',F16.5/ &
& 1X,'The minimum element length = ',F16.5/ &
& 1X,'The maximum element length = ',F16.5/ )
102 FORMAT (1X,'The averege element area = ',F16.5/ &
& 1X,'The averege element length = ',F16.5/ &
& 1X,'The minimum element length = ',F16.5/ &
& 1X,'The maximum element length = ',F16.5/ )
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE INIT_TRIAD
USE DATAPOOL
IMPLICIT NONE
TRI_ISP = INT( LOG(TWO) / XISLN )
TRI_ISP1 = TRI_ISP + 1
TRI_WISP = (2. - XIS**TRI_ISP) / (XIS**TRI_ISP1 - XIS**TRI_ISP)
TRI_WISP1 = 1. - TRI_WISP
TRI_ISM = INT( LOG(0.5) / XISLN )
TRI_ISM1 = TRI_ISM - 1
TRI_WISM = (XIS**TRI_ISM -0.5) / (XIS**TRI_ISM - XIS**TRI_ISM1)
TRI_WISM1 = 1. - TRI_WISM
TRI_ISBEGIN = MAX(1, 1-TRI_ISM1)
IF (MESTR .eq. 1) THEN
TRI_ARR(1) = 0.1
TRI_ARR(2) = 2.2
TRI_ARR(3) = 10.
TRI_ARR(4) = 0.2
TRI_ARR(5) = 0.01
IF (TRICO .GT. 0.) TRI_ARR(1) = TRICO
IF (TRIRA .GT. 0.) TRI_ARR(2) = TRIRA
IF (TRIURS .GT. 0.) TRI_ARR(5) = TRIURS
END IF
IF (MESTR .eq. 5) THEN
TRI_ARR(1) = 0.25
TRI_ARR(2) = 2.5
TRI_ARR(3) = 10.
TRI_ARR(4) = 0.2
TRI_ARR(5) = 0.01
IF (TRICO .GT. 0.) TRI_ARR(1) = TRICO
IF (TRIRA .GT. 0.) TRI_ARR(2) = TRIRA
IF (TRIURS .GT. 0.) TRI_ARR(5) = TRIURS
END IF
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE INIT_SPECTRAL_GRID()
USE DATAPOOL
IMPLICIT NONE
INTEGER :: IS, ID
INTEGER :: NUMSIG1, NUMSIG2
REAL(rkind) :: SSB, SPECTRAL_BANDWIDTH
REAL(rkind) :: TMP, CO1
ALLOCATE( SPSIG(NUMSIG), SPDIR(NUMDIR), FR(NUMSIG), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_initio, allocate error 5')
SPSIG = zero
SPDIR = zero
FR = zero
ALLOCATE( COSTH(NUMDIR), SINTH(NUMDIR), COS2TH(NUMDIR), SIN2TH(NUMDIR), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_initio, allocate error 6')
COSTH = zero
SINTH = zero
COS2TH = zero
SIN2TH = zero
ALLOCATE( SINCOSTH(NUMDIR), SIGPOW(NUMSIG,6), DS_BAND(0:NUMSIG+1), DS_INCR(0:NUMSIG+1), ID_NEXT(NUMDIR), ID_PREV(NUMDIR), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_initio, allocate error 7')
SINCOSTH = zero
SIGPOW = zero
DS_BAND = zero
DS_INCR = zero
DO ID=1,NUMDIR-1
ID_NEXT(ID)=ID+1
END DO
ID_NEXT(NUMDIR)=1
DO ID=2,NUMDIR
ID_PREV(ID)=ID-1
END DO
ID_PREV(1)=NUMDIR
SGLOW = PI2*FRLOW
SGHIGH = PI2*FRHIGH
!2do check FRINTF for LOPTSIG
IF (LOPTSIG) THEN
SFAC = 1.1_rkind
FRINTF = LOG(SFAC)
FRHIGH = FRLOW
DO IS=2,NUMSIG
FRHIGH = FRHIGH * SFAC
END DO
ELSE
FRINTF = LOG(SGHIGH/SGLOW)/MyREAL(NUMSIG-1)
SFAC = EXP(FRINTF)
END IF
SGLOW = PI2*FRLOW
SGHIGH = PI2*FRHIGH
FRATIO = SFAC
FRINTH = SQRT(SFAC)
FR(1) = FRLOW
WRITE(STAT%FHNDL,*) 'RESOLUTION IN SIGMA SPACE AND FACTORS'
WRITE(STAT%FHNDL,*) 'SGLOW', SGLOW
WRITE(STAT%FHNDL,*) 'LOPTSIG', LOPTSIG
WRITE(STAT%FHNDL,*) 'SFAC, FRINTF, FRINTH', SFAC, FRINTF, FRINTH
DO IS = 2, NUMSIG
FR(IS) = FR(IS-1) * SFAC
END DO
SPSIG = FR * PI2
WRITE(STAT%FHNDL,*) 'SFAC=', SFAC
WRITE(STAT%FHNDL,'("+TRACE...",A,F15.4)') 'REL. FREQ. Distribution is =', FRINTF
IF ( ABS(FRINTF - .1)/FRINTF * 100. .GT. 1. ) THEN
WRITE(DBG%FHNDL,*) 'Freq. resolution is not optimal for Snl4'
WRITE(DBG%FHNDL,'(3F15.4)') 1. + FRINTF, ABS(FRINTF - .1)/FRINTF * 100.
WRITE(DBG%FHNDL,*) 'rel. freq. res. should be 1.1 is now', 1. + FRINTF, 'ERROR IS:', ABS(FRINTF - .1)/FRINTF * 100.
END IF
IF (NUMSIG .GE. 2) THEN
DS_BAND(0) = SPSIG(2)- SPSIG(1)
DS_BAND(1) = DS_BAND(0)
DS_BAND(NUMSIG) = SPSIG(NUMSIG) - SPSIG(NUMSIG-1)
DS_BAND(NUMSIG+1) = DS_BAND(NUMSIG)
DS_INCR(0) = DS_BAND(0)
DS_INCR(1) = DS_BAND(0)
DS_INCR(NUMSIG) = DS_BAND(NUMSIG)
DS_INCR(NUMSIG+1) = DS_INCR(NUMSIG)
DO IS = 2, NUMSIG-1 ! Bandwith at gridpoints
DS_BAND(IS) = (SPSIG(IS)-SPSIG(IS-1))/2. + (SPSIG(IS+1)-SPSIG(IS))/2.
END DO
DO IS = 2, NUMSIG ! Stepwidth between gridpoints K and K-1
DS_INCR(IS) = SPSIG(IS) - SPSIG(IS-1)
END DO
END IF
SPECTRAL_BANDWIDTH = (FRHIGH-FRLOW)
SSB = (SPSIG(2)-SPSIG(1))/PI2
NUMSIGL = INT(SPECTRAL_BANDWIDTH/SSB)
ALLOCATE(SPSIGL(NUMSIGL));SPSIGL = ZERO
SPSIGL(1) = FRLOW * PI2
DO IS = 2, NUMSIGL
SPSIGL(IS) = SPSIGL(IS-1) + SSB
ENDDO
!
! *** the ratio of the consecutive frequency ... for quad
!
IF ( NUMSIG .GT. 3) THEN
NUMSIG2 = INT(MyREAL(NUMSIG)/TWO)
NUMSIG1 = NUMSIG2-1
XIS = SPSIG(NUMSIG2)/SPSIG(NUMSIG1)
XISLN = LOG(XIS)
IF (SMETHOD .GT. 0 .AND. MESTR .GT. 0) CALL INIT_TRIAD
ELSE
IF (SMETHOD .GT. 0 .AND. MESNL .GT. 0) CALL WWM_ABORT('TOO LESS FREQ FOR SNL4 SET MESNL = 0')
IF (SMETHOD .GT. 0 .AND. MESTR .GT. 0) CALL WWM_ABORT('TOO LESS FREQ FOR SNL3 SET MESTR = 0')
END IF
!
! *** frequency grid in [Hz]
!
FR = SPSIG/PI2
!
! *** set the distribution of the dectional domain
!
IF (NUMDIR == 1) THEN
DDIR = 1._rkind
SPDIR(1) = 0._rkind
ELSE
IF (LCIRD) THEN
DDIR = ABS(MAXDIR-MINDIR)/MyREAL(NUMDIR)
DO ID = 1, NUMDIR
IF (LSTAG) THEN
SPDIR(ID) = MINDIR + DDIR * MyREAL(ID-1) + DDIR/2.0
ELSE
SPDIR(ID) = MINDIR + DDIR * MyREAL(ID-1)
END IF
IF (SPDIR(ID) >= PI2) SPDIR(ID) = SPDIR(ID) - PI2
END DO
ELSE
IF (LNAUTIN) THEN
TMP = MAXDIR
MAXDIR = MINDIR
MINDIR = TMP
END IF
WRITE(STAT%FHNDL,*) 'MINDIR MAXDIR', MINDIR, MAXDIR, MINDIR*RADDEG, MAXDIR*RADDEG
IF (MAXDIR.LT.MINDIR) MAXDIR = MAXDIR + PI2
DDIR = (MAXDIR-MINDIR) / MyREAL(NUMDIR)
DO ID = 1, NUMDIR
SPDIR(ID) = MINDIR + DDIR * MyREAL(ID-1)
END DO
END IF
END IF
!
! *** set trig. in angular space
!
COSTH(:) = COS(SPDIR(:))
SINTH(:) = SIN(SPDIR(:))
COS2TH(:) = COS(SPDIR(:))**2
SIN2TH(:) = SIN(SPDIR(:))**2
SINCOSTH(:) = COS(SPDIR(:))*SIN(SPDIR(:))
!
! *** set POWERS OF SPSIG
!
SIGPOW(:,1) = SPSIG(:)
SIGPOW(:,2) = SPSIG(:)**2
SIGPOW(:,3) = SPSIG(:) * SIGPOW(:,2)
SIGPOW(:,4) = SPSIG(:) * SIGPOW(:,3)
SIGPOW(:,5) = SPSIG(:) * SIGPOW(:,4)
SIGPOW(:,6) = SPSIG(:) * SIGPOW(:,5)
!
FDIR = FRINTF * DDIR
DELTH = PI2/MyREAL(NUMDIR)
!
! WAM related definitions
!
ALLOCATE(DFIM(NUMSIG), DFFR(NUMSIG), DFFR2(NUMSIG), DFIMOFR(NUMSIG), FR5(NUMSIG), FRM5(NUMSIG), COFRM4(NUMSIG), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_gridcf, allocate error 1')
DFIM = ZERO
DFFR = ZERO
DFFR2 = ZERO
DFIMOFR = ZERO
FR5 = ZERO
FRM5 = ZERO
COFRM4 = ZERO
IF (ISOURCE == 2) THEN
ALLOCATE(TH(NUMDIR), stat=istat)
th = zero
DELTH = PI2/REAL(NUMDIR)
DELTR = DELTH*REARTH
DO ID=1,NUMDIR
IF (LSTAG) THEN
TH(ID) = REAL(ID-1)*DELTH
ELSE
TH(ID) = REAL(ID-1)*DELTH + 0.5*DELTH
ENDIF
COSTH(ID) = COS(TH(ID))
SINTH(ID) = SIN(TH(ID))
ENDDO
ENDIF
CO1 = 0.5*(FRATIO-1.)*DELTH
DFIM(1)= CO1*FR(1)
DO IS=2,NUMSIG-1
DFIM(IS)=CO1 * (FR(IS)+FR(IS-1))
ENDDO
DFIM(IS)=CO1*FR(IS-1)
DO IS = 1, NUMSIG
DFFR(IS) = DFIM(IS)*FR(IS)
DFFR2(IS) = DFIM(IS)*FR(IS)**2
DFIMOFR(IS) = DFIM(IS)/FR(IS)
FR5(IS) = FR(IS)**5
FRM5(IS) = ONE/FR5(IS)
COFRM4(IS) = COEF4*G9/FR(IS)**4
END DO
FLOGSPRDM1=1./LOG10(FRATIO)
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE DEALLOC_SPECTRAL_GRID()
USE DATAPOOL
IMPLICIT NONE
DEALLOCATE(DFIM, DFFR, DFFR2, DFIMOFR, FR5)
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************