WAVEWATCH III  beta 0.0.1
Functions/Subroutines | Variables
w3snl3md Module Reference

Generalized and optimized multiple DIA implementation. More...

Functions/Subroutines

subroutine w3snl3 (A, CG, WN, DEPTH, S, D)
 Multiple Discrete Interaction Parameterization for arbitrary depths with generalized quadruplet layout. More...
 
subroutine insnl3
 Initialization for generalized multiple DIA routine. More...
 

Variables

real, parameter, public lammax = 0.49999
 
real, parameter, public delthm = 90.
 

Detailed Description

Generalized and optimized multiple DIA implementation.

Expressions in terms of original F(f,theta) spectrum.

Author
H. L. Tolman
Date
13-Jul-2012
Copyright
Copyright 2009-2022 National Weather Service (NWS), National Oceanic and Atmospheric Administration. All rights reserved. WAVEWATCH III is a trademark of the NWS. No unauthorized use without permission.

Function/Subroutine Documentation

◆ insnl3()

subroutine w3snl3md::insnl3

Initialization for generalized multiple DIA routine.

Fill storage arrays as described in the main subroutine with interpolation, model and distribution data.

Author
H. L. Tolman
Date
13-Nov-2009

Definition at line 788 of file w3snl3md.F90.

788  !/
789  !/ +-----------------------------------+
790  !/ | WAVEWATCH-III NOAA/NCEP |
791  !/ | H. L. Tolman |
792  !/ | FORTRAN 90 |
793  !/ | Last update : 13-Nov-2009 |
794  !/ +-----------------------------------+
795  !/
796  !/ 21-Jul-2008 : Origination as NLX option. ( version 3.13 )
797  !/ 03-Jan-2009 : Bug fixes NTHMAX and NTHMX2. ( version 3.13 )
798  !/ 21-Aug-2009 : Conversion to F(f,theta) form. ( version 3.13 )
799  !/ 13-Nov-2009 : Harden DELTH computation. ( version 3.13 )
800  !/
801  ! 1. Purpose :
802  !
803  ! Initialization for generalized multiple DIA routine.
804  !
805  ! 2. Method :
806  !
807  ! Fill storage aryays as described in the main subroutine with
808  ! interpolation, model and distribution data.
809  !
810  ! 3. Parameters :
811  !
812  ! Variables in W3GDATMD describing model setup.
813  !
814  ! Name Type Scope Description
815  ! ----------------------------------------------------------------
816  ! SNLNQ Int. Public Number of quadruplet definitions.
817  ! SNLL R.A. Public Array with lambda for quadruplet.
818  ! SNLM R.A. Public Array with mu for quadruplet.
819  ! SNLT R.A. Public Array with Dtheta for quadruplet.
820  ! SNLCD R.A. Public Array with Cd for quadruplet.
821  ! SNLCS R.A. Public Array with Cs for quadruplet.
822  ! ----------------------------------------------------------------
823  !
824  ! 4. Subroutines used :
825  !
826  ! Name Type Module Description
827  ! ----------------------------------------------------------------
828  ! STRACE Subr. W3SERVMD Subroutine tracing.
829  ! EXTCDE Subr. W3SERVMD Program abort.
830  ! WAVNU2 Subr. W3DISPMD Solve dispersion relation.
831  ! ----------------------------------------------------------------
832  !
833  ! 5. Called by :
834  !
835  ! Name Type Module Description
836  ! ----------------------------------------------------------------
837  ! W3IOGR Subr. W3IOGRMD Process model definiton file.
838  ! ----------------------------------------------------------------
839  !
840  ! 6. Error messages :
841  !
842  ! See error escape location.
843  !
844  ! 8. Remarks :
845  !
846  ! - Allocation of arrays directly done in data structure, using
847  ! IGRID and resetting pointer of aliaases.
848  ! - In the 03-Jan-2009 bug fix !/T3 error output was fixed, and
849  ! NTHMAX is increased by 1 to assure that NTHMX2 .LE. NTHMAX
850  ! for any lambda and mu. With this, the label 810 test is
851  ! changed from equality testing to .LE. testing.
852  !
853  ! 8. Structure :
854  !
855  ! See source code.
856  !
857  ! 9. Switches :
858  !
859  ! !/S Enable subroutine tracing.
860  ! !/T General test output.
861  ! !/T1 Filling of lookup table for quadruplet and interaction
862  ! strength.
863  ! !/T2 Filling of lookup table for combining interactions.
864  ! !/T3 Display raw lookup table of second type.
865  !
866  ! 10. Source code :
867  !
868  !/ ------------------------------------------------------------------- /
869  USE constants
870  USE w3odatmd, ONLY: ndse, ndst
871  USE w3gdatmd, nfr => nk
872  !
873  USE w3dispmd, ONLY: wavnu2
874  USE w3servmd, ONLY: extcde
875 #ifdef W3_S
876  USE w3servmd, ONLY: strace
877 #endif
878  !/
879  IMPLICIT NONE
880  !/
881  !/ ------------------------------------------------------------------- /
882  !/ Parameter list
883  !/
884  !/ ------------------------------------------------------------------- /
885  !/ Local parameters
886  !/
887  INTEGER :: IFRMIN, IFRMAX, IKD, IERR, IQ, NQD, &
888  NQS, J, IFR, IQA, JJ, JF, NTHMX2, &
889  JIQ, JOF, JQR, IST
890  INTEGER :: JFR(4), JFR1(4), JTH(4), JTH1(4)
891 #ifdef W3_S
892  INTEGER, SAVE :: IENT = 0
893 #endif
894  INTEGER, ALLOCATABLE :: AST1(:,:,:), AST2(:,:,:)
895  REAL :: SITMAX, XFRLN
896  REAL :: OFF12, OFF34, TH12, DEPTH, &
897  S0, S1, S2, S3, S4, AUXFR(4), &
898  WN0, WN1, WN2, WN3, WN4, &
899  CG0, CG1, CG2, CG3, CG4, AUXF, &
900  AA, BB, CC, DELTH(4), AUX1, AUX2, &
901  WFR(4), WFR1(4), WTH(4), WTH1(4), &
902  WFROFF, SIOFF, WF
903  !
904  TYPE qst
905  INTEGER :: OFR(4), OFR1(4), OTH(4), OTH1(4)
906  REAL :: HFR(4), HFR1(4), HTH(4), HTH1(4)
907  REAL :: F1, F2, F3, F4, CQD, CQS
908  END TYPE qst
909  !
910  TYPE(QST), ALLOCATABLE :: TSTORE(:,:)
911  !/
912  !/ ------------------------------------------------------------------- /
913  !/
914 #ifdef W3_S
915  CALL strace (ient, 'INSNL3')
916 #endif
917  !
918  ! 1. Initialization ------------------------------------------------- *
919  ! 1.a Checks
920  !
921  xfrln = log(xfr)
922  !
923  IF ( lammax.LE.0. .OR. lammax.GT.0.5 .OR. delthm.LT.0. ) GOTO 800
924  !
925  ! 1.b Set up relative depths
926  !
927  ALLOCATE ( tstore(snlnq*4,1:nkd) )
928  !
929  depth = 1.
930  sitmin = sqrt( kdmin * tanh(kdmin) )
931  sitmax = sqrt( kdmax * tanh(kdmax) )
932  xsit = (sitmax/sitmin)**(1./real(nkd-1))
933  !
934 #ifdef W3_T
935  WRITE (ndst,9010) nkd, kdmin, kdmax, xsit
936 #endif
937  !
938  ! 2. Building quadruplet data base ---------------------------------- *
939  ! For quadruplet and interaction strength evaluation
940  !
941  ifrmin = 0
942  ifrmax = 0
943  nthmax = 0
944  !
945  ! 2.a Loop over relative depths
946  !
947  s0 = sitmin * sqrt( grav / depth ) / xsit
948  !
949  DO ikd=1, nkd
950  !
951  s0 = s0 * xsit
952  CALL wavnu2 ( s0, depth, wn0, cg0, 1.e-6, 25, ierr)
953  !
954  ! 2.b Loop over representative quadruplets
955  !
956  nqa = 0
957  nqd = 0
958  nqs = 0
959  !
960  DO iq=1, snlnq
961  !
962 #ifdef W3_T1
963  WRITE (ndst,9020) ikd, iq, wn0*depth, s0*tpiinv, depth
964 #endif
965  !
966  off12 = snlm(iq)
967  off34 = snll(iq)
968  th12 = snlt(iq) * dera
969  IF ( snlcd(iq) .GT. 0. ) nqd = nqd + 1
970  IF ( snlcs(iq) .GT. 0. ) nqs = nqs + 1
971  !
972  IF ( th12 .LT. 0. ) THEN
973  IF ( off12.LT.0. .OR. off12.GT.0.5 .OR. &
974  off34.LT.0. .OR. off34.GT.0.5 ) GOTO 801
975  ELSE
976  IF ( snlt(iq).GT.delthm .OR. off12.LT.0. .OR. &
977  off12.GE.1. &
978  .OR. off34.LT.minlam(off12,snlt(iq)) .OR. &
979  off34.GT.maxlam(off12,snlt(iq)) ) GOTO 802
980  END IF
981  !
982 #ifdef W3_T1
983  WRITE (ndst,9021) snlt(iq), off12, off34, &
984  snlcd(iq), snlcs(iq)
985 #endif
986  !
987  ! 2.c Offset angles
988  !
989  s1 = s0 * ( 1. + off12 )
990  CALL wavnu2 ( s1, depth, wn1, cg1, 1.e-6, 25, ierr)
991  s2 = s0 * ( 1. - off12 )
992  CALL wavnu2 ( s2, depth, wn2, cg2, 1.e-6, 25, ierr)
993  s3 = s0 * ( 1. + off34 )
994  CALL wavnu2 ( s3, depth, wn3, cg3, 1.e-6, 25, ierr)
995  s4 = s0 * ( 1. - off34 )
996  CALL wavnu2 ( s4, depth, wn4, cg4, 1.e-6, 25, ierr)
997  !
998  auxfr(1) = s1 / s0
999  auxfr(2) = s2 / s0
1000  auxfr(3) = s3 / s0
1001  auxfr(4) = s4 / s0
1002  !
1003  IF ( th12 .LT. 0. ) THEN
1004  bb = 2. * wn0
1005  ELSE
1006  bb = wn1**2 + wn2**2 + 2.*wn1*wn2*cos(th12)
1007  bb = sqrt( max( bb , 0. ) )
1008  END IF
1009  !
1010  IF ( th12.LT.0. .AND. abs(off12).LE.1.e-4 ) THEN
1011  delth(1) = 0.
1012  delth(2) = 0.
1013  ELSE
1014  cc = wn1
1015  aa = wn2
1016  aux1 = (cc**2+bb**2-aa**2) / (2.*bb*cc)
1017  aux2 = (aa**2+bb**2-cc**2) / (2.*bb*aa)
1018  delth(1) = - acos( max( 0. , min( 1. , aux1 ) ) )
1019  delth(2) = acos( max( 0. , min( 1. , aux2 ) ) )
1020  END IF
1021  cc = wn3
1022  aa = wn4
1023  aux1 = (cc**2+bb**2-aa**2) / (2.*bb*cc)
1024  aux2 = (aa**2+bb**2-cc**2) / (2.*bb*aa)
1025  delth(3) = - acos( max( 0. , min( 1. , aux1 ) ) )
1026  delth(4) = acos( max( 0. , min( 1. , aux2 ) ) )
1027  !
1028 #ifdef W3_T1
1029  WRITE (ndst,9022) delth(:) * rade
1030 #endif
1031  !
1032  ! 2.d Frequency indices
1033  !
1034  DO j=1, 4
1035  jfr(j) = int( log(auxfr(j)) / xfrln )
1036  jfr1(j) = jfr(j) + 1 * sign(1.,auxfr(j)-1.)
1037  wfr(j) = (xfr**jfr1(j)-auxfr(j))/(xfr**jfr1(j)-xfr**jfr(j))
1038  wfr1(j) = 1. - wfr(j)
1039  END DO
1040  !
1041  ifrmin = min( ifrmin , minval(jfr1) )
1042  ifrmax = max( ifrmax , maxval(jfr1) )
1043  !
1044 #ifdef W3_T1
1045  WRITE (ndst,9023) 1, jfr(1), jfr1(1), wfr(1), wfr1(1)
1046  DO, j=2, 4
1047  WRITE (ndst,9024) j, jfr(j), jfr1(j), wfr(j), wfr1(j)
1048  END DO
1049 #endif
1050  !
1051  ! 2.e Directional indices
1052  !
1053  DO j=1, 4
1054  aux1 = delth(j) / dth
1055  jth(j) = int(aux1)
1056  jth1(j) = jth(j) + 1 * sign(1.,delth(j))
1057  wth1(j) = abs(aux1) - real(abs(jth(j)))
1058  wth(j) = 1. - wth1(j)
1059  END DO
1060  !
1061  nthmax = max( nthmax , maxval(abs(jth1)) )
1062  !
1063 #ifdef W3_T1
1064  WRITE (ndst,9025) 1, jth(1), jth1(1), wth(1), wth1(1)
1065  DO, j=2, 4
1066  WRITE (ndst,9024) j, jth(j), jth1(j), wth(j), wth1(j)
1067  END DO
1068 #endif
1069  !
1070  ! 2.f Temp storage of data
1071  !
1072  IF ( snlm(iq).EQ.0. .AND. snlt(iq).LT.0. ) THEN
1073  jj = 2
1074  ELSE
1075  jj = 4
1076  END IF
1077  !
1078  DO j=1, jj
1079  SELECT CASE (j)
1080  CASE (2)
1081  jth(3) = -jth(3)
1082  jth(4) = -jth(4)
1083  jth1(3) = -jth1(3)
1084  jth1(4) = -jth1(4)
1085  CASE (3)
1086  jth = -jth
1087  jth1 = -jth1
1088  CASE (4)
1089  jth(3) = -jth(3)
1090  jth(4) = -jth(4)
1091  jth1(3) = -jth1(3)
1092  jth1(4) = -jth1(4)
1093  CASE DEFAULT
1094  END SELECT
1095  !
1096  nqa = nqa + 1
1097  tstore(nqa,ikd)%OFR = jfr
1098  tstore(nqa,ikd)%OFR1 = jfr1
1099  tstore(nqa,ikd)%HFR = wfr
1100  tstore(nqa,ikd)%HFR1 = wfr1
1101  tstore(nqa,ikd)%OTH = jth
1102  tstore(nqa,ikd)%OTH1 = jth1
1103  tstore(nqa,ikd)%HTH = wth
1104  tstore(nqa,ikd)%HTH1 = wth1
1105  IF ( jj .EQ. 2 ) THEN
1106  tstore(nqa,ikd)%CQD = snlcd(iq) * 2.
1107  tstore(nqa,ikd)%CQS = snlcs(iq) * 2.
1108  ELSE
1109  tstore(nqa,ikd)%CQD = snlcd(iq)
1110  tstore(nqa,ikd)%CQS = snlcs(iq)
1111  END IF
1112  auxf = ( wn0 * s0 ) / cg0
1113  tstore(nqa,ikd)%F1 = auxf * cg1 / ( wn1 * s1 )
1114  tstore(nqa,ikd)%F2 = auxf * cg2 / ( wn2 * s2 )
1115  tstore(nqa,ikd)%F3 = auxf * cg3 / ( wn3 * s3 )
1116  tstore(nqa,ikd)%F4 = auxf * cg4 / ( wn4 * s4 )
1117  !
1118  END DO
1119  !
1120  ! ... End loop 2.b
1121  !
1122  END DO
1123  !
1124  ! ... End loop 2.a
1125  !
1126  END DO
1127  !
1128 #ifdef W3_T1
1129  WRITE (ndst,*)
1130 #endif
1131 #ifdef W3_T
1132  WRITE (ndst,9026) nqa, snlnq*4, nqd, nqs
1133 #endif
1134  !
1135  ! 2.g Expanded spectral range
1136  !
1137  nthmax = nthmax + 1
1138  !
1139  nfrmin = 1 + ifrmin
1140  nfrmax = nfr + ifrmax - ifrmin
1141  nfrcut = nfr - ifrmin
1142  nthexp = nth + 2*nthmax
1143  !
1144  nspmin = 1 + (nfrmin-1)*nthexp - nthmax
1145  nspmax = nfrmax * nthexp - nthmax
1146  nspmx2 = nfrcut * nthexp - nthmax
1147  !
1148 #ifdef W3_T
1149  WRITE (ndst,9027) nfr, nfrmin, nfrmax, nfrcut, nth, &
1150  1-nthmax, nth+nthmax, nthexp
1151 #endif
1152  !
1153  ALLOCATE ( mpars(igrid)%SNLPS%FRQ(nfrmax), &
1154  mpars(igrid)%SNLPS%XSI(nfrmax) )
1155  frq => mpars(igrid)%SNLPS%FRQ
1156  xsi => mpars(igrid)%SNLPS%XSI
1157  !
1158  xsi(1:nfr) = sig(1:nfr)
1159  DO ifr=nfr+1, nfrmax
1160  xsi(ifr) = xsi(ifr-1) * xfr
1161  END DO
1162  frq = xsi * tpiinv
1163  !
1164  ! 2.h Final storage
1165  !
1166  ALLOCATE ( mpars(igrid)%SNLPS%QST1(16,nqa,nkd), &
1167  mpars(igrid)%SNLPS%QST3(6,nqa,nkd), &
1168  mpars(igrid)%SNLPS%QST2(16,nqa,nkd) )
1169  qst1 => mpars(igrid)%SNLPS%QST1
1170  qst2 => mpars(igrid)%SNLPS%QST2
1171  qst3 => mpars(igrid)%SNLPS%QST3
1172  !
1173  ! 2.h.1 Basic data
1174  !
1175  DO ikd=1, nkd
1176  DO iqa=1, nqa
1177  !
1178  DO j=1, 4
1179  !
1180  qst1((j-1)*4+1,iqa,ikd) = tstore(iqa,ikd)%OTH (j) + &
1181  tstore(iqa,ikd)%OFR (j) * nthexp
1182  qst1((j-1)*4+2,iqa,ikd) = tstore(iqa,ikd)%OTH1(j) + &
1183  tstore(iqa,ikd)%OFR (j) * nthexp
1184  qst1((j-1)*4+3,iqa,ikd) = tstore(iqa,ikd)%OTH (j) + &
1185  tstore(iqa,ikd)%OFR1(j) * nthexp
1186  qst1((j-1)*4+4,iqa,ikd) = tstore(iqa,ikd)%OTH1(j) + &
1187  tstore(iqa,ikd)%OFR1(j) * nthexp
1188  !
1189  qst2((j-1)*4+1,iqa,ikd) = tstore(iqa,ikd)%HFR (j) * &
1190  tstore(iqa,ikd)%HTH (j)
1191  qst2((j-1)*4+2,iqa,ikd) = tstore(iqa,ikd)%HFR (j) * &
1192  tstore(iqa,ikd)%HTH1(j)
1193  qst2((j-1)*4+3,iqa,ikd) = tstore(iqa,ikd)%HFR1(j) * &
1194  tstore(iqa,ikd)%HTH (j)
1195  qst2((j-1)*4+4,iqa,ikd) = tstore(iqa,ikd)%HFR1(j) * &
1196  tstore(iqa,ikd)%HTH1(j)
1197  !
1198  END DO
1199  !
1200  qst3(1,iqa,ikd) = tstore(iqa,ikd)%F1
1201  qst3(2,iqa,ikd) = tstore(iqa,ikd)%F2
1202  qst3(3,iqa,ikd) = tstore(iqa,ikd)%F3
1203  qst3(4,iqa,ikd) = tstore(iqa,ikd)%F4
1204  qst3(5,iqa,ikd) = tstore(iqa,ikd)%CQD
1205  qst3(6,iqa,ikd) = tstore(iqa,ikd)%CQS
1206  !
1207  END DO
1208  END DO
1209  !
1210  IF ( nqd .GT. 0 ) qst3(5,:,:) = qst3(5,:,:) / real(nqd)
1211  IF ( nqs .GT. 0 ) qst3(6,:,:) = qst3(6,:,:) / real(nqs)
1212  !
1213  DEALLOCATE ( tstore )
1214  !
1215  ! 3. Building quadruplet data base ---------------------------------- *
1216  ! For constructing interactions and diagonal from contributions
1217  !
1218  nthmx2 = 0
1219  ALLOCATE ( mpars(igrid)%SNLPS%QST4(16,nqa,nkd), &
1220  mpars(igrid)%SNLPS%QST5(16,nqa,nkd), &
1221  mpars(igrid)%SNLPS%QST6(16,nqa,nkd) )
1222  qst4 => mpars(igrid)%SNLPS%QST4
1223  qst5 => mpars(igrid)%SNLPS%QST5
1224  qst6 => mpars(igrid)%SNLPS%QST6
1225  ALLOCATE ( ast1(16,nqa,nkd), ast2(16,nqa,nkd) )
1226  !
1227  ! 3.a Loop over relative depths
1228  !
1229  s0 = sitmin * sqrt( grav / depth ) / xsit
1230  !
1231  DO ikd=1, nkd
1232  !
1233  s0 = s0 * xsit
1234  CALL wavnu2 ( s0, depth, wn0, cg0, 1.e-6, 25, ierr)
1235  !
1236  ! 3.b Loop over representative quadruplets
1237  !
1238  nqa = 0
1239  !
1240  DO iq=1, snlnq
1241  !
1242 #ifdef W3_T2
1243  WRITE (ndst,9030) ikd, iq, wn0*depth, s0*tpiinv, depth
1244 #endif
1245  !
1246  off12 = snlm(iq)
1247  off34 = snll(iq)
1248  th12 = snlt(iq) * dera
1249  !
1250 #ifdef W3_T2
1251  WRITE (ndst,9031) snlt(iq), off12, off34
1252 #endif
1253  !
1254  ! 3.c Frequency indices
1255  !
1256  auxfr(1) = ( 1. + off12 )
1257  auxfr(2) = ( 1. - off12 )
1258  auxfr(3) = ( 1. + off34 )
1259  auxfr(4) = ( 1. - off34 )
1260  !
1261  DO j=1, 4
1262  jfr(j) = int( log(auxfr(j)) / xfrln )
1263  jfr1(j) = jfr(j) + 1 * sign(1.,auxfr(j)-1.)
1264  wfr(j) = (xfr**jfr1(j)-auxfr(j))/(xfr**jfr1(j)-xfr**jfr(j))
1265  wfr1(j) = 1. - wfr(j)
1266  END DO
1267  !
1268 #ifdef W3_T2
1269  WRITE (ndst,9032) 1, jfr(1), jfr1(1), wfr(1), wfr1(1)
1270  DO, j=2, 4
1271  WRITE (ndst,9033) j, jfr(j), jfr1(j), wfr(j), wfr1(j)
1272  END DO
1273 #endif
1274  !
1275  ! 3.d Loop over quadruplet components
1276  !
1277  DO jiq=1, 4
1278  !
1279  IF ( jiq .LE. 2 ) THEN
1280  wf = -1.
1281  ELSE
1282  wf = 1.
1283  END IF
1284  !
1285  ! 3.e Loop over frequency offsets, get directional offsets
1286  !
1287  DO jof=1, 2
1288  !
1289  IF ( jof .EQ. 1 ) THEN
1290  ifr = -jfr(jiq)
1291  wfroff = wfr(jiq)
1292  ELSE
1293  ifr = -jfr1(jiq)
1294  wfroff = wfr1(jiq)
1295  END IF
1296  !
1297  sioff = s0 * xfr**ifr
1298  CALL wavnu2 ( sioff, depth, wn0, cg0, 1.e-6, 25, ierr)
1299  s1 = sioff * ( 1. + off12 )
1300  CALL wavnu2 ( s1, depth, wn1, cg1, 1.e-6, 25, ierr)
1301  s2 = sioff * ( 1. - off12 )
1302  CALL wavnu2 ( s2, depth, wn2, cg2, 1.e-6, 25, ierr)
1303  s3 = sioff * ( 1. + off34 )
1304  CALL wavnu2 ( s3, depth, wn3, cg3, 1.e-6, 25, ierr)
1305  s4 = sioff * ( 1. - off34 )
1306  CALL wavnu2 ( s4, depth, wn4, cg4, 1.e-6, 25, ierr)
1307  !
1308 #ifdef W3_T2
1309  WRITE (ndst,9034) jiq, jof, ifr, wfroff, sioff/s0
1310 #endif
1311  !
1312  IF ( th12 .LT. 0. ) THEN
1313  bb = 2. * wn0
1314  ELSE
1315  bb = wn1**2 + wn2**2 + 2.*wn1*wn2*cos(th12)
1316  bb = sqrt( max( bb , 0. ) )
1317  END IF
1318  !
1319  IF ( th12.LT.0. .AND. abs(off12).LE.1.e-4 ) THEN
1320  delth(1) = 0.
1321  delth(2) = 0.
1322  ELSE
1323  cc = wn1
1324  aa = wn2
1325  aux1 = (cc**2+bb**2-aa**2) / (2.*bb*cc)
1326  aux2 = (aa**2+bb**2-cc**2) / (2.*bb*aa)
1327  delth(1) = - acos( max( 0. , min( 1. , aux1 ) ) )
1328  delth(2) = acos( max( 0. , min( 1. , aux2 ) ) )
1329  END IF
1330  cc = wn3
1331  aa = wn4
1332  aux1 = (cc**2+bb**2-aa**2) / (2.*bb*cc)
1333  aux2 = (aa**2+bb**2-cc**2) / (2.*bb*aa)
1334  delth(3) = - acos( max( 0. , min( 1. , aux1 ) ) )
1335  delth(4) = acos( max( 0. , min( 1. , aux2 ) ) )
1336  !
1337 #ifdef W3_T2
1338  WRITE (ndst,9035) delth(:) * rade
1339 #endif
1340  !
1341  aux1 = delth(jiq) / dth
1342  jth(jiq) = int(aux1)
1343  jth1(jiq) = jth(jiq) + 1 * sign(1.,delth(jiq))
1344  wth1(jiq) = abs(aux1) - real(abs(jth(jiq)))
1345  wth(jiq) = 1. - wth1(jiq)
1346  !
1347  nthmx2 = max( nthmx2 , abs(jth1(jiq)) )
1348  !
1349 #ifdef W3_T2
1350  WRITE (ndst,9036) jiq, jth(jiq), jth1(jiq), &
1351  wth(jiq), wth1(jiq)
1352 #endif
1353  !
1354  ! 3.f Loop over quadruplet realizations
1355  !
1356  IF ( snlm(iq).EQ.0. .AND. snlt(iq).LT.0. ) THEN
1357  jj = 2
1358  ELSE
1359  jj = 4
1360  END IF
1361  !
1362  DO jqr=1, jj
1363  !
1364  SELECT CASE (jqr)
1365  CASE (2)
1366  jth(3) = -jth(3)
1367  jth(4) = -jth(4)
1368  jth1(3) = -jth1(3)
1369  jth1(4) = -jth1(4)
1370  CASE (3)
1371  jth = -jth
1372  jth1 = -jth1
1373  CASE (4)
1374  jth(3) = -jth(3)
1375  jth(4) = -jth(4)
1376  jth1(3) = -jth1(3)
1377  jth1(4) = -jth1(4)
1378  CASE DEFAULT
1379  jth = -jth
1380  jth1 = -jth1
1381  END SELECT
1382  !
1383  ist = (jiq-1)*4 + (jof-1)*2 + 1
1384  ast1(ist,nqa+jqr,ikd) = ifr
1385  ast2(ist,nqa+jqr,ikd) = jth(jiq)
1386  qst5(ist,nqa+jqr,ikd) = wf * ( wfroff * wth(jiq) )
1387  qst6(ist,nqa+jqr,ikd) = wf * ( wfroff * wth(jiq) )**2
1388  ist = ist + 1
1389  ast1(ist,nqa+jqr,ikd) = ifr
1390  ast2(ist,nqa+jqr,ikd) = jth1(jiq)
1391  qst5(ist,nqa+jqr,ikd) = wf * ( wfroff * wth1(jiq) )
1392  qst6(ist,nqa+jqr,ikd) = wf * ( wfroff * wth1(jiq) )**2
1393  !
1394  ! ... End loop 3.f
1395  !
1396  END DO
1397  !
1398  ! ... End loop 3.e
1399  !
1400  END DO
1401  !
1402  ! ... End loop 3.d
1403  !
1404  END DO
1405  !
1406 #ifdef W3_T3
1407  DO jqr=1, jj
1408  WRITE (ndst,9037) ikd, nqa+jqr
1409  DO ist=1, 16
1410  WRITE (ndst,9038) ist, ast1(ist,nqa+jqr,ikd), &
1411  ast2(ist,nqa+jqr,ikd), &
1412  qst5(ist,nqa+jqr,ikd), &
1413  qst6(ist,nqa+jqr,ikd)
1414  END DO
1415  END DO
1416 #endif
1417  !
1418  ! ... End loop 3.b
1419  !
1420  nqa = nqa + jj
1421  !
1422  END DO
1423  !
1424  ! ... End loop 3.a
1425  !
1426  END DO
1427  !
1428  ! 3.g Finalize storage
1429  !
1430  qst4 = ast1*nthexp + ast2
1431  !
1432  IF ( nthmax .LT. nthmx2 ) GOTO 810
1433  IF ( nqa .NE. SIZE(ast1(1,:,1)) ) GOTO 811
1434  !
1435  DEALLOCATE ( ast1, ast2 )
1436  !
1437  RETURN
1438  !
1439  ! Error escape locations
1440  !
1441 800 CONTINUE
1442  WRITE (ndse,1000) lammax, delthm
1443  CALL extcde ( 1000 )
1444  !
1445 801 CONTINUE
1446  WRITE (ndse,1001) off12, off34
1447  CALL extcde ( 1001 )
1448  !
1449 802 CONTINUE
1450  WRITE (ndse,1002) off12, off34, snlt(iq), &
1451  minlam(off12,snlt(iq)), maxlam(off12,snlt(iq))
1452  CALL extcde ( 1002 )
1453  !
1454 810 CONTINUE
1455  WRITE (ndse,1010) nthmax, nthmx2
1456  CALL extcde ( 1010 )
1457  !
1458 811 CONTINUE
1459  WRITE (ndse,1011) nqa, SIZE(ast1(1,:,1))
1460  CALL extcde ( 1011 )
1461  !
1462  RETURN
1463  !
1464  ! Formats
1465  !
1466 1000 FORMAT (/' *** WAVEWATCH-III ERROR IN INSNL3 :'/ &
1467  ' PARAMETER OUT OF RANGE '/ &
1468  ' LAMMAX, DELTHM :', 2e12.4/)
1469 1001 FORMAT (/' *** WAVEWATCH-III ERROR IN INSNL3 :'/ &
1470  ' PARAMETER OUT OF RANGE '/ &
1471  ' MU, LAMBDA :', 2e12.4/)
1472 1002 FORMAT (/' *** WAVEWATCH-III ERROR IN INSNL3 :'/ &
1473  ' PARAMETER OUT OF RANGE '/ &
1474  ' MU, LAMBDA, TH12 :',3e12.4/ &
1475  ' LAMBDA RANGE :',2e12.4)
1476 1010 FORMAT (/' *** WAVEWATCH-III ERROR IN INSNL3 :'/ &
1477  ' NTHMAX LESS THAN NTHMX2 :', 2i8/)
1478 1011 FORMAT (/' *** WAVEWATCH-III ERROR IN INSNL3 :'/ &
1479  ' NQA INCONSISTENT :', 2i8/)
1480  !
1481 #ifdef W3_T
1482 9010 FORMAT (/' TEST INSNL3: NKD, KDMIN/MAX/X : ',i8,3f10.4)
1483 #endif
1484  !
1485 #ifdef W3_T1
1486 9020 FORMAT (/' TEST INSNL3: IKD, IQ, KD, F, D: ',2i8,2f10.4,f10.2)
1487 9021 FORMAT (/' TEST INSNL3: TH12 : ',3x,f8.2/ &
1488  ' OFF12, OFF34 : ',3x,2f8.2/ &
1489  ' CD, CS : ',3x,2e10.2)
1490 9022 FORMAT ( ' ANGLES (DEGR) : ',1x,4f8.2)
1491 9023 FORMAT ( ' FREQUENCY IND. : ',1x,3i4,2f6.2)
1492 9024 FORMAT ( ' : ',1x,3i4,2f6.2)
1493 9025 FORMAT ( ' DIRECTIONAL IND. : ',1x,3i4,2f6.2)
1494 #endif
1495 #ifdef W3_T
1496 9026 FORMAT ( ' TEST INSNL3: FILLING FIRST DATA TABLES :'/ &
1497  ' NQA AND MAXIMUM : ',2i8/ &
1498  ' NQD AND NQS : ',2i8)
1499 9027 FORMAT ( ' NFR, MIN/MAX/CUT : ',4i8/ &
1500  ' NTH, MIN/MAX/EXP : ',4i8)
1501 #endif
1502  !
1503 #ifdef W3_T2
1504 9030 FORMAT (/' TEST INSNL3: IKD, IQ, KD, F, D: ',2i8,2f10.4,f10.2)
1505 9031 FORMAT (/' TEST INSNL3: TH12 : ',3x,f8.2/ &
1506  ' OFF12, OFF34 : ',3x,2f8.2)
1507 9032 FORMAT ( ' FREQUENCY IND. : ',1x,3i4,2f6.2)
1508 9033 FORMAT ( ' : ',1x,3i4,2f6.2)
1509 9034 FORMAT ( ' J,J,J, W, SIn : ',1x,3i4,2f6.2)
1510 9035 FORMAT ( ' ANGLES (DEGR) : ',3x,4f8.2)
1511 9036 FORMAT ( ' DIRECTIONAL IND. : ',1x,3i4,2f6.2)
1512 #endif
1513 #ifdef W3_T3
1514 9037 FORMAT (/' TEST INSNL3: STORAGE ARRAYS FOR IKD, IQA =',2i6)
1515 9038 FORMAT (23x,3i4,3f8.3)
1516 #endif
1517  !/
1518  !/ Embedded subroutines
1519  !/
1520  CONTAINS
1521  !/ ------------------------------------------------------------------- /
1522 
1533  REAL FUNCTION MINLAM ( MU, THETA )
1534  !/
1535  !/ +-----------------------------------+
1536  !/ | WAVEWATCH-III NOAA/NCEP |
1537  !/ | H. L. Tolman |
1538  !/ | FORTRAN 90 |
1539  !/ | Last update : 28-Jan-2004 |
1540  !/ +-----------------------------------+
1541  !/
1542  !/ 28-Jan-2009 : Origination.
1543  !/
1544  ! 1. Purpose :
1545  !
1546  ! Calculate minimum allowed lambda for quadruplet configuration.
1547  !
1548  ! 3. Parameters :
1549  !
1550  ! Parameter list
1551  ! ----------------------------------------------------------------
1552  ! MU, THETA Real Quadruplet parameters, theta in degree.
1553  ! ----------------------------------------------------------------
1554  !
1555  ! 10. Source code :
1556  !
1557  !/ ------------------------------------------------------------------- /
1558  IMPLICIT NONE
1559  !/
1560  !/ Parameter list
1561  !/
1562  REAL, INTENT(IN) :: MU, THETA
1563  !/
1564  !/ Local parameters
1565  !/
1566  REAL :: MULOC, THETAR, BB, AUX
1567  !/
1568  !/ ------------------------------------------------------------------- /
1569  !/
1570  IF ( theta .LT. 0. ) THEN
1571  minlam = 0.
1572  ELSE
1573  muloc = max( 0. , min( 1., mu ) )
1574  thetar = theta * atan(1.) / 45.
1575  bb = (1.+muloc)**4 + (1.-muloc)**4 + &
1576  2. * (1.+muloc)**2 * (1.-muloc)**2 * cos(thetar)
1577  bb = sqrt( max( bb , 0. ) )
1578  aux = max( 0. , 0.5*bb-1. )
1579  minlam = sqrt( aux )
1580  END IF
1581  !
1582  RETURN
1583  !/
1584  !/ End of MINLAM ----------------------------------------------------- /
1585  !/
1586  END FUNCTION minlam
1587  !/ ------------------------------------------------------------------- /
1588 
1601  REAL FUNCTION MAXLAM ( MU, THETA )
1602  !/
1603  !/ +-----------------------------------+
1604  !/ | WAVEWATCH-III NOAA/NCEP |
1605  !/ | H. L. Tolman |
1606  !/ | FORTRAN 90 |
1607  !/ | Last update : 28-Jan-2004 |
1608  !/ +-----------------------------------+
1609  !/
1610  !/ 28-Jan-2009 : Origination.
1611  !/
1612  ! 1. Purpose :
1613  !
1614  ! Calculate minimum allowed lambda for quadruplet configuration.
1615  !
1616  ! 3. Parameters :
1617  !
1618  ! Parameter list
1619  ! ----------------------------------------------------------------
1620  ! MU, THETA Real Quadruplet parameters, theta in degree.
1621  ! ----------------------------------------------------------------
1622  !
1623  ! 10. Source code :
1624  !
1625  !/ ------------------------------------------------------------------- /
1626  IMPLICIT NONE
1627  !/
1628  !/ Parameter list
1629  !/
1630  REAL, INTENT(IN) :: MU, THETA
1631  !/
1632  !/ Local parameters
1633  !/
1634  REAL :: MULOC, THETAR, BB, AUX
1635  !/
1636  !/ ------------------------------------------------------------------- /
1637  !/
1638  IF ( theta .LT. 0. ) THEN
1639  maxlam = 0.5
1640  ELSE
1641  muloc = max( 0. , min( 1., mu ) )
1642  thetar = theta * atan(1.) / 45.
1643  bb = (1.+muloc)**4 + (1.-muloc)**4 + &
1644  2. * (1.+muloc)**2 * (1.-muloc)**2 * cos(thetar)
1645  bb = sqrt( max( bb , 0. ) )
1646  maxlam = 0.25 * bb
1647  END IF
1648  !
1649  RETURN
1650  !/
1651  !/ End of MAXLAM ----------------------------------------------------- /
1652  !/
1653  END FUNCTION maxlam
1654  !/
1655  !/ End of INSNL3 ----------------------------------------------------- /
1656  !/

