subroutine waveh(ia1, 1,1WAVEH.52
*CALL ARGWVAL 
WAVEH.53
*CALL ARGWVFD WAVEH.54
& l_wvtra,ndata,ngrid,nfldmax,len1,energy,windsp,windir,wh,icode) WAVEH.55
C WAVEH.56
WAVEH.57
*CALL TYPWVAL WAVEH.58
*CALL TYPWVFD WAVEH.59
WAVEH.60
C WAVEH.61
c * set max number of wavetrains to search for * WAVEH.62
parameter(kwtmax=4) WAVEH.63
WAVEH.64
C UM constants WAVEH.65
WAVEH.66
*CALL C_MDI WAVEH.67
*CALL C_G WAVEH.68
*CALL C_PI WAVEH.69
WAVEH.70
INTEGER ICODE WAVEH.71
real wh(ngrid,nfldmax) ! array for output fields ! OUT WAVEH.72
real energy(len1),windir(ndata),windsp(ndata) ! IN WAVEH.73
C WAVEH.74
C # LOCAL arrays for use by wavetrain partitioning WAVEH.75
C WAVEH.76
CC note can rearange and not declare these arrays here by calling WAVEH.77
CC with array WH - have kwtmax consec wavehts then kwtmax periods etc. WAVEH.78
CC real pswh(ndata,kwtmax) ! height of wavetrains WAVEH.79
CC real perio(ndata,kwtmax) ! period of wavetrain WAVEH.80
CC real pdir(ndata,kwtmax) ! direction of wavetrain WAVEH.81
integer kwtot(ndata) ! number of wavetrains present WAVEH.82
C WAVEH.83
C LOCAL ARRAYS WAVEH.84
C WAVEH.85
REAL ECOS(ndata),ESIN(ndata),EMAX(ndata), WAVEH.86
- ESIN1(ndata),ECOS1(ndata),EMSW(ndata),EMWS(ndata), WAVEH.87
- ESIN2(ndata),ECOS2(ndata),EWS(ndata) WAVEH.88
WAVEH.89
WAVEH.90
C WAVEH.91
WAVEH.92
CC local work arrays replacing WH(ip,??) for development of output WAVEH.93
CC on full grid WAVEH.94
CC work7(ip) equates to WH(ip,7) in the original WAVEH routine WAVEH.95
WAVEH.96
REAL work1(ndata),work4(ndata),work5(ndata) WAVEH.97
REAL work7(ndata),work10(ndata),work11(ndata) WAVEH.98
WAVEH.99
C WAVEH.100
INTEGER IMAX(ndata),IMSW(ndata),IMWS(ndata) WAVEH.101
C WAVEH.102
LOGICAL LFLIM,LTRUE,LWS(ndata) WAVEH.103
LOGICAL l_wvtra WAVEH.104
WAVEH.105
LOGICAL IA1(ngx,ngy) WAVEH.106
C WAVEH.107
C ---------------------------------------------------------------------- WAVEH.108
C VALUES STORED IN WH(IP,I) IN THIS SUBROUTINE WAVEH.109
C I= 1 TOTAL WAVE HEIGHT WAVEH.110
C 2 MEAN DIRECTION WAVEH.111
C 3 PRINCIPAL DIRECTION WAVEH.112
C 4 ZERO UP CROSSING PERIOD WAVEH.113
C 5 MEAN PERIOD WAVEH.114
C 6 PEAK PERIOD WAVEH.115
C 7 TOTAL WINDSEA WAVE HEIGHT WAVEH.116
C 8 MEAN DIRECTION WAVEH.117
C 9 PRINCIPAL DIRECTION WAVEH.118
C 10 ZERO UP CROSSING PERIOD WAVEH.119
C 11 MEAN PERIOD WAVEH.120
C 12 PEAK PERIOD WAVEH.121
C 13 TOTAL SWELL WAVE HEIGHT WAVEH.122
C 14 MEAN DIRECTION WAVEH.123
C 15 PRINCIPAL DIRECTION WAVEH.124
C 16 ZERO UP CROSSING PERIOD WAVEH.125
C 17 MEAN PERIOD WAVEH.126
C 18 PEAK PERIOD WAVEH.127
C 19 WAVEH.128
C 20 WAVEH.129
C ---------------------------------------------------------------------- WAVEH.130
C LOCALLY USED CONSTANTS WAVEH.131
WAVEH.132
PI2 = 2.0*PI WAVEH.133
PIO2 = PI/2.0 WAVEH.134
WAVEH.135
CCC UKMO value G14 = 0.14*G : use 0.13*G for 10m wind - WAVEH.136
CCC see wavetrain code of Anne Guillaume WAVEH.137
WAVEH.138
G14 = 0.13*G WAVEH.139
WAVEH.140
CL 13.0 ZERO DATA ARRAY WAVEH.141
C --------------------- WAVEH.142
C WAVEH.143
c * this loop initialises array WH WAVEH.144
WAVEH.145
do i=1,nfldmax WAVEH.146
DO II=1,ngrid WAVEH.147
WH(II,I) = RMDI WAVEH.148
ENDDO WAVEH.149
ENDDO WAVEH.150
WAVEH.151
do ip=1,ndata WAVEH.152
if(windsp(ip).lt.0.5) then WAVEH.153
windsp(ip)=0.5 WAVEH.154
endif WAVEH.155
enddo WAVEH.156
C WAVEH.157
EMIN = 1.0E-7 WAVEH.158
DO IP=1,NDATA WAVEH.159
ECOS(IP) = 0.0 WAVEH.160
ESIN(IP) = 0.0 WAVEH.161
ECOS1(IP) = 0.0 WAVEH.162
ESIN1(IP) = 0.0 WAVEH.163
EWS (IP) = 0.0 WAVEH.164
EMAX(IP) = EMIN WAVEH.165
EMWS(IP) = EMIN WAVEH.166
EMSW(IP) = EMIN WAVEH.167
IMAX(IP) = 0 WAVEH.168
IMWS(IP) = 0 WAVEH.169
IMSW(IP) = 0 WAVEH.170
work1(ip)=0. WAVEH.171
work4(ip)=0. WAVEH.172
work5(ip)=0. WAVEH.173
work7(ip)=0. WAVEH.174
work10(ip)=0. WAVEH.175
work11(ip)=0. WAVEH.176
ENDDO WAVEH.177
C WAVEH.178
CL 13.1 FIRST GUESS AT WINDSEA ENERGY - EWS WAVEH.179
C -------------------------------------------- WAVEH.180
C WAVEH.181
DO 10 L=1,nfre WAVEH.182
CCC DTHDF = DTH*DF(L) WAVEH.183
DTHDF = DFIM(L) WAVEH.184
DO 11 M=1,nang WAVEH.185
CC NOTE need to ensure this indexing is ok for energy WAVEH.186
NST = ((L-1)*nang+M-1)*NDATA WAVEH.187
DO 12 IP=1,NDATA WAVEH.188
FPMC = 0.8*G14/WINDSP(IP) WAVEH.189
LFLIM = FR(L).GT.FPMC WAVEH.190
IF (L.EQ.nfre) LFLIM=.TRUE. WAVEH.191
ANG = ABS(TH(M)-WINDIR(IP)) WAVEH.192
IF (ANG.GT.PI) ANG = PI2-ANG WAVEH.193
ANG = PIO2-ANG WAVEH.194
LTRUE= ANG.GT.-0.01 WAVEH.195
IF (LTRUE.AND.LFLIM) THEN ! WINDSEA WAVEH.196
EWS(IP) = EWS(IP) + ENERGY(NST+IP) * DTHDF WAVEH.197
END IF WAVEH.198
12 CONTINUE WAVEH.199
11 CONTINUE WAVEH.200
10 CONTINUE WAVEH.201
C WAVEH.202
CL 13.2 CALCULATE VALUES USING NEW CUT OFF FREQUENCY WAVEH.203
C -------------------------------------------------- WAVEH.204
DO 20 L=1,nfre WAVEH.205
ccc DTHDF = DTH*DF(L) WAVEH.206
DTHDF = DFIM(L) WAVEH.207
DO 21 M=1,nang WAVEH.208
NST = ((L-1)*nang+M-1)*NDATA WAVEH.209
DO 22 IP=1,NDATA WAVEH.210
FPM = G14/WINDSP(IP) WAVEH.211
EPM = (FPM**4)*0.0001 WAVEH.212
IF (EWS(IP).GT.EPM) EWS(IP)=EPM WAVEH.213
IF (EWS(IP).GT.0.0) THEN WAVEH.214
WAVEH.215
CC CHECK THIS COEFFICIENT 0.31 elsewhere is 0.33 ?? WAVEH.216
WAVEH.217
FPMC = 0.8*FPM*((EPM/EWS(IP))**0.31) WAVEH.218
ELSE WAVEH.219
FPMC = 0.8*G14 WAVEH.220
END IF WAVEH.221
LFLIM = FR(L).GT.FPMC WAVEH.222
IF (L.EQ.nfre) LFLIM=.TRUE. WAVEH.223
ANG = ABS(TH(M)-WINDIR(IP)) WAVEH.224
IF (ANG.GT.PI) ANG = PI2-ANG WAVEH.225
ANG = PIO2-ANG WAVEH.226
LTRUE= ANG.GT.-0.01 WAVEH.227
LWS(IP) = LTRUE.AND.LFLIM WAVEH.228
C WAVEH.229
C *** WINDSEA WAVEH.230
C WAVEH.231
IF (LWS(IP)) THEN WAVEH.232
Work7(IP) = Work7 (IP)+ENERGY(NST+IP) * DTHDF WAVEH.233
Work10(IP)= Work10(IP)+ENERGY(NST+IP)*FR(L)*FR(L)*DTHDF WAVEH.234
Work11(IP)= Work11(IP)+ENERGY(NST+IP)*FR(L)*DTHDF WAVEH.235
ECOS1(IP) = ECOS1(IP) + ENERGY(NST+IP)*COSTH(M)*DTHDF WAVEH.236
ESIN1(IP) = ESIN1(IP) + ENERGY(NST+IP)*SINTH(M)*DTHDF WAVEH.237
END IF WAVEH.238
C WAVEH.239
C *** TOTAL VALUES WAVEH.240
C WAVEH.241
cccc if(energy(nst+ip).gt.0.0001) then WAVEH.242
Work1(IP) = Work1(IP)+ENERGY(NST+IP) * DTHDF WAVEH.243
Work4(IP) = Work4(IP)+ENERGY(NST+IP)*FR(L)*FR(L)*DTHDF WAVEH.244
Work5(IP) = Work5(IP)+ENERGY(NST+IP)*FR(L)*DTHDF WAVEH.245
ECOS(IP) = ECOS(IP) + ENERGY(NST+IP)*COSTH(M)*DTHDF WAVEH.246
ESIN(IP) = ESIN(IP) + ENERGY(NST+IP)*SINTH(M)*DTHDF WAVEH.247
cccc endif WAVEH.248
WAVEH.249
22 CONTINUE WAVEH.250
C WAVEH.251
DO 24 IP=1,NDATA WAVEH.252
C WAVEH.253
C *** WINDSEA WAVEH.254
C WAVEH.255
IF (LWS(IP)) THEN WAVEH.256
IF (ENERGY(NST+IP).GT.EMWS(IP)) THEN WAVEH.257
EMWS(IP) = ENERGY(NST+IP) WAVEH.258
IMWS(IP) = L WAVEH.259
END IF WAVEH.260
ELSE WAVEH.261
C WAVEH.262
C *** SWELL WAVEH.263
C WAVEH.264
IF (ENERGY(NST+IP).GT.EMSW(IP)) THEN WAVEH.265
EMSW(IP) = ENERGY(NST+IP) WAVEH.266
IMSW(IP) = L WAVEH.267
END IF WAVEH.268
END IF WAVEH.269
C WAVEH.270
C *** TOTAL VALUES WAVEH.271
C WAVEH.272
IF (ENERGY(NST+IP).GT.EMAX(IP)) THEN WAVEH.273
EMAX(IP) = ENERGY(NST+IP) WAVEH.274
IMAX(IP) = L WAVEH.275
END IF WAVEH.276
WAVEH.277
24 CONTINUE WAVEH.278
C WAVEH.279
21 CONTINUE WAVEH.280
20 CONTINUE WAVEH.281
WAVEH.282
WAVEH.283
C WAVEH.284
C 13.3 EVALUATE SWELL VALUES FROM TOTAL AND WINDSEA WAVEH.285
C ------------------------------------------------- WAVEH.286
C WAVEH.287
WAVEH.288
cc make this a loop over ngrid ?? with test on IA1 lsmask value ?? WAVEH.289
WAVEH.290
cc DO 30 IP=1,Ngrid WAVEH.291
WAVEH.292
idata=0 WAVEH.293
WAVEH.294
do j=1,ngy WAVEH.295
WAVEH.296
do i=1,ngx WAVEH.297
WAVEH.298
ip=(j-1)*ngx + i WAVEH.299
WAVEH.300
if(.not.ia1(i,j)) then WAVEH.301
WAVEH.302
idata=idata+1 WAVEH.303
WAVEH.304
if(idata.gt.ndata)then WAVEH.305
WRITE(6,*)'error idata gt ndata ',idata,ndata GIE0F403.673
WRITE(6,*)'RETURNING FROM WAVEH ERROR CODE 1' GIE0F403.674
ICODE=1 WAVEH.308
return WAVEH.309
endif WAVEH.310
WAVEH.311
wh(ip,1) = work1(idata) WAVEH.312
wh(ip,4) = work4(idata) WAVEH.313
wh(ip,5) = work5(idata) WAVEH.314
wh(ip,7) = work7(idata) WAVEH.315
wh(ip,10)= work10(idata) WAVEH.316
wh(ip,11)= work11(idata) WAVEH.317
WAVEH.318
C WAVEH.319
C *** SWELL ENERGY WAVEH.320
C WAVEH.321
WH(IP,13) = WH(IP,1) - WH(IP,7) WAVEH.322
IF (WH(IP,13).LT.0.0) WH(IP,13) = 0.0 WAVEH.323
WAVEH.324
CC re-use work1 to hold WH(ip,13) for later use WAVEH.325
work1(idata)=wh(ip,13) WAVEH.326
WAVEH.327
C WAVEH.328
C *** SWELL INTEGRATED VALUES FOR PERIODS WAVEH.329
C WAVEH.330
WH(IP,16) = WH(IP,4) - WH(IP,10) WAVEH.331
IF (WH(IP,16).LT.0.0) WH(IP,16) = 0.0 WAVEH.332
WH(IP,17) = WH(IP,5) - WH(IP,11) WAVEH.333
IF (WH(IP,17).LT.0.0) WH(IP,17) = 0.0 WAVEH.334
C WAVEH.335
C *** TOTAL SEA DIRECTIONS - MEAN WAVEH.336
C WAVEH.337
C note for UM wave that WAM directions are from zero=North increasing WAVEH.338
C anticlockwise ( ie wind direction is atan2(u,v) which is tan -1 (u/v) WAVEH.339
C WAVEH.340
C SO here we have that ECOS is the northward component and WAVEH.341
C ESIN is the eastward component WAVEH.342
C WAVEH.343
C SO taking ATAN2(ESIN, ECOS) gives tan -1 (esin/ecos) = tan -1 (u/v) WAVEH.344
C SO the direction evaluated here is with ZERO=NORTH WAVEH.345
C and increasing clockwise WAVEH.346
C WAVEH.347
C (in radians) so only need to convert to degrees WAVEH.348
C WAVEH.349
IF(WH(IP,1).GT.0.00001) THEN WAVEH.350
WH(IP,2) = ATAN2(ESIN(Idata),ECOS(Idata)) WAVEH.351
END IF WAVEH.352
C WAVEH.353
C *** WINDSEA DIRECTIONS WAVEH.354
C WAVEH.355
c ?? 7 not 8 WAVEH.356
IF(WH(IP,7).GT.0.00001) THEN WAVEH.357
WH(IP,8) = ATAN2(ESIN1(Idata),ECOS1(Idata)) WAVEH.358
END IF WAVEH.359
C WAVEH.360
C *** NOTE CHECKS TOTAL SWELL ENERGY NONZERO WAVEH.361
C WAVEH.362
AA = WH(IP,13) WAVEH.363
IF (AA.GT.0.00001) THEN WAVEH.364
C WAVEH.365
C *** SWELL DIRECTIONS WAVEH.366
C WAVEH.367
ESIN2(Idata) = ESIN(Idata)-ESIN1(Idata) WAVEH.368
ECOS2(Idata) = ECOS(Idata)-ECOS1(Idata) WAVEH.369
WH(IP,14) = ATAN2(ESIN2(Idata),ECOS2(Idata)) WAVEH.370
END IF WAVEH.371
C WAVEH.372
C *** RESET WORK ARRAYS TO ZERO WAVEH.373
C WAVEH.374
ECOS(idata) = 0.0 WAVEH.375
ESIN(idata) = 0.0 WAVEH.376
ECOS1(idata) = 0.0 WAVEH.377
ESIN1(idata) = 0.0 WAVEH.378
ECOS2(idata) = 0.0 WAVEH.379
ESIN2(idata) = 0.0 WAVEH.380
C WAVEH.381
WAVEH.382
endif WAVEH.383
enddo WAVEH.384
enddo WAVEH.385
cc 30 CONTINUE WAVEH.386
C WAVEH.387
C 13.4 CALCULATE PRINCIPAL DIRECTIONS WAVEH.388
C ------------------------------------ WAVEH.389
C WAVEH.390
DO 45 M=1,nang WAVEH.391
DO 42 IP=1,NDATA WAVEH.392
LL=IMAX(IP) WAVEH.393
IF (LL.GT.0) THEN WAVEH.394
NST = ((LL-1)*nang+M-1)*NDATA WAVEH.395
ECOS(IP) = ECOS(IP) + ENERGY(NST+IP)*COSTH(M) WAVEH.396
ESIN(IP) = ESIN(IP) + ENERGY(NST+IP)*SINTH(M) WAVEH.397
END IF WAVEH.398
L1=IMWS(IP) WAVEH.399
IF (L1.GT.0) THEN WAVEH.400
NS1 = ((L1-1)*nang+M-1)*NDATA WAVEH.401
ECOS1(IP) = ECOS1(IP) + ENERGY(NS1+IP)*COSTH(M) WAVEH.402
ESIN1(IP) = ESIN1(IP) + ENERGY(NS1+IP)*SINTH(M) WAVEH.403
END IF WAVEH.404
CC need to sort out this use of WH(,13) WAVEH.405
CC WH(IP,13) is copied into data work array work1 in previous loop WAVEH.406
CC IF (WH(IP,13).GT.0.00001) THEN WAVEH.407
IF (Work1(IP).GT.0.00001) THEN WAVEH.408
L2=IMSW(IP) WAVEH.409
NS2 = ((L2-1)*nang+M-1)*NDATA WAVEH.410
ECOS2(IP) = ECOS2(IP) + ENERGY(NS2+IP)*COSTH(M) WAVEH.411
ESIN2(IP) = ESIN2(IP) + ENERGY(NS2+IP)*SINTH(M) WAVEH.412
END IF WAVEH.413
42 CONTINUE WAVEH.414
45 CONTINUE WAVEH.415
C WAVEH.416
WAVEH.417
CC MAKE THIS LOOP over ngx ngy WAVEH.418
WAVEH.419
idata=0 ! must reset idata before next use WAVEH.420
WAVEH.421
cc DO IP=1,NDATA WAVEH.422
do j=1,ngy WAVEH.423
do i=1,ngx WAVEH.424
WAVEH.425
ip=(j-1)*ngx + i ! counter for points on full grid WAVEH.426
WAVEH.427
if(.not.ia1(i,j)) then WAVEH.428
WAVEH.429
idata=idata+1 ! counter for data points only WAVEH.430
WAVEH.431
if(idata.gt.ndata)then WAVEH.432
WRITE(6,*)'error idata gt ndata ',idata,ndata GIE0F403.675
WRITE(6,*)'setting error code and returning' GIE0F403.676
icode=99 WAVEH.435
return WAVEH.436
endif WAVEH.437
WAVEH.438
IF(IMAX(idata).GT.0)WH(IP,3) = ATAN2(ESIN(idata),ECOS(idata)) WAVEH.439
IF(IMWS(idata).GT.0)WH(IP,9) = ATAN2(ESIN1(idata),ECOS1(idata)) WAVEH.440
IF (WH(IP,13).GT.0.00001.AND.IMSW(idata).GT.0) THEN WAVEH.441
WH(IP,15) = ATAN2(ESIN2(idata),ECOS2(idata)) WAVEH.442
ENDIF WAVEH.443
C WAVEH.444
C WAVEH.445
C *** TOTAL SPECTRUM - ZERO UP CROSSING PERIOD WAVEH.446
C - & MEAN PERIOD WAVEH.447
C - & PEAK PERIOD WAVEH.448
C WAVEH.449
IF(WH(IP,4).GT.0.00001) THEN WAVEH.450
WH(IP,4) = SQRT(WH(IP,1)/WH(IP,4)) WAVEH.451
ENDIF WAVEH.452
WAVEH.453
IF(WH(IP,5).GT.0.00001) THEN WAVEH.454
WH(IP,5) = WH(IP,1)/WH(IP,5) WAVEH.455
ENDIF WAVEH.456
WAVEH.457
IF(IMAX(idata).GT.0) WH(IP,6) = 1./FR(IMAX(idata)) WAVEH.458
C WAVEH.459
C *** WIND SEA PERIODS WAVEH.460
C WAVEH.461
IF (WH(IP,10).GT.0.0) THEN WAVEH.462
WH(IP,10) = SQRT(WH(IP,7)/WH(IP,10)) WAVEH.463
ELSE WAVEH.464
WH(IP,10) = 0.0 WAVEH.465
END IF WAVEH.466
IF (WH(IP,11).GT.0.0) THEN WAVEH.467
WH(IP,11) = WH(IP,7)/WH(IP,11) WAVEH.468
ELSE WAVEH.469
WH(IP,11) = 0.0 WAVEH.470
ENDIF WAVEH.471
IF(IMWS(idata).GT.0) WH(IP,12) = 1./FR(IMWS(idata)) WAVEH.472
C WAVEH.473
C *** NOTE CHECKS TOTAL SWELL ENERGY NONZERO WAVEH.474
C *** SWELL PERIODS WAVEH.475
C WAVEH.476
IF (WH(IP,16).GT.0.000001) THEN WAVEH.477
WH(IP,16) = SQRT(WH(IP,13)/WH(IP,16)) WAVEH.478
ELSE WAVEH.479
WH(IP,16) = 0.0 WAVEH.480
END IF WAVEH.481
IF (WH(IP,17).GT.0.000001) THEN WAVEH.482
WH(IP,17) = WH(IP,13)/WH(IP,17) WAVEH.483
ELSE WAVEH.484
WH(IP,17) = 0.0 WAVEH.485
END IF WAVEH.486
AA = WH(IP,13) WAVEH.487
IF (AA.GT.0.00001) THEN WAVEH.488
IF(IMSW(idata).GT.0) WH(IP,18) = 1./FR(IMSW(idata)) WAVEH.489
ELSE WAVEH.490
WH(IP,18) = 0.0 WAVEH.491
END IF WAVEH.492
C WAVEH.493
C *** TOTAL WAVE HEIGHT WAVEH.494
C WAVEH.495
WH(IP,1) = 4.*SQRT(WH(IP,1)) WAVEH.496
C WAVEH.497
C *** WINDSEA WAVE HEIGHT WAVEH.498
C WAVEH.499
WH(IP,7) = 4.*SQRT(WH(IP,7)) WAVEH.500
C WAVEH.501
C *** SWELL WAVE HEIGHT WAVEH.502
C WAVEH.503
WH(IP,13) = 4.*SQRT(WH(IP,13)) WAVEH.504
C WAVEH.505
WAVEH.506
c convert directions to degrees in range 0-360 WAVEH.507
WAVEH.508
WH(ip,2)=wh(ip,2)*RECIP_PI_OVER_180 WAVEH.509
WH(ip,3)=wh(ip,3)*RECIP_PI_OVER_180 WAVEH.510
WH(ip,8)=wh(ip,8)*RECIP_PI_OVER_180 WAVEH.511
WH(ip,9)=wh(ip,9)*RECIP_PI_OVER_180 WAVEH.512
WH(ip,14)=wh(ip,14)*RECIP_PI_OVER_180 WAVEH.513
WH(ip,15)=wh(ip,15)*RECIP_PI_OVER_180 WAVEH.514
WAVEH.515
if(wh(ip,2).lt.0.) wh(ip,2)=wh(ip,2)+360. WAVEH.516
if(wh(ip,3).lt.0.) wh(ip,3)=wh(ip,3)+360. WAVEH.517
if(wh(ip,8).lt.0.) wh(ip,8)=wh(ip,8)+360. WAVEH.518
if(wh(ip,9).lt.0.) wh(ip,9)=wh(ip,9)+360. WAVEH.519
if(wh(ip,14).lt.0.) wh(ip,14)=wh(ip,14)+360. WAVEH.520
if(wh(ip,15).lt.0.) wh(ip,15)=wh(ip,15)+360. WAVEH.521
WAVEH.522
if(wh(ip,2).gt.360.) wh(ip,2)=wh(ip,2)-360. WAVEH.523
if(wh(ip,3).gt.360.) wh(ip,3)=wh(ip,3)-360. WAVEH.524
if(wh(ip,8).gt.360.) wh(ip,8)=wh(ip,8)-360. WAVEH.525
if(wh(ip,9).gt.360.) wh(ip,9)=wh(ip,9)-360. WAVEH.526
if(wh(ip,14).gt.360.) wh(ip,14)=wh(ip,14)-360. WAVEH.527
if(wh(ip,15).gt.360.) wh(ip,15)=wh(ip,15)-360. WAVEH.528
WAVEH.529
WAVEH.530
endif ! check for sea points if(ia1.... WAVEH.531
enddo WAVEH.532
enddo WAVEH.533
WAVEH.534
C WAVEH.535
if(l_wvtra) then WAVEH.536
WAVEH.537
cccccc len1=nfre*nang*ndata WAVEH.538
WAVEH.539
ccc call wavetr(energy,pswh,perio,pdir,kwtot,fr,dfim,th, WAVEH.540
ccc +len1,ndata,kwtmax,nfre,nang,rmdi,icode) WAVEH.541
WAVEH.542
CC pointers to wavetrain data in array WH WAVEH.543
CC usage altered from UKMO model, so that now array WH has WAVEH.544
CC all wvtrain heights together then all periods together then WAVEH.545
CC all directions together WAVEH.546
CC (means we don't need array pswh / perio / pdir in this routine) WAVEH.547
WAVEH.548
i_SWH=19 WAVEH.549
i_PER=i_SWH + kwtmax WAVEH.550
i_DIR=i_PER + kwtmax WAVEH.551
i_NUM=i_DIR + kwtmax WAVEH.552
WAVEH.553
call wavetr(energy,WH(1,i_swh),WH(1,i_per),WH(1,i_dir), WAVEH.554
+kwtot,fr,dfim,th, WAVEH.555
+len1,ndata,kwtmax,nfre,nang,rmdi,icode) WAVEH.556
WAVEH.557
C * here add test return code if ne 0 then error trap * WAVEH.558
C WAVEH.559
C WAVEH.560
do ip=1,ndata WAVEH.561
C wave heights WAVEH.562
ccc wh(ip,25)=pswh(ip,1) WAVEH.563
ccc wh(ip,28)=pswh(ip,2) WAVEH.564
ccc wh(ip,31)=pswh(ip,3) WAVEH.565
WAVEH.566
C wave periods WAVEH.567
ccc wh(ip,26)=perio(ip,1) WAVEH.568
ccc wh(ip,29)=perio(ip,2) WAVEH.569
ccc wh(ip,32)=perio(ip,3) WAVEH.570
WAVEH.571
C wave directions WAVEH.572
ccc wh(ip,27)=pdir(ip,1) WAVEH.573
ccc wh(ip,30)=pdir(ip,2) WAVEH.574
ccc wh(ip,33)=pdir(ip,3) WAVEH.575
WAVEH.576
C number of trains WAVEH.577
wh(ip,i_num) = 1.0*kwtot(ip) WAVEH.578
C WAVEH.579
enddo WAVEH.580
c WAVEH.581
endif WAVEH.582
WAVEH.583
RETURN WAVEH.584
END WAVEH.585
*ENDIF WAVEH.586