References delthm, constants::dera, w3gdatmd::dth, w3servmd::extcde(), w3gdatmd::frq, constants::grav, w3gdatmd::igrid, lammax, maxlam(), minlam(), w3gdatmd::mpars, w3odatmd::ndse, w3odatmd::ndst, w3gdatmd::nfrcut, w3gdatmd::nfrmax, w3gdatmd::nfrmin, w3gdatmd::nk, w3gdatmd::nqa, w3gdatmd::nspmax, w3gdatmd::nspmin, w3gdatmd::nspmx2, w3gdatmd::nth, w3gdatmd::nthexp, w3gdatmd::nthmax, w3gdatmd::qst1, w3gdatmd::qst2, w3gdatmd::qst3, w3gdatmd::qst4, w3gdatmd::qst5, w3gdatmd::qst6, constants::rade, w3gdatmd::sig, w3gdatmd::snlcd, w3gdatmd::snlcs, w3gdatmd::snll, w3gdatmd::snlm, w3gdatmd::snlnq, w3gdatmd::snlt, w3servmd::strace(), constants::tpiinv, w3dispmd::wavnu2(), w3gdatmd::xfr, and w3gdatmd::xsi.

Referenced by w3iogrmd::w3iogr().

◆ w3snl3()

subroutine w3snl3md::w3snl3 ( real, dimension(nth,nfr), intent(in)  A,
real, dimension(nfr), intent(in)  CG,
real, dimension(nfr), intent(in)  WN,
real, intent(in)  DEPTH,
real, dimension(nth,nfr), intent(out)  S,
real, dimension(nth,nfr), intent(out)  D 
)

Multiple Discrete Interaction Parameterization for arbitrary depths with generalized quadruplet layout.

This is a direct implementation of the Discrete Interaction Paramterization (DIA) with multiple representative quadruplets (MDIA) for arbitrary water depths.

The outer loop of the code is over quadruplet realizations, which implies two realizations for a conventional quadruplet definitions and four for extended definitions (with rescaling of the contants for consistency). Within this loop the compu- tations are performed in two stages. First, interactions contributions are computed for the entire spectral space, second all contributions are combined into the actual inter- actions and diagonal contributions.

Arbitrary depths are addressed by generating a lookup table for the relative depth. These tables are used for each discrete frequency separately. Efficient memory usages requires relative addressing to reduce the size of the lookup tables. To use this the spectral space is expanded to higher and lower frequencies as well as directional space is expanded/volded. This is done for the input (pseudo-) spectrum (action spectrum devided by the wavenumber) to determine spectral densities at the quadruplet components, and the spectral space describing individual contri- butions before they are combined into the actual interactions.

Parameters
[in]AAction spectrum A(ITH,IK) as a function of direction (rad) and wavenumber.
[in]CGGroup velocities (dimension NK).
[in]WNWavenumbers (dimension NK).
[in]DEPTHWater depth in meters.
[out]SSource term.
[out]DDiagonal term of derivative.
Author
H. L. Tolman
Date
01-Dec-2009

Definition at line 181 of file w3snl3md.F90.

181  !/
182  !/ +-----------------------------------+
183  !/ | WAVEWATCH-III NOAA/NCEP |
184  !/ | H. L. Tolman |
185  !/ | FORTRAN 90 |
186  !/ | Last update : 01-Dec-2009 |
187  !/ +-----------------------------------+
188  !/
189  !/ 21-Jul-2008 : Origination as NLX option. ( version 3.13 )
190  !/ 25-Aug-2009 : Conversion to F(f,theta) form. ( version 3.13 )
191  !/ 01-Dec-2009 : Bug fix frequency filtering. ( version 3.13 )
192  !/
193  ! 1. Purpose :
194  !
195  ! Multiple Discrete Interaction Parameterization for arbitrary
196  ! depths with generalized quadruplet layout.
197  !
198  ! 2. Method :
199  !
200  ! This is a direct implementation of the Discrete Interaction
201  ! Paramterization (DIA) with multiple representative quadruplets
202  ! (MDIA) for arbitrary water depths.
203  !
204  ! The outer loop of the code is over quadruplet realizations,
205  ! which implies two realizations for a conventional quadruplet
206  ! definitions and four for extended definitions (with rescaling
207  ! of the contants for consistency). Within this loop the compu-
208  ! tations are performed in two stages. First, interactions
209  ! contributions are computed for the entire spectral space,
210  ! second all contributions are combined into the actual inter-
211  ! actions and diagonal contributions.
212  !
213  ! Arbitrary depths are addressed by generating a lookup table
214  ! for the relative depth. These tables are used for each discrete
215  ! frequency separately. Efficient memory usages requires relative
216  ! addressing to reduce the size of the lookup tables. To use this
217  ! the spectral space is expanded to higher and lower frequencies
218  ! as well as directional space is expanded/volded. This is done
219  ! for the input (pseudo-) spectrum (action spectrum devided by the
220  ! wavenumber) to determine spectral densities at the quadruplet
221  ! components, and the spectral space describing individual contri-
222  ! butions before they are combined into the actual interactions.
223  !
224  ! 3. Parameters :
225  !
226  ! Parameter list
227  ! ----------------------------------------------------------------
228  ! A R.A. I Action spectrum A(ITH,IK) as a function of
229  ! direction (rad) and wavenumber.
230  ! CG R.A. I Group velocities (dimension NK).
231  ! WN R.A. I Wavenumbers (dimension NK).
232  ! DEPTH Real I Water depth in meters.
233  ! S R.A. O Source term.
234  ! D R.A. O Diagonal term of derivative.
235  ! ----------------------------------------------------------------
236  !
237  ! Variables describing the expanded frequency space from the
238  ! dynamic storage in w3gdatmd.
239  !
240  ! Name Type Scope Description
241  ! ----------------------------------------------------------------
242  ! NFR Int. Public Number of frequencies or wavenumbers
243  ! in discrete spectral space (NFR=>NK).
244  ! NFRMIN Int. Public Minimum discrete frequency in the
245  ! expanded frequency space.
246  ! NFRMAX Int. Public Idem maximum for first part.
247  ! NFRCUT Int. Public Idem maximum for second part.
248  ! NTHMAX Int. Public Extension of directional space.
249  ! NTHEXP Int Public Number of bins in extended dir. space.
250  ! NSPMIN, NSPMAX, NSPMX2
251  ! Int. Public 1D spectral space range.
252  ! FRQ R.A. Public Expanded frequency range (Hz).
253  ! XSI R.A. Public Expanded frequency range (rad/s).
254  ! ----------------------------------------------------------------
255  !
256  ! Variables describing lookup tables.
257  !
258  ! Name Type Scope Description
259  ! ----------------------------------------------------------------
260  ! NQA Int. Public Number of actual quadruplets.
261  ! QST1 I.A. Public Spectral offsets for compuation of
262  ! quadruplet spectral desnities.
263  ! QST2 R.A. Public Idem weights.
264  ! QST3 R.A. Public Norm. factors in product term and
265  ! in diagonal strength.
266  ! QST4 I.A. Public Spectral offsets for combining of
267  ! interactions and diagonal.
268  ! QST5 R.A. Public Idem weights for interactions.
269  ! QST6 R.A. Public Idem weights for diagonal.
270  ! ----------------------------------------------------------------
271  !
272  ! Variables describing model setup.
273  !
274  ! Name Type Scope Description
275  ! ----------------------------------------------------------------
276  ! SNLMSC Real Public Tuning power 'deep' scaling.
277  ! SNLNSC Real Public Tuning power 'shallow' scaling.
278  ! SNLSFD Real Public 'Deep' nondimensional filer freq.
279  ! SNLSFS Real Public 'Shallow' nondimensional filer freq.
280  ! ----------------------------------------------------------------
281  !
282  ! 4. Subroutines used :
283  !
284  ! Name Type Module Description
285  ! ----------------------------------------------------------------
286  ! STRACE Subr. W3SERVMD Subroutine tracing.
287  ! ----------------------------------------------------------------
288  !
289  ! 5. Called by :
290  !
291  ! Name Type Module Description
292  ! ----------------------------------------------------------------
293  ! W3SRCE Subr. W3SRCEMD Source term integration.
294  ! W3EXPO Subr. N/A Point output post-processor.
295  ! GXEXPO Subr. N/A GrADS point output post-processor.
296  ! ----------------------------------------------------------------
297  !
298  ! 6. Error messages :
299  !
300  ! None.
301  !
302  ! 7. Remarks :
303  !
304  ! - Note that this code uses explicit unroling of potential loop
305  ! structures for optimization purposes.
306  ! - Normalization with respect to the number of quadruplets is
307  ! included in the proportionality constant.
308  ! - Note that the outer loop in the routine considers one actual
309  ! quadruplet realization per loop cycle. For the traditional
310  ! quadruplet layout two realizations occure, for the expanded
311  ! four realizations occur. For consistency, strength of a
312  ! traditional layout is therefore doubled.
313  ! - 1D representation is used of 2D spectral space for optimization
314  ! purposes.
315  ! - Contributions are first computed in the convetional spectral
316  ! space and are then expancded "in place" into the expanded
317  ! spectral space in EXPND2.
318  !
319  ! 8. Structure :
320  !
321  ! See source code.
322  !
323  ! 9. Switches :
324  !
325  ! !/S Enable subroutine tracing.
326  !
327  ! 10. Source code :
328  !
329  !/ ------------------------------------------------------------------- /
330  USE constants
331  USE w3gdatmd, ONLY: nfr => nk, nth, sig, fachfe, facti1, facti2,&
332  nfrmin, nfrmax, nfrcut, nthmax, nthexp, &
333  nspmin, nspmax, nspmx2, frq, xsi, nqa, &
334  qst1, qst2, qst3, qst4, qst5, qst6, snlmsc, &
335  snlnsc, snlsfd, snlsfs
336  USE w3odatmd, ONLY: ndse, ndst
337  !
338  USE w3servmd, ONLY: extcde
339  USE w3dispmd, ONLY: wavnu1, wavnu3
340 #ifdef W3_S
341  USE w3servmd, ONLY: strace
342 #endif
343  !/
344  IMPLICIT NONE
345  !/
346  !/ ------------------------------------------------------------------- /
347  !/ Parameter list
348  !/
349  REAL, INTENT(IN) :: A(NTH,NFR), CG(NFR), WN(NFR), DEPTH
350  REAL, INTENT(OUT) :: S(NTH,NFR), D(NTH,NFR)
351  !/
352  !/ ------------------------------------------------------------------- /
353  !/ Local parameters
354  !/
355  INTEGER :: IFR, IERR, IKD, JKD(NFRCUT), IQA, IF1MIN, &
356  IF1MAX, IF2MIN, IF2MAX, ISP0, ISPX0, ITH, &
357  ISP, ISPX
358 #ifdef W3_S
359  INTEGER, SAVE :: IENT = 0
360 #endif
361  INTEGER :: LQST1(16), LQST4(16)
362  REAL :: XSITLN, SIT, FPROP, FQ1, FQ2, FQ3, FQ4, &
363  AUX1, AUX2
364  REAL :: XWN(NFRMAX), XCG(NFRMAX), SCALE1(NFRCUT), &
365  SCALE2(NFRCUT), LQST2(16), FACT(6), &
366  LQST5(16), LQST6(16)
367  REAL :: UE(NSPMIN:NSPMAX), DSB(NSPMIN:NSPMX2), &
368  DD1(NSPMIN:NSPMX2), DD2(NSPMIN:NSPMX2), &
369  DD3(NSPMIN:NSPMX2), DD4(NSPMIN:NSPMX2)
370  !/
371  !/ ------------------------------------------------------------------- /
372  !/
373 #ifdef W3_S
374  CALL strace (ient, 'W3SNL3')
375 #endif
376  !
377  ! 1. Initialization ------------------------------------------------- *
378  ! 1.a Constants and arrays
379  !
380  xsitln = log(xsit)
381  !
382  s = 0.
383  d = 0.
384  ! DSB = 0.
385  ! DD1 = 0.
386  ! DD2 = 0.
387  ! DD3 = 0.
388  ! DD4 = 0.
389  !
390  ! 1.a Extended frequency range
391  !
392  xwn(1:nfr) = wn
393  xcg(1:nfr) = cg
394  !
395  DO ifr = nfr+1, nfrmax
396 #ifdef W3_PDLIB
397  CALL wavnu3(xsi(ifr), depth, xwn(ifr), xcg(ifr))
398 #else
399  CALL wavnu1(xsi(ifr), depth, xwn(ifr), xcg(ifr))
400 #endif
401  END DO
402  !
403  ! 1.b Expanded pseudo spetrum
404  !
405  CALL expand ( ue )
406  !
407  ! 1.c Set up scaling functions
408  !
409  aux1 = 1. / ( tpi**11 * grav**(4.-snlmsc) )
410  aux2 = grav**2 / tpi**11
411  !
412  DO ifr=1, nfrcut
413  scale1(ifr) = aux1 * xwn(ifr)**(4.+snlmsc) * &
414  xsi(ifr)**(13.-2.*snlmsc) / xcg(ifr)**2
415  scale2(ifr) = aux2 * xwn(ifr)**11 * &
416  (xwn(ifr)*depth)**snlnsc / xcg(ifr)
417  END DO
418  !
419  ! 1.d Set up depth scaling counters
420  !
421  DO ifr=1, nfrcut
422  sit = xsi(ifr) * sqrt(depth/grav)
423  ikd = 1 + nint( ( log(sit) - log(sitmin) ) / xsitln )
424  jkd(ifr) = max( 1 , min(ikd,nkd) )
425  END DO
426  !
427  ! 2. Base loop over quadruplet realizations ------------------------- *
428  !
429  DO iqa=1 , nqa
430  !
431  ! 3. Obtain quadruplet energies for all spectral bins --------------- *
432  ! 3.a Set frequency ranges
433  !
434  aux1 = qst3(5,iqa,1)
435  aux2 = qst3(6,iqa,1)
436  !
437  if1min = 1
438  if1max = nfrcut
439  if2min = 1
440  if2max = nfr
441  !
442  IF ( aux1 .LE. 0. .AND. aux2 .LE. 0. ) THEN
443  !
444  cycle
445  !
446  ELSE IF ( aux2 .LE. 0. ) THEN
447  !
448  sit = snlsfd * sqrt(grav/depth)
449  ifr = nint( facti2 + facti1*log(sit) )
450  IF ( ifr .GT. nfr ) cycle
451  !
452  IF ( ifr .GT. 1 ) THEN
453  if1min = max( 1 , ifr )
454  if2min = max( 1 , if1min + nfrmin )
455  dsb(1:(if1min-1)*nth) = 0.
456  dd1(1:(if1min-1)*nth) = 0.
457  dd2(1:(if1min-1)*nth) = 0.
458  dd3(1:(if1min-1)*nth) = 0.
459  dd4(1:(if1min-1)*nth) = 0.
460  END IF
461  !
462  ELSE IF ( aux1 .LE. 0. ) THEN
463  !
464  sit = snlsfs * sqrt(grav/depth)
465  ifr = nint( facti2 + facti1*log(sit) )
466  IF ( ifr .LT. 1 ) cycle
467  !
468  IF ( ifr .LT. nfrcut ) THEN
469  if1max = min( nfrcut, ifr )
470  ! IF2MAX = NFR
471  dsb(if1max*nth+1:nfrcut*nth) = 0.
472  dd1(if1max*nth+1:nfrcut*nth) = 0.
473  dd2(if1max*nth+1:nfrcut*nth) = 0.
474  dd3(if1max*nth+1:nfrcut*nth) = 0.
475  dd4(if1max*nth+1:nfrcut*nth) = 0.
476  END IF
477  !
478  END IF
479  !
480  ! 3.b Loop over frequencies
481  !
482  DO ifr=if1min, if1max
483  !
484  ! 3.c Find discrete depths
485  !
486  ikd = jkd(ifr)
487  !
488  ! 3.d Get offsets and weights
489  !
490  lqst1 = qst1(:,iqa,ikd)
491  lqst2 = qst2(:,iqa,ikd)
492  fact = qst3(:,iqa,ikd)
493  fact(1:4) = fact(1:4) * xcg(ifr) / ( xwn(ifr) *xsi(ifr) )
494  fprop = scale1(ifr)*fact(5) + scale2(ifr)*fact(6)
495  !
496  ! 3.e Loop over directions
497  !
498  isp0 = (ifr-1)*nth
499  ispx0 = (ifr-1)*nthexp
500  !
501  DO ith=1, nth
502  !
503  isp = isp0 + ith
504  ispx = ispx0 + ith
505  !
506  fq1 = ( ue(ispx+lqst1( 1)) * lqst2( 1) + &
507  ue(ispx+lqst1( 2)) * lqst2( 2) + &
508  ue(ispx+lqst1( 3)) * lqst2( 3) + &
509  ue(ispx+lqst1( 4)) * lqst2( 4) ) * fact(1)
510  fq2 = ( ue(ispx+lqst1( 5)) * lqst2( 5) + &
511  ue(ispx+lqst1( 6)) * lqst2( 6) + &
512  ue(ispx+lqst1( 7)) * lqst2( 7) + &
513  ue(ispx+lqst1( 8)) * lqst2( 8) ) * fact(2)
514  fq3 = ( ue(ispx+lqst1( 9)) * lqst2( 9) + &
515  ue(ispx+lqst1(10)) * lqst2(10) + &
516  ue(ispx+lqst1(11)) * lqst2(11) + &
517  ue(ispx+lqst1(12)) * lqst2(12) ) * fact(3)
518  fq4 = ( ue(ispx+lqst1(13)) * lqst2(13) + &
519  ue(ispx+lqst1(14)) * lqst2(14) + &
520  ue(ispx+lqst1(15)) * lqst2(15) + &
521  ue(ispx+lqst1(16)) * lqst2(16) ) * fact(4)
522  !
523  aux1 = fq1 * fq2 * ( fq3 + fq4 )
524  aux2 = fq3 * fq4 * ( fq1 + fq2 )
525  dsb(isp) = fprop * ( aux1 - aux2 )
526  !
527  aux1 = fq3 + fq4
528  aux2 = fq3 * fq4
529  dd1(isp) = fprop * fact(1) * ( fq2 * aux1 - aux2 )
530  dd2(isp) = fprop * fact(2) * ( fq1 * aux1 - aux2 )
531  !
532  aux1 = fq1 + fq2
533  aux2 = fq1 * fq2
534  dd3(isp) = fprop * fact(3) * ( aux2 - fq4*aux1 )
535  dd4(isp) = fprop * fact(4) * ( aux2 - fq3*aux1 )
536  !
537  ! ... End loop 3.e
538  !
539  END DO
540  !
541  ! ... End loop 3.b
542  !
543  END DO
544  !
545  ! 3.e Expand arrays
546  !
547  CALL expnd2 ( dsb(1:nth*nfrcut), dsb )
548  CALL expnd2 ( dd1(1:nth*nfrcut), dd1 )
549  CALL expnd2 ( dd2(1:nth*nfrcut), dd2 )
550  CALL expnd2 ( dd3(1:nth*nfrcut), dd3 )
551  CALL expnd2 ( dd4(1:nth*nfrcut), dd4 )
552  !
553  ! 4. Put it all together -------------------------------------------- *
554  ! 4.a Loop over frequencies
555  !
556  DO ifr=if2min, if2max
557  !
558  ! 4.b Find discrete depths and storage
559  !
560  ikd = jkd(ifr)
561  !
562  ! 4.c Get offsets and weights
563  !
564  lqst4 = qst4(:,iqa,ikd)
565  lqst5 = qst5(:,iqa,ikd)
566  lqst6 = qst6(:,iqa,ikd)
567  !
568  ! 4.d Loop over directions
569  !
570  ispx0 = (ifr-1)*nthexp
571  !
572  DO ith=1, nth
573  !
574  ispx = ispx0 + ith
575  !
576  s(ith,ifr) = s(ith,ifr) + dsb(ispx+lqst4( 1)) * lqst5( 1) &
577  + dsb(ispx+lqst4( 2)) * lqst5( 2) &
578  + dsb(ispx+lqst4( 3)) * lqst5( 3) &
579  + dsb(ispx+lqst4( 4)) * lqst5( 4) &
580  + dsb(ispx+lqst4( 5)) * lqst5( 5) &
581  + dsb(ispx+lqst4( 6)) * lqst5( 6) &
582  + dsb(ispx+lqst4( 7)) * lqst5( 7) &
583  + dsb(ispx+lqst4( 8)) * lqst5( 8) &
584  + dsb(ispx+lqst4( 9)) * lqst5( 9) &
585  + dsb(ispx+lqst4(10)) * lqst5(10) &
586  + dsb(ispx+lqst4(11)) * lqst5(11) &
587  + dsb(ispx+lqst4(12)) * lqst5(12) &
588  + dsb(ispx+lqst4(13)) * lqst5(13) &
589  + dsb(ispx+lqst4(14)) * lqst5(14) &
590  + dsb(ispx+lqst4(15)) * lqst5(15) &
591  + dsb(ispx+lqst4(16)) * lqst5(16)
592  !
593  d(ith,ifr) = d(ith,ifr) + dd1(ispx+lqst4( 1)) * lqst6( 1) &
594  + dd1(ispx+lqst4( 2)) * lqst6( 2) &
595  + dd1(ispx+lqst4( 3)) * lqst6( 3) &
596  + dd1(ispx+lqst4( 4)) * lqst6( 4) &
597  + dd2(ispx+lqst4( 5)) * lqst6( 5) &
598  + dd2(ispx+lqst4( 6)) * lqst6( 6) &
599  + dd2(ispx+lqst4( 7)) * lqst6( 7) &
600  + dd2(ispx+lqst4( 8)) * lqst6( 8) &
601  + dd3(ispx+lqst4( 9)) * lqst6( 9) &
602  + dd3(ispx+lqst4(10)) * lqst6(10) &
603  + dd3(ispx+lqst4(11)) * lqst6(11) &
604  + dd3(ispx+lqst4(12)) * lqst6(12) &
605  + dd4(ispx+lqst4(13)) * lqst6(13) &
606  + dd4(ispx+lqst4(14)) * lqst6(14) &
607  + dd4(ispx+lqst4(15)) * lqst6(15) &
608  + dd4(ispx+lqst4(16)) * lqst6(16)
609  !
610  ! ... End loop 4.d
611  !
612  END DO
613  !
614  ! ... End loop 4.a
615  !
616  END DO
617  !
618  ! ... End of loop 2.
619  !
620  END DO
621  !
622  ! 5. Convert back to wave action ------------------------------------ *
623  !
624  DO ifr=if2min, if2max
625  s(:,ifr) = s(:,ifr) / xsi(ifr) * xcg(ifr) * tpiinv
626  END DO
627  !
628  RETURN
629  !/
630  !/ Embedded subroutines
631  !/
632  CONTAINS
633  !/ ------------------------------------------------------------------- /
634 
643  SUBROUTINE expand ( SPEC )
644  !/
645  !/ +-----------------------------------+
646  !/ | WAVEWATCH-III NOAA/NCEP |
647  !/ | H. L. Tolman |
648  !/ | FORTRAN 90 |
649  !/ | Last update : 21-Aug-2009 |
650  !/ +-----------------------------------+
651  !/
652  !/ 03-Jul-2008 : Origination. ( version 3.13 )
653  !/ 21-Aug-2009 : Conversion to F(f,theta) form. ( version 3.13 )
654  !/
655  ! 1. Purpose :
656  !
657  ! Expand spectrum, subroutine used to simplify addressing.
658  !
659  ! 3. Parameters :
660  !
661  ! Parameter list
662  ! ----------------------------------------------------------------
663  ! SPEC R.A. O Expanded spectrum.
664  ! ----------------------------------------------------------------
665  !
666  ! 10. Source code :
667  !
668  !/ ------------------------------------------------------------------- /
669  IMPLICIT NONE
670  !/
671  !/ Parameter list
672  !/
673  REAL, INTENT(OUT) :: SPEC(1-NTHMAX:NTH+NTHMAX,NFRMIN:NFRMAX)
674  !/
675  !/ Local parameters
676  !/
677  INTEGER :: IFR, ITH
678  !/
679  !/ ------------------------------------------------------------------- /
680  !
681  spec(:,nfrmin:0) = 0.
682  !
683  spec(1:nth,1:nfr) = a * tpi
684  !
685  DO ifr=1, nfr
686  spec(1:nth,ifr) = spec(1:nth,ifr) * xsi(ifr) / xcg(ifr)
687  END DO
688  !
689  DO ifr=nfr+1, nfrmax
690  spec(1:nth,ifr) = spec(1:nth,ifr-1) * fachfe
691  END DO
692  !
693  DO ith=1, nthmax
694  spec(nth+ith,1:nfrmax) = spec( ith ,1:nfrmax)
695  spec( 1 -ith,1:nfrmax) = spec(nth+1-ith,1:nfrmax)
696  END DO
697  !
698  RETURN
699  !/
700  !/ End of EXPAND ----------------------------------------------------- /
701  !/
702  END SUBROUTINE expand
703  !/ ------------------------------------------------------------------- /
704 
716  SUBROUTINE expnd2 ( ARIN, AROUT )
717  !/
718  !/ +-----------------------------------+
719  !/ | WAVEWATCH-III NOAA/NCEP |
720  !/ | H. L. Tolman |
721  !/ | FORTRAN 90 |
722  !/ | Last update : 16-Jul-2008 |
723  !/ +-----------------------------------+
724  !/
725  !/ 16-Jul-2008 : Origination. ( version 3.13 )
726  !/
727  ! 1. Purpose :
728  !
729  ! Expand spectrum to simplify indirect addressing.
730  ! Done 'in place' with temporary array ( ARIN = AROUT )
731  !
732  ! 3. Parameters :
733  !
734  ! Parameter list
735  ! ----------------------------------------------------------------
736  ! SPIN R.A. I Input array.
737  ! SPOUT R.A. I Output array.
738  ! ----------------------------------------------------------------
739  !
740  ! 10. Source code :
741  !
742  !/ ------------------------------------------------------------------- /
743  IMPLICIT NONE
744  !/
745  !/ Parameter list
746  !/
747  REAL, INTENT(IN) :: ARIN(NTH,NFRCUT)
748  REAL, INTENT(OUT) :: AROUT(1-NTHMAX:NTH+NTHMAX,NFRMIN:NFRCUT)
749  !/
750  !/ Local parameters
751  !/
752  INTEGER :: IFR, ITH
753  REAL :: TEMP(NTH,NFRCUT)
754  !/
755  !/ ------------------------------------------------------------------- /
756  !
757  temp = arin
758  !
759  arout(:,nfrmin:0) = 0.
760  !
761  arout(1:nth,1:nfrcut) = temp
762  !
763  DO ith=1, nthmax
764  arout(nth+ith,1:nfrcut) = arout( ith ,1:nfrcut)
765  arout( 1 -ith,1:nfrcut) = arout(nth+1-ith,1:nfrcut)
766  END DO
767  !
768  RETURN
769  !/
770  !/ End of EXPND2 ----------------------------------------------------- /
771  !/
772  END SUBROUTINE expnd2
773  !/
774  !/ End of W3SNL3 ----------------------------------------------------- /
775  !/

References expand(), expnd2(), w3servmd::extcde(), w3gdatmd::fachfe, w3gdatmd::facti1, w3gdatmd::facti2, w3gdatmd::frq, constants::grav, w3odatmd::ndse, w3odatmd::ndst, w3gdatmd::nfrcut, w3gdatmd::nfrmax, w3gdatmd::nfrmin, w3gdatmd::nk, w3gdatmd::nqa, w3gdatmd::nspmax, w3gdatmd::nspmin, w3gdatmd::nspmx2, w3gdatmd::nth, w3gdatmd::nthexp, w3gdatmd::nthmax, w3gdatmd::qst1, w3gdatmd::qst2, w3gdatmd::qst3, w3gdatmd::qst4, w3gdatmd::qst5, w3gdatmd::qst6, w3gdatmd::sig, w3gdatmd::snlmsc, w3gdatmd::snlnsc, w3gdatmd::snlsfd, w3gdatmd::snlsfs, w3servmd::strace(), constants::tpi, constants::tpiinv, w3dispmd::wavnu1(), w3dispmd::wavnu3(), and w3gdatmd::xsi.

Referenced by gxexpo(), w3exnc(), w3expo(), and w3srcemd::w3srce().

Variable Documentation

◆ delthm

real, parameter, public w3snl3md::delthm = 90.

Definition at line 132 of file w3snl3md.F90.

132  REAL, PUBLIC, PARAMETER :: DELTHM = 90.

Referenced by insnl3().

◆ lammax

real, parameter, public w3snl3md::lammax = 0.49999

Definition at line 131 of file w3snl3md.F90.

131  REAL, PUBLIC, PARAMETER :: LAMMAX = 0.49999

Referenced by insnl3().

w3gdatmd::nk
integer, pointer nk
Definition: w3gdatmd.F90:1230
w3gdatmd::snlcs
real, dimension(:), pointer snlcs
Definition: w3gdatmd.F90:1361
w3gdatmd::dth
real, pointer dth
Definition: w3gdatmd.F90:1232
w3gdatmd::frq
real, dimension(:), pointer frq
Definition: w3gdatmd.F90:1361
w3gdatmd::snlmsc
real, pointer snlmsc
Definition: w3gdatmd.F90:1360
w3gdatmd::nfrmax
integer, pointer nfrmax
Definition: w3gdatmd.F90:1356
constants::dera
real, parameter dera
DERA Conversion factor from degrees to radians.
Definition: constants.F90:77
w3dispmd::wavnu3
pure subroutine wavnu3(SI, H, K, CG)
Definition: w3dispmd.F90:347
w3gdatmd::qst3
real, dimension(:,:,:), pointer qst3
Definition: w3gdatmd.F90:1361
w3gdatmd::sig
real, dimension(:), pointer sig
Definition: w3gdatmd.F90:1234
w3gdatmd::qst2
real, dimension(:,:,:), pointer qst2
Definition: w3gdatmd.F90:1361
w3gdatmd::nspmax
integer, pointer nspmax
Definition: w3gdatmd.F90:1356
w3gdatmd::fachfe
real, pointer fachfe
Definition: w3gdatmd.F90:1232
constants::rade
real, parameter rade
RADE Conversion factor from radians to degrees.
Definition: constants.F90:76
w3gdatmd::snll
real, dimension(:), pointer snll
Definition: w3gdatmd.F90:1361
w3gdatmd::qst6
real, dimension(:,:,:), pointer qst6
Definition: w3gdatmd.F90:1361
w3gdatmd::snlsfd
real, pointer snlsfd
Definition: w3gdatmd.F90:1360
expand
subroutine expand(SPEC)
Expand spectrum, subroutine used to simplify addressing.
Definition: w3snl3md.F90:644
w3gdatmd::nthexp
integer, pointer nthexp
Definition: w3gdatmd.F90:1356
w3gdatmd::nspmx2
integer, pointer nspmx2
Definition: w3gdatmd.F90:1356
w3odatmd::ndse
integer, pointer ndse
Definition: w3odatmd.F90:456
w3gdatmd::qst1
integer, dimension(:,:,:), pointer qst1
Definition: w3gdatmd.F90:1359
w3gdatmd::nkd
integer, parameter nkd
Definition: w3gdatmd.F90:636
w3gdatmd::xsi
real, dimension(:), pointer xsi
Definition: w3gdatmd.F90:1361
w3gdatmd::snlnsc
real, pointer snlnsc
Definition: w3gdatmd.F90:1360
w3dispmd::wavnu2
subroutine wavnu2(W, H, K, CG, EPS, NMAX, ICON)
Definition: w3dispmd.F90:204
w3servmd
Definition: w3servmd.F90:3
w3gdatmd::nthmax
integer, pointer nthmax
Definition: w3gdatmd.F90:1356
maxlam
real function maxlam(MU, THETA)
Calculate maximum allowed lambda for quadruplet configuration.
Definition: w3snl3md.F90:1602
constants::tpiinv
real, parameter tpiinv
TPIINV Inverse of 2*Pi.
Definition: constants.F90:74
w3gdatmd::nth
integer, pointer nth
Definition: w3gdatmd.F90:1230
w3gdatmd::qst5
real, dimension(:,:,:), pointer qst5
Definition: w3gdatmd.F90:1361
w3odatmd
Definition: w3odatmd.F90:3
w3gdatmd::igrid
integer igrid
Definition: w3gdatmd.F90:618
w3gdatmd::snlsfs
real, pointer snlsfs
Definition: w3gdatmd.F90:1360
w3gdatmd::snlnq
integer, pointer snlnq
Definition: w3gdatmd.F90:1356
minlam
real function minlam(MU, THETA)
Calculate minimum allowed lambda for quadruplet configuration.
Definition: w3snl3md.F90:1534
constants::tpi
real, parameter tpi
TPI 2*Pi.
Definition: constants.F90:72
w3servmd::strace
subroutine strace(IENT, SNAME)
Definition: w3servmd.F90:148
w3gdatmd::facti2
real, pointer facti2
Definition: w3gdatmd.F90:1232
w3adatmd::nfr
integer, pointer nfr
Definition: w3adatmd.F90:657
w3gdatmd::xfr
real, pointer xfr
Definition: w3gdatmd.F90:1232
w3odatmd::ndst
integer, pointer ndst
Definition: w3odatmd.F90:456
constants
Define some much-used constants for global use (all defined as PARAMETER).
Definition: constants.F90:20
w3gdatmd
Definition: w3gdatmd.F90:16
w3dispmd::wavnu1
subroutine wavnu1(SI, H, K, CG)
Definition: w3dispmd.F90:85
w3servmd::extcde
subroutine extcde(IEXIT, UNIT, MSG, FILE, LINE, COMM)
Definition: w3servmd.F90:736
w3gdatmd::snlt
real, dimension(:), pointer snlt
Definition: w3gdatmd.F90:1361
w3gdatmd::nfrcut
integer, pointer nfrcut
Definition: w3gdatmd.F90:1356
expnd2
subroutine expnd2(ARIN, AROUT)
Expand spectrum to simplify indirect addressing.
Definition: w3snl3md.F90:717
w3gdatmd::qst4
integer, dimension(:,:,:), pointer qst4
Definition: w3gdatmd.F90:1359
w3gdatmd::nqa
integer, pointer nqa
Definition: w3gdatmd.F90:1356
w3gdatmd::facti1
real, pointer facti1
Definition: w3gdatmd.F90:1232
w3gdatmd::mpars
type(mpar), dimension(:), allocatable, target mpars
Definition: w3gdatmd.F90:1090
w3dispmd
Definition: w3dispmd.F90:3
w3gdatmd::snlm
real, dimension(:), pointer snlm
Definition: w3gdatmd.F90:1361
constants::grav
real, parameter grav
GRAV Acc.
Definition: constants.F90:61
w3gdatmd::nspmin
integer, pointer nspmin
Definition: w3gdatmd.F90:1356
w3gdatmd::nfrmin
integer, pointer nfrmin
Definition: w3gdatmd.F90:1356
w3gdatmd::snlcd
real, dimension(:), pointer snlcd
Definition: w3gdatmd.F90:1361