146 subroutine z_fileio(filename,qual,iufind,iunit,iostat)
234 character(len=*),
intent(inout) :: filename
235 character(len=2),
intent(in) :: qual
236 integer,
intent(in) :: iufind
237 integer,
intent(inout) :: iunit
238 integer,
intent(out) :: iostat
278 character(len=7) :: cstat
279 character(len=11) :: cform
289 if(iufind==1) iunit = -1
293 if(iufind/=0 .and. iufind/=1)
then 295 'Z_FILEIO: Incorrect value for IUFIND:',iufind
300 if (
parll .and. iunit <= 0 )
then 301 ilpos = index( filename,
' ' )-1
302 write(filename(ilpos+1:ilpos+4),33)
inode 308 if(
i_print>=1)
write(
i_out,*)
'Z_FILEIO/A:',trim(filename),
' ', &
311 if (index(
'ORSUD',qual(1:1)) ==0 .or. &
312 index(
'FUB',qual(2:2)) ==0)
then 313 if(
i_print > 0)
write(
i_out,*)
'Incorrect file qualifier' 316 if(qual(1:1) ==
'O') cstat =
'old' 317 if(qual(1:1) ==
'R') cstat =
'replace' 318 if(qual(1:1) ==
'S') cstat =
'scratch' 319 if(qual(1:1) ==
'U') cstat =
'unknown' 320 if(qual(1:1) ==
'D') cstat =
'delete' 322 if(qual(2:2) ==
'F') cform =
'formatted' 323 if(qual(2:2) ==
'U') cform =
'unformatted' 324 if(qual(2:2) ==
'B') cform =
'binary' 325 if(qual(2:2) ==
'R') cform =
'unformatted' 329 inquire(file=filename,exist=lexist)
330 if(
i_print >=2)
write(
i_out,*)
'Z_FILEIO file exists?:', &
331 trim(filename),
':',lexist
335 if(lexist .and. qual(1:1)==
'D')
then 336 inquire(file=filename,opened=lopen)
338 inquire(file=filename,number=junit)
340 if(iufind == 1)
call z_flunit(iunit,iuerr)
343 open(file=filename,unit=junit,form=cform,iostat=iostat)
350 close(junit,status=cstat)
357 if(
i_print >=2)
write(
i_out,*)
'Z_FILEIO: File exists:', &
359 inquire(file=filename,opened=lopen)
361 if(
i_print >=2)
write(
i_out,*)
'Z_FILEIO: File is opened:', &
367 inquire(file=filename,number=junit)
369 'Z_FILEIO: File is connected to unit:', junit
376 'Z_FILEIO: File is not connected to a unit number' 379 'Z_FILEIO: Assign user defined unit number:',iunit
380 elseif(iufind==1)
then 383 'Z_FILEIO: New unit number IUNIT:',iunit
388 open(file=filename,unit=junit,form=cform,status=cstat)
399 write(
i_out,*)
'Z_FILEIO: File does not exist !' 400 write(
i_out,*)
'Z_FILEIO: Qual:',qual(1:1)
403 if(index(
'SRU',qual(1:1)) > 0)
then 407 'Z_FILEIO: New unit number IUNIT:',iunit
414 open(file=filename,unit=junit,form=cform,iostat=iuerr)
428 elseif(
'O'==qual(1:1))
then 430 'Z_FILEIO: File cannot be opened because it does not exist' 438 if(
i_print>=1)
write(
i_out,*)
'Z_FILEIO/Z:',trim(filename),
' ', &
492 integer,
intent(inout) :: iunit
566 integer,
intent(out) :: iunit
567 integer,
intent(out) :: ierr
599 integer,
parameter :: lu_min=60
600 integer,
parameter :: lu_max=200
604 integer,
parameter :: lu_nr=3
605 integer lu_not(lu_nr)
615 write(
i_out,*)
'Z_FLUNIT: forbidden :',lu_not
616 write(
i_out,*)
'Z_FLUNIT: lu_min lu_max :',lu_min,lu_max
621 if(lu_min >= lu_max)
then 624 'Z_FLUNIT: Incorrect boundaries for LU_MIN & LU_MAX:', &
632 do while (iunit ==-1)
636 inquire(unit=junit,opened=lopen)
642 if(lu_not(i_not)==junit)
then 645 'Z_FLUNIT: a forbidden unit number was encountered:',junit
649 if(lopen.or.lnot)
then 654 if(junit > lu_max)
exit 659 'ERROR in Z_FLUNIT: No free unit number could be found' 671 INTEGER,
intent(in) :: n
672 real,
intent(in) :: x1
673 real,
intent(in) :: x2
674 real,
intent(out) :: x(n)
675 real,
intent(out) :: w(n)
678 parameter(eps=3.d-14)
680 DOUBLE PRECISION p1,p2,p3,pp,xl,xm,z,z1
686 z=cos(3.141592654d0*(i-.25d0)/(n+.5d0))
693 p1=((2.d0*j-1.d0)*z*p2-(j-1.d0)*p3)/j
695 pp=n*(z*p1-p2)/(z*z-1.d0)
698 if(
abs(z-z1).gt.eps)
goto 1
701 w(i)=2.d0*xl/((1.d0-z*z)*pp*pp)
709 subroutine z_cmpcg(sigma,depth,grav_w,cg)
759 real,
intent(in) :: sigma
760 real,
intent(in) :: depth
761 real,
intent(in) :: grav_w
762 real,
intent(out) :: cg
769 k =
z_wnumb(sigma,depth,grav_w)
771 if(depth <= 0. .or. sigma <= 0.)
then 774 if(depth*k > 30.)
then 775 cg = grav_w/(2.*sigma)
777 cg = sigma/k*(0.5+depth*k/sinh(2.*depth*k))
785 subroutine z_intp1(x1,y1,x2,y2,n1,n2,ierr)
843 integer,
intent(in) :: n1
844 integer,
intent(in) :: n2
845 real,
intent(in) :: x1(n1)
846 real,
intent(in) :: y1(n1)
847 real,
intent(in) :: x2(n2)
848 real,
intent(out) :: y2(n2)
849 integer,
intent(out) :: ierr
890 real,
parameter :: eps=1.e-20
910 if (
abs(xmin1-xmax1) < eps .or.
abs(x1(1)-x1(n1)) < eps)
then 915 if ((
abs(xmin2-xmax2) < eps .or.
abs(x2(1)-x2(n2)) < eps) .and. &
923 if(x1(1) < x1(n1))
then 925 if(x1(i1) > x1(i1+1))
then 927 write(*,*)
'z_intp1: i1 x1(i1) x1(i1+1):',i1,x1(i1),x1(i1+1)
933 if(x2(i2) > x2(i2+1))
then 935 write(*,*)
'z_intp1: i2 x2(i2) x2(i2+1):',i2,x2(i2),x2(i2+1)
942 if(x1(i1) < x1(i1+1))
then 944 write(*,*)
'z_intp1: i1 x1(i1) x1(i1+1):',i1,x1(i1),x1(i1+1)
950 if(x2(i2) < x2(i2+1))
then 952 write(*,*)
'z_intp1: i2 x2(i2) x2(i2+1):',i2,x2(i2),x2(i2+1)
967 do while (s1 <= s2 .and. i1 < n1)
979 fac = ds/(x1(i1)-x1(i1-1))
980 y2(i2) = y1(i1-1) + fac*(y1(i1)-y1(i1-1))
985 if(i2==n2 .and. s2>s1) y2(n2) = y1(n1)
1043 integer,
intent(in) :: npol
1044 real,
intent(in) :: xpol(npol)
1045 real,
intent(in) :: ypol(npol)
1046 real,
intent(out) :: area
1055 real xmin,xmax,ymin,ymax
1078 xmean = 0.5*(xmin + xmax)
1079 ymean = 0.5*(ymin + ymax)
1090 if(ipol==npol) ipol1 = 1
1096 darea = 0.5*((xa-xmean)*(yb-ymean) - (xb-xmean)*(ya-ymean))
1098 sumx = sumx + darea*(xa+xb+xmean)/3.
1099 sumy = sumy + darea*(ya+yb+ymean)/3.
1148 integer,
intent(in) :: nx
1149 real,
intent(in) :: x(nx)
1150 real,
intent(out) :: dx(nx)
1161 dx(ix) = 0.5 * (x(ix+1) - x(ix-1))
1165 dx(1) = dx(2)*dx(2)/dx(3)
1166 dx(nx) = dx(nx-1)*dx(nx-1)/dx(nx-2)
1176 real function z_root2(func,x1,x2,xacc,iprint,ierr)
1265 real,
intent (in) :: x1
1266 real,
intent (in) :: x2
1267 real,
intent (in) :: xacc
1268 integer,
intent (in) :: iprint
1269 integer,
intent (out) :: ierr
1278 parameter(maxit = 20)
1300 if(luprint > 0)
then 1301 inquire(unit=luprint,opened=lopen)
1304 write(*,
'(a,i4)')
'Z_ROOT2: invalid unit number:',iprint
1322 if((fl > 0. .and. fh < 0.) .or. (fl < 0. .and. fh > 0.))
then 1330 s = sqrt(fm**2-fl*fh)
1331 if(s == 0.)
goto 9000
1332 xnew = xm+(xm-xl)*(sign(1.,fl-fh)*fm/s)
1337 if (
abs(xnew-zriddr) <= xacc)
then 1344 if (fnew == 0.)
goto 9000
1346 if(sign(fm,fnew) /= fm)
then 1351 elseif(sign(fl,fnew) /= fl)
then 1354 elseif(sign(fh,fnew) /= fh)
then 1362 if(
abs(xh-xl) <= xacc)
goto 9000
1364 if(luprint > 0)
write(luprint,
'(a,i4,5e14.6)') &
1365 'Z_ROOT2: iter,x1,x2,|x1-x2|,xacc,z:', iter,xl,xh, &
1366 abs(xl-xh),xacc,fnew
1370 if(luprint > 0)
write(luprint,
'(a)')
'Z_ROOT2: -> ierr=2' 1372 else if (fl == 0.)
then 1374 else if (fh == 0.)
then 1386 if(luprint > 0)
write(luprint,
'(a,2i3,5e13.5)') &
1387 'Z_ROOT2: ierr,iter,xl,xh,acc,x0,z0:', ierr,iter,xl,xh,xacc, &
1388 z_root2,func(z_root2)
1439 character(len=*),
intent(inout) :: str
1448 integer i,ial,iau,izl
1457 if(ichar(str(i:i)) >= ial.and. ichar(str(i:i)) <= izl)
then 1458 str(i:i) = char(ichar(str(i:i))-ial+iau)
1466 real function z_wnumb(w,d,grav_w)
1516 real,
intent(in) :: w
1517 real,
intent(in) :: d
1518 real,
intent(in) :: grav_w
1530 if(d<=0 .or. w<= 0.)
then 1535 xx = y*(y+1./(1.+y*(0.66667+y*(0.35550+y*(0.16084+y* &
1536 (0.06320+y*(0.02174+y*(0.00654+y*(0.00171+y* &
1537 (0.00039+y*0.00011))))))))))
1910 real,
allocatable ::
a(:,:)
1991 subroutine xnl_init(sigma,dird,nsigma,ndir,pftail,x_grav,depth, &
1992 ! ndepth,jquad,iqgrid,iproc,ierror)
1993 ndepth,jquad,iqgrid,ierror)
2072 integer,
intent(in) :: nsigma
2073 integer,
intent(in) :: ndir
2074 integer,
intent(in) :: ndepth
2075 real,
intent(in) :: sigma(nsigma)
2076 real,
intent(in) :: dird(ndir)
2077 real,
intent(in) :: pftail
2078 real,
intent(in) :: depth(ndepth)
2079 real,
intent(in) :: x_grav
2080 integer,
intent(in) :: jquad
2081 integer,
intent(in) :: iqgrid
2083 integer,
intent(out) :: ierror
2194 dstep = dird(2)-dird(1)
2195 dgap = 180.- abs(180.- abs(dird(1)-dird(ndir)))
2202 if(abs(dstep-dgap) < 0.001)
then 2204 write(
iscreen,
'(a,i4,2f10.2)')
'XNL_INIT iq_grid dstep dgap:', &
2211 if(abs(dird(1)+dird(ndir)) > 0.01)
then 2212 write(
iscreen,
'(a,i4,2f10.2)') &
2213 'XNL_INIT iq_grid dird(1) dird(n):',
iq_grid,dird(1),dird(ndir)
2248 call q_error(
'e',
'FILEIO',
'Problem in deleting error file *.ERR')
2269 write(
luq_log,
'(2a,i4)')
'XNL_INIT: ', &
2271 write(
luq_log,
'(2a,i4)')
'XNL_INIT: ', &
2273 write(
luq_log,
'(2a,i4)')
'XNL_INIT: ', &
2276 write(
luq_log,
'(2a,i4)')
'XNL_INIT: ', &
2278 write(
luq_log,
'(2a,i4)')
'XNL_INIT: ', &
2280 write(
luq_log,
'(2a,i4)')
'XNL_INIT: ', &
2284 '---------------------------------------------------------------' 2287 'Solution of Boltzmann integral using Webb/Resio/Tracy method' 2289 '---------------------------------------------------------------' 2291 write(
luq_prt,
'(a)')
'Initialisation' 2335 call q_error(
'w',
'DEPTH',
'Invalid depth')
2336 write(
luq_err,
'(a,e12.5,f10.2)')
'Incorrect depth & minimum:', &
2347 if(jquad==1 .and. ndepth > 0)
then 2349 'XNL_INIT: For deep water only one grid suffices' 2383 subroutine xnl_main(aspec,sigma,angle,nsig,ndir,depth,iquad,xnl, &
2457 integer,
intent(in) :: nsig
2458 integer,
intent(in) :: ndir
2459 integer,
intent(in) :: iquad
2460 integer,
intent(in) :: iproc
2462 real,
intent(in) :: aspec(nsig,ndir)
2463 real,
intent(in) :: sigma(nsig)
2464 real,
intent(in) :: angle(ndir)
2465 real,
intent(in) :: depth
2466 real,
intent(out) :: xnl(nsig,ndir)
2468 real,
intent(out) :: diag(nsig,ndir)
2469 integer,
intent(out) :: ierror
2496 integer,
save :: i_qmain
2517 i_qmain = i_qmain + 1
2521 write(
luq_prt,
'(a,i4,f16.3,i4)') &
2522 'XNL_MAIN: Input arguments: iquad depth iproc:', &
2537 call q_error(
'w',
'DEPTH',
'Zero transfer returned')
2552 if(iquad>=1 .and. iquad <=3)
then 2555 call q_xnl4v4(aspec,sigma,angle,nsig,ndir,depth,xnl,diag,ierror)
2558 call q_error(
'e',
'wrtvv',
'Problem in Q_XNL4V4')
2577 'XNL_MAIN depth scale factor:',q_dfac
2583 call q_chkcons(xnl,nsig,ndir,sum_e,sum_a,sum_mx,sum_my)
2586 write(
luq_prt,
'(a)')
'XNL_MAIN: Conservation checks' 2587 write(
luq_prt,
'(a,4e13.5)')
'XNL_MAIN: E/A/MOMX/MOMY:', &
2588 sum_e,sum_a,sum_mx,sum_my
2598 write(
luq_log,
'(a,i4)')
'XNL_MAIN: Number of errors :',
iq_err 2608 subroutine q_addtail(xnl,diag,nsig,na,pf_tail)
2659 integer,
intent(in) :: nsig
2660 integer,
intent(in) :: na
2661 real,
intent(inout) :: xnl(nsig,na)
2662 real,
intent(inout) :: diag(nsig,na)
2663 real,
intent(in) :: pf_tail
2689 real diag_tail(nsig)
2701 if(
q_sig(isig)<sig_tail) isig_tail = isig
2709 x_tail = 11+3*pf_tail-1
2710 d_tail = 12+2*pf_tail-1
2713 write(
iscreen,
'(a,2i4,2f8.2)') &
2714 'Q_ADDTAIL nsig,isig_tail,qf_tail:',nsig,isig_tail,pf_tail
2715 write(
iscreen,
'(a,3f8.2)')
'Q_ADDTAIL qf_tail, x_tail d_tail:', &
2716 pf_tail,x_tail,d_tail
2719 do isig=isig_tail,nsig
2720 xnl_tail(isig) = (
q_sig(isig)/
q_sig(isig_tail))**x_tail
2721 diag_tail(isig) = (
q_sig(isig)/
q_sig(isig_tail))**d_tail
2726 do isig=isig_tail,nsig
2728 xnl(isig,ia) = xnl(isig_tail,ia) *xnl_tail(isig)
2729 diag(isig,ia) = diag(isig_tail,ia)*diag_tail(isig)
2839 'Q_ALLOCATE: mkq maq mlocus klocus:',mkq,maq,
mlocus,
klocus 2841 if (
allocated (
q_xk))
deallocate (
q_xk)
2843 if (
allocated (
q_sk))
deallocate (
q_sk)
3008 if (
allocated(
r_zz))
deallocate (
r_zz)
3024 if (
allocated(
r_ws))
deallocate (
r_ws)
3053 if (
allocated(
t_zz))
deallocate (
t_zz)
3059 if (
allocated(
t_ws))
deallocate (
t_ws)
3064 if (
allocated(
dt13))
deallocate (
dt13)
3069 if (
allocated(
q_k))
deallocate (
q_k)
3071 if (
allocated(
q_dk))
deallocate (
q_dk)
3075 if (
allocated(
q_f))
deallocate (
q_f)
3077 if (
allocated(
q_df))
deallocate (
q_df)
3083 if (
allocated(
q_cg))
deallocate (
q_cg)
3085 if (
allocated(
q_a))
deallocate (
q_a)
3087 if (
allocated(
q_ad))
deallocate (
q_ad)
3093 if (
allocated(
a))
deallocate (
a)
3095 if (
allocated(
nk1d))
deallocate (
nk1d)
3101 write(
luq_log,
'(a)')
'Q_ALLOCATE: size of arrays' 3102 write(
luq_log,
'(a,i4)')
'Q_ALLOCATE: mkq :',mkq
3103 write(
luq_log,
'(a,i4)')
'Q_ALLOCATE: maq :',maq
3104 write(
luq_log,
'(a,i4)')
'Q_ALLOCATE: nkq :',
nkq 3105 write(
luq_log,
'(a,i4)')
'Q_ALLOCATE: naq :',
naq 3202 call q_error(
'e',
'CONFIG',
'Incorrect power of spectral: qf_tail')
3206 'Invalid option for coupling coefficient iq_cple')
3209 call q_error(
'e',
'CONFIG',
'iq_compact /= 0,1')
3212 call q_error(
'e',
'CONFIG',
'iq_filt /= 0,1')
3215 call q_error(
'e',
'CONFIG',
'iq_gauleg <0')
3218 call q_error(
'e',
'CONFIG',
'iq_geom /= 0,1')
3221 call q_error(
'e',
'CONFIG',
'iq_interp /= 1,2')
3225 'Invalid combination of iq_disp & iq_geom')
3226 write(
luq_err,
'(a,2i4)')
'Q_CHKCONFIG: iq_disp iq_geom:', &
3231 call q_error(
'e',
'CONFIG',
'Invalid value for IQ_DSCALE')
3237 'Lumping and Gauss-Legendre interpolation not together')
3238 write(
luq_err,
'(a,2i4)')
'Q_CHKCONFIG: iq_lump iq_gauleg:', &
3244 'Incorrect value for IQ_DSCALE, (0,1)')
3248 'Incorrect value for IQ_DISP [DISP],(1,2) ')
3252 'CONFIG',
'Incorrect value for IQ_GRID, (1,2,3)')
3255 call q_error(
'e',
'CONFIG',
'Invalid value for iq_integ')
3260 call q_error(
'e',
'CONFIG',
'Incorrect value for IQ_LOG, (>=0) ')
3264 'CONFIG',
'Incorrect specifier for locus method')
3267 call q_error(
'e',
'CONFIG',
'Invalid value for iq_lump')
3272 call q_error(
'e',
'CONFIG',
'Invalid value for iq_make')
3278 'CONFIG',
'Incorrect value for IQ_MOD [MOD] (0,1)')
3281 call q_error(
'e',
'CONFIG',
'klocus < nlocus0')
3283 'Q_CHKCONFIG: Lumping or Gauss-Integration enabled when IMOD=0' 3287 call q_error(
'e',
'CONFIG',
'Incorrect value for IQ_PRT, (>=0) ')
3293 call q_error(
'e',
'CONFIG',
'Incorrect value of IQ_SYM /=[0,1]')
3296 call q_error(
'e',
'CONFIG',
'Incorrect value for IQ_TEST, (>=0) ')
3299 call q_error(
'e',
'CONFIG',
'Incorrect value for IQ_TRF ')
3305 call q_error(
'e',
'CONFIG',
'Incorrect value for FQMIN')
3308 call q_error(
'e',
'CONFIG',
'Incorrect value for FQMAX')
3311 call q_error(
'e',
'CONFIG',
'fmax <= fmin')
3314 call q_error(
'e',
'CONFIG',
'Number of wave numbers NKQ < 0')
3317 call q_error(
'e',
'CONFIG',
'Number of directions NKQ < 0')
3321 'CONFIG',
'Preferred number of points on locus NLOCUS0 < 6')
3325 'Sector too small (<40) or too large (>180)')
3329 'Value of Q_DSTEP outside range 0.1 - 1000')
3338 subroutine q_chkcons(xnl,nk,ndir,sum_e,sum_a,sum_mx,sum_my)
3395 integer,
intent(in) :: nk
3396 integer,
intent(in) :: ndir
3397 real,
intent(in) :: xnl(nk,ndir)
3398 real,
intent(out) :: sum_e
3399 real,
intent(out) :: sum_a
3400 real,
intent(out) :: sum_mx
3401 real,
intent(out) :: sum_my
3451 momx = aa*kk*cos(
q_a(ia))
3452 momy = aa*kk*sin(
q_a(ia))
3454 sum_a = sum_a + aa*qq
3455 sum_e = sum_e + ee*qq
3456 sum_mx = sum_mx + momx*qq
3457 sum_my = sum_my + momy*qq
3466 subroutine q_chkres(k1x,k1y,k2x,k2y,k3x,k3y,k4x,k4y,dep, &
3467 sum_kx,sum_ky,sum_w)
3524 real,
intent(in) :: k1x
3525 real,
intent(in) :: k1y
3526 real,
intent(in) :: k2x
3527 real,
intent(in) :: k2y
3528 real,
intent(in) :: k3x
3529 real,
intent(in) :: k3y
3530 real,
intent(in) :: k4x
3531 real,
intent(in) :: k4y
3532 real,
intent(in) :: dep
3533 real,
intent(out) :: sum_kx
3534 real,
intent(out) :: sum_ky
3535 real,
intent(out) :: sum_w
3559 sum_kx = (k1x + k2x) - (k3x + k4x)
3560 sum_ky = (k1y + k2y) - (k3y + k4y)
3564 w1 =
x_disper(sqrt(k1x**2 + k1y**2),dep)
3565 w2 =
x_disper(sqrt(k2x**2 + k2y**2),dep)
3566 w3 =
x_disper(sqrt(k3x**2 + k3y**2),dep)
3567 w4 =
x_disper(sqrt(k4x**2 + k4y**2),dep)
3569 sum_w = w1 + w2 - (w3 + w4)
3637 real,
intent(out) :: ka,kb
3638 real,
intent(out) :: km
3639 real,
intent(out) :: kw
3640 real,
intent(out) :: loclen
3700 eps = 10.*epsilon(1.)
3711 k1m = sqrt(
k1x**2 +
k1y**2)
3712 k3m = sqrt(
k3x**2 +
k3y**2)
3731 write(
luq_tst,
'(a,4f11.5)')
'Q_CMPLOCUS: k1x/y k3x/y :', &
3733 write(
luq_tst,
'(a,3f11.5)')
'Q_CMPLOCUS: Px Py Pmag :', &
3735 write(
luq_tst,
'(a,3f11.4)')
'Q_CMPLOCUS: Pang Sang Xang :', &
3737 write(
luq_tst,
'(a,2f11.5)')
'Q_CMPLOCUS: k1m,k3m :', &
3739 write(
luq_tst,
'(a,f11.5)')
'Q_CMPLOCUS: q :',
q 3780 call q_polar2(ka,kb,kx_beg,ky_beg,kx_end,ky_end,loclen,ierr)
3785 area2 = 2.*
pih*(kb-ka)*0.5*kw
3786 ratio = max(area1/area2,area2/area1)
3788 if(
iq_test>=2 .and.ratio>1.5)
then 3789 write(
luq_tst,
'(a,4i4,2e12.4,4f12.5)') &
3801 if(ratio>1.5 .and. k3m/k1m < 100.)
then 3802 call q_error(
'e',
'LOCUS',
'Severe problem in POLAR2')
3803 write(
luq_err,
'(a)')
'Q_CMPLOCUS: ratio > 1.5' 3856 ds_loc(iloc) = 0.5*(dsp + dsm)
3879 write(
luq_tst,
'(a,i4,4f10.5,2e12.5)') &
3880 'Q_CMPLOCUS: k2x k2y ds s jac cp:', &
3885 if(itest >= 1)
write(
luq_tst,
'(a,2f10.5)') &
3886 'Q_CMPLOCUS: length of locus:',sum, &
3975 integer,
intent(in) :: itask
3979 integer,
intent(out) :: igrid
4042 integer iaz,ikz,jaz,jkz
4043 integer iz_geom,iz_disp,iz_cple
4052 real,
allocatable :: z_ad(:),z_sig(:)
4061 'Q_CTRGRID: input arguments: itask depth:',itask,
q_depth 4075 write(
luq_prt,
'(a,3i5)')
'Q_CTRGRID: naq nkq nlocus0:',
naq,
nkq, &
4077 write(
luq_prt,
'(a,3i3)')
'Q_CTRGRID: iq_grid,iq_geom,iq_disp:', &
4100 jdep = min(99999,jdep)
4105 'Q_CTRGRID: q_depth q_dstep s_depth idep jdep:', &
4112 call q_error(
'e',
'DISPER',
'Incorrect value for IQ_DISP')
4123 'Q_CTRGRID: Last and actual BQF file: ', &
4128 'Q_CTRGRID: Rereading of bqfile skipped: ',
lastquadfile 4130 'Q_CTRGRID: Rereading of bqfile skipped: ',
lastquadfile 4136 write(
luq_prt,
'(2a)')
'Q_CTRGRID: Header line of grid file:', &
4138 write(
luq_prt,
'(2a)')
'Q_CTRGRID: Name of BINARY grid file:', &
4160 if(
luq_bqf > 0 .and. iuerr ==0)
then 4163 'Q_CTRGRID: Binary grid file detected: ',trim(
bqname)
4164 write(
luq_prt,
'(a,i4)')
'Q_CTRGRID: Connected to unit:', &
4176 call q_error(
'w',
'READBQF',
'Read error for header in BQF file')
4177 write(
luq_err,
'(a)')
'BQF file deleted' 4206 'Q_CTRGRID: Header in binary quad file :', &
4209 'Q_CTRGRID: Expected header of binary quad file:', &
4211 write(
luq_prt,
'(a)')
'Q_CTRGRID: The file headers disagree' 4213 'Q_CTRGRID: A new grid will be generated' 4221 allocate (z_sig(nkz),z_ad(naz))
4224 read(
luq_bqf) iz_geom,iz_disp,iz_cple
4228 write(
luq_prt,
'(a)')
'Q_CTRGRID: Contents of BQF file' 4230 write(
luq_prt,
'(a,i4)')
'Q_CTRGRID: NK:',nkz
4231 write(
luq_prt,
'(a,i4)')
'Q_CTRGRID: NA:',naz
4239 if(abs(
q_ad(iaz)-z_ad(iaz)) > 0.01)
then 4240 write(
luq_prt,
'(a)')
'Q_CTRGRID: Directions do not agree' 4242 write(
luq_prt,
'(1x,a,i4,2f10.3)')
'iaz q_ad z_ad:',jaz, &
4254 if(abs(
q_sig(ikz)-z_sig(ikz)) > 0.01)
then 4255 write(
luq_prt,
'(a)')
'Q_CTRGRID: Wave numbers do not agree' 4257 write(
luq_prt,
'(1x,a,i4,2f10.3)')
'ikz q_k z_sig:',jkz, &
4258 q_sig(jkz),z_sig(jkz)
4269 if(abs(z_depth-s_depth) > 0.09 .and.
iq_disp > 1 .and. &
4271 write(
luq_prt,
'(a)')
'Q_CTRGRID: Water depths do not agree' 4272 write(
luq_prt,
'(a,2f16.2)')
'Q_CTRGRID: q_depth z_depth:', &
4281 'Q_CTRGRID: Existing BQF-file invalid, it will be closed' 4307 write(
luq_log,
'(a)')
'Q_CTRGRID: New grid will be generated' 4308 write(
luq_log,
'(a,a)')
'Q_CTRGRID: Name of BQF file:', &
4315 'Q_CTRGRID: Generating wave number grid for '// &
4316 'quadruplet interactions: ',trim(
bqname)
4337 'Q_CTRGRID: Grid generation completed succesfully' 4346 'Q_CTRGRID: Reading existing grid: ',trim(
bqname)
4348 'Q_CTRGRID: Existing grid will be read:',trim(
bqname)
4350 'Q_CTRGRID: Existing grid will be read:',trim(
bqname)
4373 'Read error for test data in BQF file')
4375 'BQF file probably generated with test option off' 4396 'Q_CTRGRID: on exit igrid & q_depth:',igrid,
q_depth 4398 if (
allocated(z_ad))
deallocate(z_ad,z_sig)
4401 'Q_CTRGRID: On exit LASTQUADFILE & q_depth: ',
lastquadfile, &
4404 'Q_CTRGRID: On exit LASTQUADFILE & q_depth: ',
lastquadfile, &
4412 subroutine q_dscale(n,sigma,angle,nsig,nang,depth,grav_w,q_dfac)
4472 integer,
intent (in) :: nsig
4473 integer,
intent (in) :: nang
4474 real,
intent(in) :: n(nsig,nang)
4475 real,
intent(in) :: sigma(nsig)
4476 real,
intent(in) :: angle(nang)
4477 real,
intent(in) :: depth
4478 real,
intent(in) :: grav_w
4479 real,
intent(out) :: q_dfac
4526 call z_steps(sigma,dsigma,nsig)
4527 delta = angle(2)-angle(1)
4539 dnn = n(isig,iang)*dsigma(isig)*delta
4541 sumk = sumk + 1./sqkk*dnn
4549 kms = (sum0/sumk)**2
4550 kd = max(0.5,0.75*kms*depth)
4551 q_dfac = 1+5.5/kd*(1.-5./6.*kd)*exp(-5./4.*kd)
4567 subroutine q_error(err_type,err_name,err_msg)
4622 character(len=1),
intent(in) :: err_type
4625 character(len=*),
intent(in) :: err_name
4626 character(len=*),
intent(in) :: err_msg
4645 character(len=80) qline
4659 'Q_ERROR: '//trim(
qbase)//
'.ERR connected to unit:',
luq_err 4666 '------------------------------------------------------------' 4671 if(index(
'wW',err_type) > 0)
then 4673 write(
luq_err,
'(a,i4)')
'Warning or non-terminating error:', &
4675 write(
luq_err,
'(a,a)')
'Name of error:',trim(err_name)
4677 elseif(index(
'eE',err_type) > 0)
then 4680 write(
luq_err,
'(a,a)')
'Name of error:',trim(err_name)
4681 write(*,
'(1x,a,i4)')
'Terminating error:',
iq_err 4682 write(*,
'(1x,a,a)')
'Name of error:',trim(err_name)
4688 ntext = len_trim(err_name)
4697 ' does not exist in current directory' 4701 'Q_ERROR: File Q_ERROR.TXT connected to unit:',
luq_txt 4707 read(
luq_txt,
'(a)',iostat=iend) qline
4712 if(qline(1:ntext) == err_name(1:ntext))
then 4715 'Explanation of error, and recommended action' 4717 '--------------------------------------------' 4718 write(
luq_err,
'(a)') trim(qline)
4723 do while (ispace ==1)
4724 read(
luq_txt,
'(a)',iostat=iend) qline
4729 if(qline(1:1) ==
' ')
then 4730 write(
luq_err,
'(a)') trim(qline)
4746 'Q_ERROR: File ',trim(
qf_error),
' disconnected from unit:', &
4751 if(len_trim(err_msg) > 0)
then 4753 write(
luq_err,
'(a)')
'Additional message from point of occurrence:' 4754 write(
luq_err,
'(a)') trim(err_msg)
4761 write(
luq_err,
'(a)')
'Trace of error' 4762 write(
luq_err,
'(a)')
'--------------' 4769 if(
iq_warn > 10) stop
'Too many warnings' 4856 integer,
intent(in) :: ik1
4857 integer,
intent(in) :: ia1
4858 integer,
intent(in) :: ik3
4859 integer,
intent(in) :: ia3
4860 integer,
intent(out) :: ifnd
4900 real xt2(nlocus),yt2(nlocus)
4901 real xt4(nlocus),yt4(nlocus)
4907 integer ikmin,ja1,ja3,jk1,jk3,itmin
4926 ikmin = min(ik1,ik3)
4927 ikdif = abs(ik1-ik3)
4937 itmin = min(it1,it3)
4938 iadif = abs(it1-it3)
4943 'Q_GETLOCUS: it1,it3,itmin,iadif,ja1,ja3:',it1,it3,itmin,iadif, &
4951 if (iadif > nhalf)
then 4952 if(it1 > nhalf) it1 = it1 -
naq 4953 if(it3 > nhalf) it3 = it3 -
naq 4955 itmin = min(it1,it3)
4956 ibdif = (
naq - abs(
naq-2*abs(it1-it3)))/2
4971 kmem = (jk3-jk1+1) - (jk1-2*
nkq-2)*(jk1-1)/2
4977 'Q_GETLOCUS: ik1 ia1 ik3 ia3 kmem amem:',ik1,ia1,ik3,ia3,kmem,amem
4983 if (amem >
iamax)
then 4985 call q_error(
'e',
'MEMORY',
'Incorrect addres')
4986 write(
luq_err,
'(a,2i4)')
'Q_GETLOCUS: iamax,amem:',
iamax,amem
5035 lambda =
q_kfac**(ikmin-1.)
5038 j_lambda = 1./sqrt(lambda)
5039 c_lambda = lambda**6
5043 zz_lambda = lambda*c_lambda/j_lambda
5052 if(ik3 > ik1 .and. it3 >= it1)
then 5069 elseif(ik3 > ik1 .and. it3 < it1)
then 5087 elseif(ik1 > ik3 .and. it3 >= it1)
then 5104 elseif(ik1 > ik3 .and. it1 > it3)
then 5121 elseif(ik1==ik3 .and. it3 > it1)
then 5138 elseif(ik1==ik3 .and. it1 > it3)
then 5159 t_zz(1:nloc) = lambda*c_lambda/j_lambda *
r_zz(1:nloc)
5174 write(
luq_trf,
'(a)')
'! ik1 ia1 ik3 ia3' 5175 write(
luq_trf,
'(a)')
'! k1x k1y k3x k3y' 5177 '! itrans kmem amem kdif iaref ibeta imirror it1 it3' 5178 write(
luq_trf,
'(a)')
'! lambda depth' 5179 write(
luq_trf,
'(a)')
'! nlocus' 5180 write(
luq_trf,
'(a)')
'! k2x k2y k4x k4y ds jac cple sym zz' 5182 write(
luq_trf,
'(a1,2i3.3,a1,2i3.3,a1)')
'(',ik1,ia1,
'-',ik3, &
5185 q_k(ik1)*sin(
q_a(ia1)), &
5186 q_k(ik3)*cos(
q_a(ia3)), &
5188 write(
luq_trf,
'(9i5)') itrans,kmem,amem,kdif,
iaref,ibeta, &
5211 write(
luq_trf,
'(3i6,4f10.4,e13.5)') iloc,
r_ik2(iloc), &
5220 write(
luq_trf,
'(3i6,4f10.4,e13.5)') iloc,
t_ik2(iloc), &
5235 xt2(iloc) = vk*cos(va)
5236 yt2(iloc) = vk*sin(va)
5242 xt4(iloc) = vk*cos(va)
5243 yt4(iloc) = vk*sin(va)
5244 write(
luq_trf,
'(5f10.4,5e13.5)') xt2(iloc),yt2(iloc), &
5245 xt4(iloc),yt4(iloc), &
5401 write(
luq_prt,
'(a)')
'Basic wave numbers, frequencies' 5413 write(
luq_prt,
'(a,i4,3f10.5,e12.4)') &
5414 'Q_INIT: ikq f sigma k k^p:', &
5423 write(
luq_prt,
'(a)')
'Extended wave numbers and spacing' 5430 elseif(ikq==
nkq)
then 5454 write(
luq_prt,
'(a)')
'Q_INIT: Additional information' 5456 write(
luq_prt,
'(a,i3)')
'Number of frequencies:',
nkq 5457 write(
luq_prt,
'(a,f8.4)')
'Geometric f-spacing factor:',
q_ffac 5458 write(
luq_prt,
'(a,f8.4)')
'Geometric k-spacing factor:',
q_kfac 5463 write(
luq_prt,*)
' i f df sig '// &
5467 write(
luq_prt,
'(1x,i4,7f10.4)') &
5487 'Q_INIT: Index of first direction for reference:',
iaref 5503 'Q_INIT: Range of indices for loop over directions:',
iaq1,
iaq2 5514 'Q_INIT: take care of q_dird1 and check if sector is OK' 5526 write(
luq_prt,
'(a,3f10.3)')
'Q_INIT: d(1),d(n),dsector:', &
5528 write(
luq_prt,
'(a,f6.2,a)')
'Q_INIT: Angular step :', &
5530 write(
luq_prt,
'(a,2f8.2,i4,a)')
'Q_INIT: ang1 ang2 nang :', &
5532 write(
luq_prt,
'(a,i4)')
'Q_INIT: #Angles on circle:', &
5543 write(
luq_prt,
'(a,i4,f10.4,f10.2)')
'Q_INIT: iaq q_a q_ad:', &
5566 write(
luq_tst,
'(a,3i4)')
'Q_INIT: iq_grid iaref iamax:', &
5568 write(
luq_tst,
'(a,4i4)')
'Q_INIT: iaq1 iaq2 iag1 iag2:', &
5573 write(
luq_trf,
'(a)')
'#GRIDINFO#' 5586 real function x_locus1(k2)
5639 real,
intent(in) :: k2
5672 w2 =
sqrtg * sqrt(k2)
5745 real,
intent(in) :: lambda
5777 kk2m = sqrt(kk2x**2 + kk2y**2)
5781 kk4m = sqrt(kk4x**2 + kk4y**2)
5785 w2 =
sqrtg * sqrt(kk2m)
5786 w4 =
sqrtg * sqrt(kk4m)
5806 subroutine q_locpos(ka,kb,km,kw,loclen)
5861 real,
intent (out) :: ka
5862 real,
intent (out) :: kb
5863 real,
intent (out) :: km
5864 real,
intent (out) :: kw
5865 real,
intent (out) :: loclen
5947 kp = sqrt(kpx**2 + kpy**2)
5957 ka = 0.5*(-qs+sqrt(2.0*
pmag-qsq))
5959 kb = (
pmag+qsq)/(2.*qs)
5964 ka = 0.5*(-qs+sqrt(2.0*
pmag-qsq))
5966 kb = (
pmag-qsq)/(2.*qs)
5972 if(itest >= 1)
write(
luq_tst,
'(a,6e12.5)') &
5973 'Q_LOCPOS: q,pmag,ka,kb,za,zb:',qs,
pmag,ka,kb,za,zb
5990 kacc = 10.*max(kk1,kk2)*eps
5994 if(itest>=2)
write(
luq_tst,
'(a,i4)') &
5995 'Q_LOCPOS/Z_ROOT2/IERR/1=',ierr
5998 if(itest >= 1)
write(
luq_tst,
'(a,4f12.5)') &
5999 'Q_LOCPOS: q kk1 kk2 kmin:',
q,kk1,kk2,ka
6011 do while (zz1*zz2 >= 0 .and. iter < maxiter)
6015 if(itest >= 2)
write(
luq_tst,
'(a,i4,3f12.5,2e13.5)') &
6016 'Q_LOCPOS iter q kk1/2 zz1/2:',iter,
q,kk1,kk2,zz1,zz2
6019 if(iter>=maxiter)
then 6020 call q_error(
'e',
'Start kb',
'Too many iterations needed')
6026 kacc = 10.*max(kk1,kk2)*eps
6028 if(ierr>0 .and. itest>=2)
write(
luq_tst,
'(a,i4)') &
6029 'Q_LOCPOS/Z_ROOT2/IERR/2=',ierr
6046 do while (zz1*zz2 >= 0 .and. iter < maxiter)
6050 if(itest >= 2)
write(
luq_tst,
'(a,i4,3f12.5,2e12.5)') &
6051 'Q_LOCPOS: iter q kk1/2 zz1/2:',iter,
q,kk1,kk2,zz1,zz2
6054 if(iter>=maxiter)
then 6055 call q_error(
'e',
'Start ka',
'Too many iterations needed')
6061 kacc = 10.*max(abs(kk1),abs(kk2))*eps
6063 if(ierr > 0 .and. itest>=2)
write(
luq_tst,
'(a,i4)') &
6064 'Q_LOCPOS/Z_ROOT2/IERR/3=',ierr
6076 do while (zz1*zz2 >= 0 .and. iter < maxiter)
6080 if(itest >= 2)
write(
luq_tst,
'(a,i4,3f12.5,2e12.5)') &
6081 'Q_LOCPOS: iter q kk1/2 zz1/2:',iter,
q,kk1,kk2,zz1,zz2
6084 if(iter>=maxiter)
then 6085 call q_error(
'e',
'Start kb',
'Too many iterations needed')
6091 kacc = 10.*max(kk1,kk2)*eps
6093 if(ierr>0.and.itest>=2)
write(
luq_tst,
'(a,i4)') &
6094 'Q_LOCPOS/Z_ROOT2/IERR/4=',ierr
6103 if(itest >= 1)
write(
luq_tst,
'(a,6e12.5)') &
6104 'Q_LOCPOS: q,pmag,ka,kb,za,zb:',
q,
pmag,ka,kb,za,zb
6120 if(itest >= 1)
write(
luq_tst,
'(a,3f12.6)') &
6136 if(itest >= 1)
write(
luq_tst,
'(a,4f10.5,2e12.5)') &
6137 'Q_LOCPOS: k1 k2 z1/2:', &
6138 kk1x,kk1y,kk2x,kk2y,zz1,zz2
6141 do while (zz1*zz2 > 0 .and. iter < maxiter)
6147 if(itest >= 1)
write(
luq_tst,
'(a,i4,4f12.5,2e12.5)') &
6148 'Q_LOCPOS: iter beta1/2 kk2x/y zz1/2:',iter,beta1,beta2, &
6156 write(
luq_tst,
'(a,2f10.4,2e13.5)') &
6157 'Q_LOCPOS: beta1/2 xlocus2(beta1/2):',beta1,beta2, &
6161 bacc = 10.*max(beta1,beta2)*eps
6165 if(itest>=2)
write(
luq_tst,
'(a,i4)') &
6166 'Q_LOCPOS/Z_ROOT2/IERR_W/5=',ierr
6167 call q_error(
'e',
'ROOT2',
'beta')
6176 if(itest >= 1)
write(
luq_tst,
'(a,4f12.6,e12.5)') &
6177 'Q_LOCPOS: betaw kwx kwy kw zw:',betaw,kwx,kwy,kw,zw
6192 a1 = 0.4630151; a2 = 0.1077812;
6193 b1 = 0.2452727; b2 = 0.0412496;
6196 loclen = 4.*max(aa,bb)
6198 loclen = 4.*max(aa,bb)*((1. + a1*mm1 + a2*mm1**2) + &
6199 (b1*mm1 + b2*mm1**2)*log(1/mm1))
6203 write(
luq_tst,
'(a,4f10.5)')
'Q_LOCPOS: aa,bb,mm,mm1:',aa,bb,mm,mm1
6204 write(
luq_tst,
'(a,4f10.5)')
'Q_LOCPOS: length of ellipse:',loclen
6318 integer iaq3,ikq1,ikq3,nkq1
6321 real aa1,aa3,kk1,kk3
6324 integer nztot1,nztot2
6329 real w1k2,w2k2,w3k2,w4k2
6330 real w1k4,w2k4,w3k4,w4k4
6341 real,
allocatable :: xloc(:),yloc(:)
6354 if(
allocated(xloc))
deallocate(xloc) ;
allocate (xloc(
mlocus))
6355 if(
allocated(yloc))
deallocate(yloc) ;
allocate (yloc(
mlocus))
6386 'Q_MAKEGRID: ik1,krefx,krefy:',ikq1,
krefx,
krefy 6394 'Q_MAKEGRID: kk1 kk3 kk3/kk1:',kk1,kk3,kk3/kk1
6398 if(iaq3 ==
iag1 .and. ikq3 == ikq1) cycle
6407 write(
luq_tst,
'(a,2f10.4,2i4)') &
6408 'Q_MAKEGRID: k3 a3 ikq3 iaq3: ',kk3,aa3,ikq3,iaq3
6431 kmem = (ikq3-ikq1+1) - (ikq1-2*
nkq-2)*(ikq1-1)/2;
6455 ik2 = floor(
wk_k2(iloc))
6456 ia2 = floor(
wa_k2(iloc))
6457 wk =
wk_k2(iloc)-real(ik2)
6458 wa =
wa_k2(iloc)-real(ia2)
6459 w1k2 = (1.-wk)*(1.-wa)
6464 ik4 = floor(
wk_k4(iloc))
6465 ia4 = floor(
wa_k4(iloc))
6466 wk =
wk_k4(iloc)-real(ik4)
6467 wa =
wa_k4(iloc)-real(ia4)
6468 w1k4 = (1.-wk)*(1.-wa)
6474 call q_nearest(ik2,ia2,w1k2,w2k2,w3k2,w4k2)
6475 call q_nearest(ik4,ia4,w1k4,w2k4,w3k4,w4k4)
6511 abs(
quad_zz(kmem,amem,iloc)) > 1.e-15)
then 6611 if(
allocated(xloc))
deallocate(xloc,yloc)
6616 if(.not. lwrite)
then 6617 close(
luq_bqf,status=
'delete')
6621 'Q_MAKEGRID: Grid files ',trim(
aqname),
' and ',trim(
bqname), &
6624 'Q_MAKEGRID: Since an error occurred during the generation' 6625 write(
luq_log,
'(a)')
'Q_MAKEGRID: of the interaction grid' 6634 'Q_MAKEGRID: Total number of points on loci :',nztot2
6636 'Q_MAKEGRID: Total number of stored points on locus:',nztot1
6638 'Q_MAKEGRID: Total number of zero points on locus :', &
6641 'Q_MAKEGRID: Reduction factor (%):',real(nztot2-nztot1)/ &
6763 real,
allocatable :: sold(:)
6764 real,
allocatable :: snew(:)
6806 if(abs(
q)>q_eps) nold = nold+1
6820 if(itest>=1)
write(
luq_tst,
'(a,2i4)') &
6821 'Q_MODIFY nold nnew:',
nlocus1,nnew
6823 allocate (sold(nold),snew(nnew))
6828 if(abs(
q)<q_eps)
then 6834 sold(iold) =
s_loc(iold)
6835 slen = slen +
ds_loc(iold)
6857 if(
iq_gauleg > nnew) stop
'Q_MODIFY: iq_gauleg > nlocus0' 6862 write(
luq_tst,
'(a,2f10.4,i4)') &
6863 'Q_MODIFY: GAULEG x1,x2,n:',zero,slen,nnew
6864 write(
luq_tst,
'(a)')
'Q_MODIFY: Gauss-Legendre spacing' 6865 write(
luq_tst,
'(10f12.4)') (snew(inew),inew=1,nnew)
6866 write(
luq_tst,
'(a)')
'Q_MODIFY: Gauss-Legendre weights' 6870 if(abs(
q)>q_eps)
then 6871 dsnew = slen/real(nnew)
6873 snew(inew) = (inew-1.)*dsnew
6876 dsnew = slen/real(nnew-1.)
6878 snew(inew) = (inew-1)*dsnew
6885 write(
luq_tst,
'(a,2f12.5)')
'Q_MODIFY: Slen q:',slen,
q 6886 write(
luq_tst,
'(a,i4)')
'Q_MODIFY: nold /sold:',nold
6887 write(
luq_tst,
'(10f12.6)') sold
6888 write(
luq_tst,
'(a,i4)')
'Q_MODIFY: nnew /snew:',nnew
6889 write(
luq_tst,
'(10f12.6)') snew
6890 write(
luq_tst,
'(a)')
'Q_MODIFY: x2_loc' 6892 write(
luq_tst,
'(a)')
'Q_MODIFY: y2_loc' 6901 if(abs(
q)<1.e-5)
then 6903 if(ierr > 0)
write(
luq_err,*)
'Z_INTP1 x_loc, ierr=',ierr
6907 if(ierr > 0)
write(
luq_err,*)
'Z_INTP1 x_loc, ierr=',ierr
6911 if(ierr > 0)
write(
luq_err,*)
'Z_INTP1 x_loc, ierr=',ierr
6915 if(ierr > 0)
write(
luq_err,*)
'Z_INTP1 x_loc, ierr=',ierr
6919 if(ierr > 0)
write(
luq_err,*)
'Z_INTP1 x_loc, ierr=',ierr
6927 diold = slen/real(nold)
6928 dinew = slen/real(nnew)
6933 jloc = floor((iloc-1.)*diold/dinew)+1
6937 write(
luq_tst,
'(a,2i4,f8.3,3e12.4,f4.0,e12.4)') &
6938 'Q_MODIFY Q=0 iloc,jloc s jac cple ds sym ds_mod:',&
6951 if(ierr > 0)
write(
luq_err,*)
'Z_INTP1 jac_loc, ierr=',ierr
6955 if(ierr > 0)
write(
luq_err,*)
'Z_INTP1 cp_loc, ierr=',ierr
6963 if(ierr > 0)
write(
luq_err,*)
'Z_INTP1 x_loc, ierr=',ierr
6967 if(ierr > 0)
write(
luq_err,*)
'Z_INTP1 y_loc, ierr=',ierr
6971 if(ierr > 0)
write(
luq_err,*)
'Z_INTP1 x_loc, ierr=',ierr
6975 if(ierr > 0)
write(
luq_err,*)
'Z_INTP1 y_loc, ierr=',ierr
6979 if(ierr > 0)
write(
luq_err,*)
'Z_INTP1 s_loc, ierr=',ierr
6983 write(
luq_tst,
'(a)')
'Q_MODIFY: s_loc' 6985 write(
luq_tst,
'(a)')
'Q_MODIFY: s_mod' 6994 diold = slen/real(nold-1)
6995 dinew = slen/real(nnew)
7000 jloc = floor((iloc-1.)*diold/dinew + 1.49999)
7001 jloc = mod(jloc-1+nnew,nnew)+1
7005 write(
luq_tst,
'(a,2i4,f8.3,3e12.4,f4.0,e12.4)') &
7006 'Q_MODIFY: q>0: iloc,jloc, s jac cple ds sym ds_mod:', &
7019 if(ierr > 0)
write(
luq_err,*)
'Z_INTP1 jac_loc, ierr=',ierr
7023 if(ierr > 0)
write(
luq_err,*)
'Z_INTP1 cp_loc, ierr=',ierr
7029 call q_error(
'e',
'INTER',
'Problem in interpolation process')
7041 write(
luq_tst,
'(a)')
'Q_MODIFY: x2_mod' 7043 write(
luq_tst,
'(a)')
'Q_MODIFY: y2_mod' 7045 write(
luq_tst,
'(a)')
'Q_MODIFY: x4_mod' 7047 write(
luq_tst,
'(a)')
'Q_MODIFY: y4_mod' 7049 write(
luq_tst,
'(a)')
'Q_MODIFY: s_mod' 7051 write(
luq_tst,
'(a)')
'Q_MODIFY: ds_loc' 7053 write(
luq_tst,
'(a)')
'Q_MODIFY: ds_mod' 7069 k2m = sqrt(
k2x**2 +
k2y**2)
7070 k4m = sqrt(
k4x**2 +
k4y**2)
7084 write(
luq_tst,
'(a)')
'Q_MODIFY: k2m_mod' 7086 write(
luq_tst,
'(a)')
'Q_MODIFY: k2a_mod' 7088 write(
luq_tst,
'(a)')
'Q_MODIFY: k4m_mod' 7090 write(
luq_tst,
'(a)')
'Q_MODIFY: k4a_mod' 7092 write(
luq_tst,
'(a)')
'Q_MODIFY: sym_mod' 7098 if(
allocated(sold))
deallocate(sold,snew)
7171 integer,
intent(inout) :: ik
7172 integer,
intent(inout) :: ia
7173 real,
intent(inout) :: w1
7174 real,
intent(inout) :: w2
7175 real,
intent(inout) :: w3
7176 real,
intent(inout) :: w4
7202 'Q_NEAREST-A:',ik,ia,w1,w2,w3,w4
7209 if(w1 >= w_max)
then 7216 if(w2 >= w_max)
then 7223 if(w3 >= w_max)
then 7230 if(w4 >= w_max)
then 7250 if(iw_max==1) w1 = 1.0
7251 if(iw_max==2) w2 = 1.0
7252 if(iw_max==3) w3 = 1.0
7253 if(iw_max==4) w4 = 1.0
7256 write(
luq_tst,
'(a,2i4,4f10.5)')
'Q_NEAREST-B:',ik,ia,w1,w2,w3,w4
7263 subroutine q_polar2(kmin,kmax,kx_beg,ky_beg,kx_end,ky_end,loclen, &
7325 real,
intent(in) :: kmin
7326 real,
intent(in) :: kmax
7327 real,
intent(in) :: kx_beg
7328 real,
intent(in) :: ky_beg
7329 real,
intent(in) :: kx_end
7330 real,
intent(in) :: ky_end
7331 real,
intent(in) :: loclen
7332 integer,
intent (out) :: ierr
7416 dk = (kmax-kmin)/real(npol-1)
7418 k_pol(ipol) = kmin + (ipol-1)*dk
7437 dk0 = (kmax - kmin)/real(npol)
7445 'Q_POLAR2: loclen dsz dk:',loclen,dsz,dk
7447 do while (
k_pol(ipol) < kmax .and. iend==0 .and. ipol < mpol)
7449 knew = min(kmax,
k_pol(ipol)+dk)
7450 dkold = knew -
k_pol(ipol)
7452 ang1 =
pang + acos(cosold)
7453 ang2 =
pang + acos(cosnew)
7456 arg = kk1**2 + kk2**2 -2.*kk1*kk2*cos(ang1-ang2)
7457 dsnew = sqrt(abs(arg))
7458 if(dsnew>0) dke = dk*dsz/dsnew
7466 c_pol(ipol) = cosnew
7469 if (abs(dkold) < 0.0005*(kmax-kmin)) iend=1
7474 if(
k_pol(ipol) < kmax .and. ipol < mpol)
then 7488 kratio = (kmax/kmin)**(1./(npol-1.))
7490 'Q_POLAR2: npol kmin kmax kratio:',npol,kmin,kmax,kratio
7492 k_pol(ipol) = kmin*kratio**(ipol-1.)
7527 write(
luq_tst,
'(a,2f12.6,i4)') &
7528 'Q_POLAR2: kmin kmax iq_locus :',kmin,kmax,
iq_locus 7532 'Q_POLAR2: kratio :',kratio
7534 write(
luq_tst,
'(a,i4,4f13.7)') &
7535 'Q_POLAR2: i k(i) a(i) x(i) y(i):',ipol,
k_pol(ipol),&
7628 integer,
intent(in) :: iquad
7658 character(len=10) cpar
7710 elseif(iquad==2)
then 7718 elseif(iquad==3)
then 7725 'Q_SETCONFIG: iquad=',iquad
7727 'No valid value of iquad has been given, default settings')
7728 write(
luq_err,
'(a,i4)')
'Q_SETCONFIG: Value of IQUAD:',iquad
7755 write(
luq_log,
'(a)')
'Q_SETCONFIG: Configuration file '// &
7756 trim(
qbase)//
'.cfg has been found' 7757 write(
luq_log,
'(a,i4)')
'Q_SETCONFIG: '//trim(
qbase)// &
7764 read(
luq_cfg,*,iostat=iend) cpar,rpar
7771 if(trim(cpar)==
'DSTEP')
q_dstep = rpar
7772 if(trim(cpar)==
'F_DMAX')
qf_dmax = rpar
7773 if(trim(cpar)==
'F_KRAT')
qf_krat = rpar
7774 if(trim(cpar)==
'FF_TAIL')
ff_tail = rpar
7775 if(trim(cpar)==
'F_FRAC')
qf_frac = rpar
7778 if(trim(cpar)==
'NLOCUS')
nlocus0 = int(rpar)
7779 if(trim(cpar)==
'SECTOR')
q_sector = rpar
7781 if(trim(cpar)==
'GEOM')
then 7786 'Q_SETCONFIG: geometric scaling disabled' 7788 'Q_SETCONFIG: geometric scaling disabled' 7791 if(trim(cpar)==
'COMPACT')
iq_compact = int(rpar)
7792 if(trim(cpar)==
'COUPLING')
iq_cple = int(rpar)
7793 if(trim(cpar)==
'DISPER')
iq_disp = int(rpar)
7794 if(trim(cpar)==
'DSCALE')
iq_dscale = int(rpar)
7795 if(trim(cpar)==
'FILT')
iq_filt = int(rpar)
7796 if(trim(cpar)==
'GAULEG')
iq_gauleg = int(rpar)
7797 if(trim(cpar)==
'GRID')
iq_grid = int(rpar)
7798 if(trim(cpar)==
'INTEG')
iq_integ = int(rpar)
7799 if(trim(cpar)==
'INTERP')
iq_interp = int(rpar)
7800 if(trim(cpar)==
'LOCUS')
iq_locus = int(rpar)
7801 if(trim(cpar)==
'LOGGING')
iq_log = int(rpar)
7802 if(trim(cpar)==
'LUMPING')
iq_lump = int(rpar)
7803 if(trim(cpar)==
'MAKE')
iq_make = int(rpar)
7804 if(trim(cpar)==
'MODIFY')
iq_mod = int(rpar)
7805 if(trim(cpar)==
'PRINT')
iq_prt = int(rpar)
7806 if(trim(cpar)==
'SCREEN')
iq_screen = int(rpar)
7807 if(trim(cpar)==
'SEARCH')
iq_search = int(rpar)
7808 if(trim(cpar)==
'SYM')
iq_sym = int(rpar)
7809 if(trim(cpar)==
'T13')
iq_t13 = int(rpar)
7810 if(trim(cpar)==
'TAIL')
iq_tail = int(rpar)
7811 if(trim(cpar)==
'TEST')
iq_test = int(rpar)
7812 if(trim(cpar)==
'TRACE')
iq_trace = int(rpar)
7813 if(trim(cpar)==
'TRANSF')
iq_trf = int(rpar)
7814 if(trim(cpar)==
'XDIA')
iq_xdia = int(rpar)
7821 'Q_SETCONFIG: '//trim(
qbase)//
'.cfg disconnected from :', &
7828 write(
luq_log,
'(a)')
'Q_SETCONFIG: Configuration file '// &
7829 trim(
qbase)//
'.CFG has not been found' 7896 real,
intent(in) :: depth
7897 integer,
intent(out) :: igrid
7953 'Q_SEARCHGRID: Input target depth:',depth
7960 'Q_SEARCHGRID: First call of Q_CTRGGRID exit code & q_depth:', &
7965 write(
luq_prt,
'(a,f10.2)')
'Q_SEARCHGRID: target depth:', &
7968 'Q_SEARCHGRID: grid accepted, read whole database' 7971 'Q_SEARCHGRID: grid accepted, read whole database' 7977 idepth = int(s_depth*10+eps)
7981 id_upper = min(
id_facmax*idepth,id_upper)
7989 'Q_SEARCHGRID: idepth,id_lower/upper:', &
7990 idepth,id_lower,id_upper
7995 do id = idepth-1,id_lower,-1
7998 'Q_SEARCHGRID: downwards id q_depth:',id,
q_depth 8003 'Q_SEARCHGRID: igrid:',igrid
8007 'Q_SEARCHGRID: valid grid found for depth:',
q_depth 8017 do id = idepth+1,id_upper
8021 'Q_SEARCHGRID: upwards id q_depth:',id,
q_depth 8026 'Q_SEARCHGRID: igrid:',igrid
8030 'Q_SEARCHGRID: valid grid found for depth:',
q_depth 8041 if(d_lower > 0)
then 8042 r_lower = s_depth/d_lower
8047 if(d_upper > 0)
then 8048 r_upper = d_upper/s_depth
8055 'Q_SEARCHGRID: d_lower d_target d_upper :', &
8056 d_lower,s_depth,d_upper
8058 'Q_SEARCHGRID: r_lower r_upper :', &
8064 if(r_lower>0 .and. r_upper>0)
then 8065 if(r_lower < r_upper)
then 8071 elseif(r_lower > 0 .and. r_upper <0 )
then 8073 elseif(r_lower < 0 .and. r_upper > 0)
then 8076 call q_error(
'e',
'SEARCHGRID', &
8077 'No valid nearest grid could be found')
8082 'Q_SEARCHGRID: target and nearest water depth :',s_depth,
q_depth 8084 'Q_SEARCHGRID: nearest valid BQF depth:',
q_depth 8096 'Q_SEARCHGRID: target and nearest scale factors:',dfac1,dfac2
8098 'Q_SEARCHGRID: compound scale factor :',
q_scale 8106 'Q_SEARCHGRID: Q_CTRGRID called with depth:',
q_depth 8108 'Q_SEARCHGRID: igrid of nearest grid operation:',igrid
8116 write(*,*)
'q_searchgrid q_depth on exit:',
q_depth 8131 'GurboQuad Version: 5.03 Build: 120 Date: 2004/05/07 [ST]' 8198 character(len=*),
intent(in) :: mod_name
8218 character(len=1) mod_task
8223 if(
iq_prt>0)
write(
luq_prt,
'(2a)')
'TRACE -> ',trim(mod_name)
8232 mod_len = len_trim(mod_name)
8233 mod_task = mod_name(1:1)
8236 if(mod_task(1:1) ==
'+')
then 8240 call q_error(
'e',
'STACKMAX',
' ')
8248 elseif(mod_task(1:1) ==
'-')
then 8253 write(
luq_err,
'(a)')
'Module name:',mod_name
8254 call q_error(
'e',
'STACKNAME',
' ')
8258 call q_error(
'e',
'STACKCALL',
' ')
8359 write(
luq_prt,
'(a)')
'Summary of settings for QUAD computation' 8361 '------------------------------------------------' 8362 write(
luq_prt,
'(a,i4)')
'Number of wave numbers :', &
8364 write(
luq_prt,
'(a,i4)')
'Number of directions :', &
8366 write(
luq_prt,
'(a,f10.5)')
'Minimum frequency (Hz) :', &
8368 write(
luq_prt,
'(a,f10.5)')
'Maximum frequency (Hz) :', &
8370 write(
luq_prt,
'(a,f10.2)')
'Water depth (m) :', &
8372 write(
luq_prt,
'(a,i4)')
'Preferred number of locus points:', &
8376 write(
luq_prt,
'(a,f10.3)')
'Gravitational acceleration:',
q_grav 8383 'IQUAD = 1: Deep water' 8385 'IQUAD = 2: Deep water & WAM depth scaling' 8387 'IQUAD = 3: Direct finite depth calculation' 8390 write(
luq_prt,
'(a,f5.2)')
'Step size in m of BQF coding:', &
8397 'Non-symmetric full circle grid' 8401 'No compacting of data along locus' 8403 'Compact data along locus by eliminating zero contributions' 8408 'WAM depth scaling of transfer' 8413 'Intermediate output to screen' 8415 'Intermediate output to screen + subroutine tracing' 8420 'No search is carried out for nearest QUAD grid' 8422 'A search is carried out for nearest QUAD grid' 8427 'Gauss-Legendre integration with N=',
iq_gauleg 8431 'Deep water coupling coefficient of Webb' 8433 'Finite depth coupling coefficient of H&H' 8435 'Finite depth coupling coefficient of Gorman' 8437 'Deep water coefficient of Zakharov' 8439 'Finite depth coefficient of Zakharov' 8443 'Deep water dispersion relation' 8445 'Finite depth linear dispersion relation' 8447 'Non linear finite depth dispersion' 8451 'Filtering of quadruplets off' 8453 write(
luq_prt,
'(a)')
'Filtering of quadruplets on' 8456 'Maximum ratio of k1 and k3 :',
qf_krat 8458 'Maximum directional difference :',
qf_dmax 8460 'Fraction of maximum energy density:',
qf_frac 8469 'Compute locus with polar method with fixed k-step' 8471 'Compute locus with polar method using adaptive k-step' 8473 'Compute locus with polar method using geometric k-step' 8477 'Handling of symmetries disabled' 8479 'Handling of symmetries enabled' 8483 'Make quadruplet grid when necessary' 8485 'Always make quadruplet grid' 8487 'Stop after generation of quadruplet grid' 8491 'Apply bi-linear interpotion to retrieve action density' 8493 'Take nearest bin to retrieve action density' 8497 'Lumping of coefficients along locus disabled' 8499 'Lumping of coefficients along locus enabled' 8503 '?X? Spacing of point along locus as initially computed' 8505 'Equidistant spacing of points along locus' 8509 'No parametric tail is added' 8511 'Parametric tail is added from ', &
8512 ff_tail,
' times maximum frequency' 8516 'Subroutine tracing disabled' 8518 'Subroutine tracing enabled' 8524 'No test output of integration' 8526 'Summary output of integration per locus' 8528 'Extended output of integration along locus' 8533 'No test output of T13 integration' 8535 'Summary output of T13 integration' 8545 write(
luq_prt,
'(a,i4)')
'Level of printed output :', &
8547 write(
luq_prt,
'(a,i4)')
'Level of logging output :', &
8549 write(
luq_prt,
'(a,i4)')
'Level of test output :', &
8551 write(
luq_prt,
'(a,i4)')
'Level of trace output :', &
8553 write(
luq_prt,
'(a,i4)')
'Level of transformation output :', &
8556 '----------------------------------------------' 8566 subroutine q_symmetry(k1x,k1y,k3x,k3y,k4x,k4y,symfac,nloc)
8615 integer,
intent(in) :: nloc
8616 real,
intent(in) :: k1x
8617 real,
intent(in) :: k1y
8618 real,
intent(in) :: k3x
8619 real,
intent(in) :: k3y
8620 real,
intent(in) :: k4x(nloc)
8621 real,
intent(in) :: k4y(nloc)
8622 real,
intent(out) :: symfac(nloc)
8655 dk13 = (k1x-k3x)**2 + (k1y-k3y)**2
8657 dk14 = (k1x-k4x(iloc))**2 + (k1y-k4y(iloc))**2
8658 if (dk13 >= dk14) symfac(iloc) = 0.
8667 subroutine q_t13v4(ik1,ia1,ik3,ia3,t13,diagk1,diagk3)
8736 integer,
intent(in) :: ik1
8737 integer,
intent(in) :: ia1
8738 integer,
intent(in) :: ik3
8739 integer,
intent(in) :: ia3
8740 real,
intent(out) :: t13
8741 real,
intent(out) :: diagk1
8742 real,
intent(out) :: diagk3
8789 real qn1,qn2,qn3,qn4
8811 'Q_T13V4: ik1,ia1 ik3 ia3:',ik1,ia1,ik3,ia3
8813 if(ik1==ik3 .and. ia1==ia3)
goto 9999
8821 if(ifnd==0 .or.
nlocusx==0)
then 8827 qn1 =
nspec(ik1,ia1)
8828 qn3 =
nspec(ik3,ia3)
8852 write(
luq_int,
'(a1,2i3.3,a1,2i3.3,a1)')
'(',ik1,ia1,
'-',ik3, &
8856 q_k(ik1)*sin(
q_a(ia1)), &
8857 q_k(ik3)*cos(
q_a(ia3)), &
8868 qn2 =
nspec(jk2,ja2)
8869 qn4 =
nspec(jk4,ja4)
8870 nprod = qn13p*(qn4-qn2) + qn2*qn4*qn13d
8872 t13 = t13 + nprod*rterm
8876 qd1 = qn3*(qn4-qn2) + qn2*qn4*qn3
8877 qd3 = qn1*(qn4-qn2) - qn2*qn4*qn1
8878 diagk1 = diagk1 + qd1*rterm
8879 diagk3 = diagk3 + qd3*rterm
8890 jk2p = min(jk2+1,
nkq)
8917 ja2p = min(ja2p,
naq)
8926 jk4p = min(jk4+1,
nkq)
8948 ja4p = min(ja4p,
naq)
8958 nprod = qn13p*(qn4-qn2) + qn2*qn4*qn13d
8960 t13 = t13 + rterm*nprod
8965 dt13(iloc) = nprod *rterm
8972 qd1 = qn3*(qn4-qn2) - qn2*qn4
8973 qd3 = qn1*(qn4-qn2) + qn2*qn4
8974 diagk1 = diagk1 + qd1*rterm
8975 diagk3 = diagk3 + qd3*rterm
8979 if(jq_int == 1)
then 8995 write(
luq_int,
'(1x,i4,5f8.3,11e12.4,2f11.6)') &
8997 t_ws(iloc),qn1,qn2,qn3,qn4,nprod, &
9010 if(
iq_integ==3 .and. ik1>=30 .and. ik1<=35 .and. ik3>=30 .and. &
9012 write(
luq_int,
'(a,4i3,i5,e13.5)')
'Q_T13V4:',ik1,ia1,ik3,ia3, &
9197 do while (k2m >
q_k(jpos))
9202 if(k2m <=
q_k(1))
then 9205 elseif(k2m <
q_k(
nkq) .and. k2m >
q_k(1))
then 9206 dk =
q_k(jpos)-
q_k(jpos-1)
9207 wk_k2(iloc) = real(jpos-1) + (k2m-
q_k(jpos-1))/dk
9209 elseif(k2m >=
q_k(
nkq))
then 9221 do while (k4m >
q_k(jpos))
9226 if(k4m <=
q_k(1))
then 9229 elseif(k4m <
q_k(
nkq) .and. k4m >
q_k(1))
then 9230 dk =
q_k(jpos)-
q_k(jpos-1)
9231 wk_k4(iloc) = real(jpos-1) + (k4m-
q_k(jpos-1))/dk
9233 elseif(k4m >=
q_k(
nkq))
then 9240 if(itest >= 1)
write(
luq_tst,
'(a,i3,10f12.4)') &
9241 'Q_WEIGHT: i k2m k2a wk2 wa2 wt2(+4):', &
9254 subroutine q_loc_w1w3(k1x,k1y,k3x,k3y,npts,k2x,k2y,k4x,k4y,s)
9303 integer,
intent(in) :: npts
9304 real,
intent(in) :: k1x
9305 real,
intent(in) :: k1y
9306 real,
intent(in) :: k3x
9307 real,
intent(in) :: k3y
9309 real,
intent(out) :: k2x(npts)
9310 real,
intent(out) :: k2y(npts)
9311 real,
intent(out) :: k4x(npts)
9312 real,
intent(out) :: k4y(npts)
9313 real,
intent(out) :: s(npts)
9360 dir1 = atan2(k1y,k1x)
9361 dir3 = atan2(k3y,k3x)
9362 dirs = 0.5*(180-abs(180-abs(dir3-dir1)))
9363 k1m = sqrt(k1x**2 + k1y**2)
9368 xk0 = k1m * cos(dirs)
9369 yk0 = k1m * sin(dirs)
9375 dk0 = 3./real(npts-1.)
9387 w2 = 2.*real(ipt-npts/2)*dk0
9390 k2x(ipt) = xx2*cos(dirs) - yy2*sin(dirs)
9391 k2y(ipt) = yy2*cos(dirs) + xx2*sin(dirs)
9394 k4x(ipt) = xx4*cos(dirs) - yy4*sin(dirs)
9395 k4y(ipt) = yy4*cos(dirs) + xx4*sin(dirs)
9396 s(ipt) = real(ipt-1)*dk0*xk0
9402 subroutine q_xnl4v4(aspec,sigma,angle,nsig,nang,depth,xnl,diag, &
9477 integer,
intent(in) :: nsig
9478 integer,
intent(in) :: nang
9479 real,
intent(in) :: aspec(nsig,nang)
9480 real,
intent(in) :: sigma(nsig)
9481 real,
intent(in) :: angle(nang)
9482 real,
intent(in) :: depth
9483 real,
intent(out) :: xnl(nsig,nang)
9485 real,
intent(out) :: diag(nsig,nang)
9486 integer,
intent(out) :: ierror
9590 'Q_XNL4V4: Checking interaction grid ' 9603 call q_error(
'e',
'NOGRID',
'No proper grid exists')
9609 call q_error(
'e',
'MAKEGRID',
'Only computation of grid')
9624 'Q_XNL4V4: Q_SEARCHGRID called with q_depth: ',
q_depth 9628 'Q_XNL4V4: Q_SEARCHGRID exited with q_depth: ',
q_depth 9631 call q_error(
'e',
'NOGRID',
'No proper grid exists')
9645 nspec(ikq,iaq) = aspec(ikq,iaq)/
q_k(ikq)*cg(ikq)
9666 qn_max = maxval(
nspec)
9677 qn1 =
nspec(ik1,ia1)
9681 qn3 =
nspec(ik3,ia3)
9684 'Q_XNL4V4: ik1 ia1 ik3 ia3:',ik1,ia1,ik3,ia3
9688 a_dif = 180. - abs(180. - abs(
q_ad(ia1) -
q_ad(ia3)))
9714 if(qn1 < qn_min .and. qn3 < qn_min)
then 9719 if(ifil_dir==0 .and. ifil_krat==0 .and. ifil_dens==0 .or. &
9725 call q_t13v4(ik1,ia1,ik3,ia3,t13,diagk1,diagk3)
9765 idir1 = int(
q_ad(ia1))
9766 idir3 = int(
q_ad(ia3))
9767 if (idir1>180) idir1 = idir1-360
9768 if (idir3>180) idir3 = idir3-360
9769 icx1 = (idir1/15-1)*20+ik1
9770 icx3 = (idir3/15-1)*20+ik3
9771 xt13 = alog10(max(1.e-20,abs(t13)))
9772 if(xt13<-19.99) xt13=0
9773 write(
luq_t13,
'(4i6,e13.5,2i6,f10.4)') ik1,idir1,ik3, &
9774 idir3,t13,icx1,icx3,xt13
9797 xnl(ik1,ia1) = xnl(ik1,ia1) + t13*
q_k(ik3)* &
9799 xnl(ik3,ja3) = xnl(ik3,ja3) - t13*
q_k(ik1)* &
9804 diag(ik1,ia1) = diag(ik1,ia1) + diagk1*
q_k(ik3)* &
9806 diag(ik3,ja3) = diag(ik3,ja3) - diagk3*
q_k(ik1)* &
9842 xnl(ikq,jaq) = xnl(ikq,iaq)
9843 diag(ikq,jaq) = diag(ikq,iaq)
9850 'Q_XNL4V4: Main computation ended' 9855 jacobian =
q_k(ikq)/cg(ikq)
9857 xnl(ikq,iaq) = xnl(ikq,iaq)*jacobian
9868 real function x_cosk(k)
9920 real,
intent(in) :: k
10025 real,
intent(in) ::
k1x 10026 real,
intent(in) ::
k1y 10027 real,
intent(in) ::
k2x 10028 real,
intent(in) ::
k2y 10029 real,
intent(in) ::
k3x 10030 real,
intent(in) ::
k3y 10031 real,
intent(in) ::
k4x 10032 real,
intent(in) ::
k4y 10033 integer,
intent(in) ::
iq_cple 10034 real,
intent(in) :: depth
10035 real,
intent(in) :: grav_w
10143 real,
intent(in) :: kxx
10144 real,
intent(in) :: kyy
10171 real (kind=8) w2,w4
10178 w2 =
sqrtg * (kxx**2 + kyy**2)**(0.25)
10179 w4 =
sqrtg * ((kxx+
px)**2 + (kyy+
py)**2)**(0.25)
10183 k2m = sqrt(kxx**2+kyy**2)
10184 k4m = sqrt((kxx+
px)**2 + (kyy+
py)**2)
10261 real,
intent(in) :: x2
10262 real,
intent(in) :: y2
10263 real,
intent(in) :: x4
10264 real,
intent(in) :: y4
10290 k2m = sqrt(x2**2 + y2**2)
10291 k4m = sqrt(x4**2 + y4**2)
10293 ang2 = atan2(x2,y2)
10294 ang4 = atan2(x4,y4)
10305 cg2 = sig2/k2m*(0.5+k2md/sinh(2*k2md))
10311 cg4 = sig4/k4m*(0.5+k4md/sinh(2*k4md))
10314 x_jacobian = sqrt(cg2**2+cg4**2-2*cg2*cg4*cos(ang2-ang4))
10374 real,
intent(in) :: k
10375 real,
intent(in) :: d
10409 if (kd > 20.) id = 1
10427 real function
xc_hh(w1x0,w1y0,w2x0,w2y0,w3x0,w3y0,z4x,z4y,h)
10436 real w1x0,w1y0,w2x0,w2y0,w3x0,w3y0,h,dsq
10437 real om1,om2,om3,om4,scpl1,scpl2,scpl3,stot
10438 real t1,t2,t3,t4,t5,tot1,tot2,tot3,tot4,tot5
10439 real som1,som2,som3
10440 real s1,s2,s3,z1,z2,z3,z4,z5
10441 real p1,p2,p3,p4,di,tnz1,tnz2,tnz3,tnz23
10442 real csz1,csz2,csz3,csz23
10443 real e,g,
gsq,omsq23,pi4
10451 real k1,k2,k3,
k1x,
k2x,
k3x,
k1y,
k2y,
k3y,k23x,k23y,k23,k1x0,k1y0, &
10452 k2x0,k2y0,k3x0,k3y0,k1zx,k1zy
10453 data pi4/0.785398163/
10504 om1=sqrt(g*k1*tnz1)
10505 om2=sqrt(g*k2*tnz2)
10506 om3=sqrt(g*k3*tnz3)
10515 k23=sqrt(k23x**2+k23y**2)
10520 dot123=
k1x*k23x+
k1y*k23y
10523 di=-(som2+som3)*(k2*k3*tnz2*tnz3-dot23) &
10524 +0.5*(som2*k3**2/(csz3)**2+som3*k2**2/(csz2)**2)
10526 e=0.5/g *(dot23-som2*som3/
gsq*(om2**2+om3**2+som2*som3))
10528 p1=2.*(som1+som2+som3)*(om1**2.*omsq23/
gsq-dot123)
10530 p2=-som1*(k23)**2/(csz23)**2
10532 p3=-(som2+som3)*k1**2/(csz1)**2
10535 z2=z2+omsq23-(som2+som3)**2
10539 t1=di/(omsq23-(som2+som3)**2 + eps ) * (p1+p2+p3)
10541 t2=-di*som1/
gsq *(om1**2+omsq23)
10543 p4=g*k1**2/(csz1)**2
10545 t3=e*(som1**3*(som2+som3)/g - g*dot123 - p4)
10547 t4=0.5*som1/
gsq*dot23*((som1+som2+som3)*(om2**2+om3**2) &
10548 +som2*som3*(som2+som3))
10550 t5=-0.5*som1*om2**2*k3**2/
gsq*(som1+som2+2.*som3) &
10551 -0.5*som1*om3**2*k2**2/
gsq*(som1+2.*som2+som3)
10553 scpl1=t1+t2+t3+t4+t5
10589 om1=sqrt(g*k1*tnz1)
10590 om2=sqrt(g*k2*tnz2)
10591 om3=sqrt(g*k3*tnz3)
10598 k23=sqrt(k23x**2+k23y**2)
10602 dot123=
k1x*k23x+
k1y*k23y
10605 di=-(som2+som3)*(k2*k3*tnz2*tnz3-dot23) &
10606 +0.5*(som2*k3**2/(csz3)**2+som3*k2**2/(csz2)**2)
10608 e=0.5/g *(dot23-som2*som3/
gsq *(om2**2+om3**2+som2*som3))
10610 p1=2.*(som1+som2+som3)*(om1**2.*omsq23/
gsq-dot123)
10612 p2=-som1*(k23)**2/(csz23)**2
10614 p3=-(som2+som3)*k1**2/(csz1)**2
10616 z2=z2+omsq23-(som2+som3)**2
10621 t1=di/(omsq23-(som2+som3)**2) * (p1+p2+p3)
10623 t2=-di*som1/
gsq *(om1**2+omsq23)
10625 p4=g*k1**2/(csz1)**2
10627 t3=e*(som1**3*(som2+som3)/g - g*dot123 - p4)
10629 t4=0.5*som1/
gsq*dot23*((som1+som2+som3)*(om2**2+om3**2) &
10630 +som2*som3*(som2+som3))
10632 t5=-0.5*som1*om2**2*k3**2/
gsq*(som1+som2+2.*som3) &
10633 -0.5*som1*om3**2*k2**2/
gsq*(som1+2.*som2+som3)
10635 scpl2=t1+t2+t3+t4+t5
10671 om1=sqrt(g*k1*tnz1)
10672 om2=sqrt(g*k2*tnz2)
10673 om3=sqrt(g*k3*tnz3)
10680 k23=sqrt(k23x**2+k23y**2)
10684 dot123=
k1x*k23x+
k1y*k23y
10687 di=-(som2+som3)*(k2*k3*tnz2*tnz3-dot23) &
10688 +0.5*(som2*k3**2/(csz3)**2+som3*k2**2/(csz2)**2)
10690 e=0.5/g *(dot23-som2*som3/
gsq *(om2**2+om3**2+som2*som3))
10692 p1=2.*(som1+som2+som3)*(om1**2.*omsq23/
gsq-dot123)
10694 p2=-som1*(k23)**2/(csz23)**2
10696 p3=-(som2+som3)*k1**2/(csz1)**2
10698 z2=z2+omsq23-(som2+som3)**2
10703 t1=di/(omsq23-(som2+som3)**2) * (p1+p2+p3)
10705 t2=-di*som1/
gsq*(om1**2+omsq23)
10707 p4=g*k1**2/(
cosz(k1*h))**2
10709 t3=e*(som1**3*(som2+som3)/g - g*dot123 - p4)
10711 t4=0.5*som1/
gsq*dot23*((som1+som2+som3)*(om2**2+om3**2) &
10712 +som2*som3*(som2+som3))
10714 t5=-0.5*som1*om2**2*k3**2/
gsq*(som1+som2+2.*som3) &
10715 -0.5*som1*om3**2*k2**2/
gsq*(som1+2.*som2+som3)
10717 scpl3=t1+t2+t3+t4+t5
10724 stot=(scpl1+scpl2+scpl3)
10726 dsq=stot*stot*pi4*
gsq/(om1*om2*om3*om4+eps)
10736 real function
tanz(x)
10739 if (x.gt.20.) x=25.
10745 real function
cosz(x)
10747 if (x.gt.20.) x=25.
10805 real,
intent(in) ::
k1x 10806 real,
intent(in) ::
k1y 10807 real,
intent(in) ::
k2x 10808 real,
intent(in) ::
k2y 10809 real,
intent(in) ::
k3x 10810 real,
intent(in) ::
k3y 10811 real,
intent(in) ::
k4x 10812 real,
intent(in) ::
k4y 10813 real,
intent(in) :: grav_w
10834 double precision wsqp12
10835 double precision wsqm13
10836 double precision wsq13
10837 double precision wsqm14
10838 double precision wsq14
10839 double precision wsq12
10842 real p1,p2,p3,p4,p5,p6,p7,p8,p9
10880 wsq12 = (w1+w2)*(w1+w2)
10882 wsq13 = (w1-w3)*(w1-w3)
10884 wsq14 = (w1-w4)*(w1-w4)
10888 z = 2.*wsq12*(k1*k2-dot12)*(k3*k4-dot34)
10890 z = 2.*wsq13*(k1*k3+dot13)*(k2*k4+dot24)
10892 z = 2.*wsq14*(k1*k4+dot14)*(k2*k3+dot23)
10894 p4 = 0.5 *(dot12*dot34 + dot13*dot24 + dot14*dot23)
10895 p5 = 0.25*(dot13+dot24) * wsq13 * wsq13
10896 p6 = -0.25*(dot12+dot34) * wsq12 * wsq12
10897 p7 = 0.25*(dot14+dot23) * wsq14 * wsq14
10898 p8 = 2.5*k1*k2*k3*k4
10899 p9 = wsq12*wsq13*wsq14* (k1 + k2 + k3 + k4)
10901 dwebb = p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9
10902 xc_webb = grav_w**2*pi4*dwebb*dwebb/(w1*w2*w3*w4+eps)
real, dimension(:), allocatable s_mod
real, dimension(:), allocatable t_w1k2
character(len=21), dimension(mq_stack) cstack
real, dimension(:), allocatable cple_mod
subroutine init_constants
real, dimension(:), allocatable q_k
subroutine q_chkcons(xnl, nk, ndir, sum_e, sum_a, sum_mx, sum_my)
subroutine z_flunit(iunit, ierr)
real, dimension(:), allocatable t_tail2
real, dimension(:), allocatable wk_k4
real, dimension(:), allocatable x4_loc
real, dimension(:), allocatable r_w3k2
integer, dimension(:), allocatable t_ia4
real, dimension(:), allocatable wa_k4
subroutine z_polyarea(xpol, ypol, npol, area)
real, dimension(:), allocatable q_sig
real, dimension(:,:,:), allocatable quad_t2
real, dimension(:), allocatable r_w2k2
subroutine z_cmpcg(sigma, depth, grav_w, cg)
subroutine q_addtail(xnl, diag, nsig, na, pf_tail)
real, dimension(:), allocatable dt13
real, dimension(:,:,:), allocatable quad_jac
real, dimension(:), allocatable t_w4k2
real, dimension(:), allocatable t_w4k4
real, dimension(:), allocatable y4_mod
real, dimension(:), allocatable r_w1k2
real, dimension(:), allocatable t_w3k4
real, dimension(:), allocatable wt_k4
real, dimension(:), allocatable r_w4k4
integer, parameter i_print
real function xc_webb(k1x, k1y, k2x, k2y, k3x, k3y, k4x, k4y, grav_w)
subroutine q_dscale(n, sigma, angle, nsig, nang, depth, grav_w, q_dfac)
real, dimension(:), allocatable ds_mod
real function z_root2(func, x1, x2, xacc, iprint, ierr)
real, dimension(:), allocatable q_kpow
real, dimension(:,:,:), allocatable quad_t4
real, dimension(:), allocatable k_pol
real, dimension(:), allocatable k4m_mod
real, dimension(:), allocatable q_dsig
real, dimension(:), allocatable s_loc
real, dimension(:), allocatable a_pol
integer, dimension(:), allocatable r_ik4
real, dimension(:), allocatable r_w4k2
real, dimension(:), allocatable wk_k2
real, dimension(:), allocatable y2_mod
real, dimension(:,:,:), allocatable quad_w2k4
subroutine q_getlocus(ik1, ia1, ik3, ia3, ifnd)
real, dimension(:), allocatable r_sym
real, dimension(:), allocatable x2_loc
real, dimension(:), allocatable t_tail4
integer, dimension(:,:,:), allocatable quad_ia2
integer, dimension(:), allocatable t_ik4
real, dimension(:), allocatable z_loc
real, dimension(:), allocatable r_w1k4
real, dimension(:,:,:), allocatable quad_w4k4
real, dimension(:), allocatable wt_k2
integer, dimension(:,:,:), allocatable quad_ik4
real, dimension(:,:,:), allocatable quad_zz
integer, dimension(:,:,:), allocatable quad_ia4
real function x_disper(k, d)
real, dimension(:), allocatable q_ad
integer, dimension(:), allocatable t_ia2
subroutine q_symmetry(k1x, k1y, k3x, k3y, k4x, k4y, symfac, nloc)
subroutine q_chkres(k1x, k1y, k2x, k2y, k3x, k3y, k4x, k4y, dep, sum_kx, sum_ky, sum_w)
real, dimension(:), allocatable t_zz
character(len=20) qf_error
real, dimension(:,:,:), allocatable quad_cple
subroutine z_intp1(x1, y1, x2, y2, n1, n2, ierr)
real, dimension(:,:), allocatable nspec
real, dimension(:), allocatable r_tail2
subroutine z_fileio(filename, qual, iufind, iunit, iostat)
real, dimension(:), allocatable t_w1k4
subroutine z_steps(x, dx, nx)
integer, dimension(:,:,:), allocatable quad_ik2
real, dimension(:), allocatable k2a_mod
real, dimension(:), allocatable t_sym
subroutine q_stack(mod_name)
real, dimension(:,:,:), allocatable quad_sym
real, dimension(:,:,:), allocatable quad_w1k4
subroutine q_locpos(ka, kb, km, kw, loclen)
real, dimension(:), allocatable y2_loc
character(len=21) q_header
real, dimension(:), allocatable q_a
real, dimension(:), allocatable r_w3k4
real, dimension(:,:,:), allocatable quad_w1k2
subroutine q_cmplocus(ka, kb, km, kw, loclen)
real, dimension(:), allocatable t_cple
integer, parameter mq_stack
real, dimension(:), allocatable z_mod
real, dimension(:,:,:), allocatable quad_w3k4
integer, dimension(:), allocatable r_ia4
real, dimension(:), allocatable x2_mod
subroutine q_xnl4v4(aspec, sigma, angle, nsig, nang, depth, xnl, diag, ierror)
real, dimension(:), allocatable r_tail4
real, dimension(:,:,:), allocatable quad_ws
subroutine xnl_init(sigma, dird, nsigma, ndir, pftail, x_grav, depth, ndepth, jquad, iqgrid, ierror)
real, dimension(:), allocatable r_zz
subroutine q_t13v4(ik1, ia1, ik3, ia3, t13, diagk1, diagk3)
real, dimension(:), allocatable sym_mod
real, dimension(:), allocatable y4_loc
real, dimension(:), allocatable c_pol
real, dimension(:), allocatable k2m_mod
subroutine q_loc_w1w3(k1x, k1y, k3x, k3y, npts, k2x, k2y, k4x, k4y, s)
character(len=21) r_header
real, dimension(:), allocatable ds_loc
real, dimension(:), allocatable q_df
real function x_cple(k1x, k1y, k2x, k2y, k3x, k3y, k4x, k4y, iq_cple, depth, grav_w)
real, dimension(mwind) pwind
real, dimension(:), allocatable jac_mod
integer, dimension(:,:), allocatable quad_nloc
subroutine z_fclose(iunit)
real, dimension(:), allocatable t_jac
subroutine q_ctrgrid(itask, igrid)
real, dimension(:), allocatable t_w2k2
real, dimension(:), allocatable r_cple
real, dimension(:), allocatable sym_loc
character(len=80) tempfile
subroutine q_searchgrid(depth, igrid)
subroutine xnl_main(aspec, sigma, angle, nsig, ndir, depth, iquad, xnl, diag, iproc, ierror)
real function z_wnumb(w, d, grav_w)
subroutine q_polar2(kmin, kmax, kx_beg, ky_beg, kx_end, ky_end, loclen, ierr)
real, dimension(:), allocatable r_w2k4
real function x_jacobian(x2, y2, x4, y4)
real, dimension(:), allocatable q_f
real, dimension(:), allocatable k4a_mod
real, dimension(:), allocatable r_ws
subroutine y_gauleg(x1, x2, x, w, n)
real, dimension(:,:), allocatable qnl
real, dimension(:,:), allocatable a
real, dimension(:,:,:), allocatable quad_w2k2
character(len=60) q_version
subroutine q_setconfig(iquad)
real, dimension(:), allocatable jac_loc
real, dimension(:), allocatable q_dk
subroutine q_nearest(ik, ia, w1, w2, w3, w4)
real, dimension(:), allocatable wa_k2
real function xc_hh(w1x0, w1y0, w2x0, w2y0, w3x0, w3y0, z4x, z4y, h)
real, dimension(:), allocatable q_cg
integer, dimension(:), allocatable r_ik2
real, dimension(:), allocatable r_jac
real, dimension(:), allocatable t_w3k2
real function x_locus2(lambda)
character(len=20) sub_name
real, dimension(:), allocatable x4_mod
real, dimension(:), allocatable t_w2k4
real, dimension(:), allocatable t_ws
real function x_locus1(k2)
real, dimension(:,:,:), allocatable quad_w4k2
real, dimension(:), allocatable q_sk
real, dimension(:), allocatable q_xk
real, dimension(:,:,:), allocatable quad_w3k2
subroutine q_error(err_type, err_name, err_msg)
integer, dimension(:), allocatable r_ia2
real function x_flocus(kxx, kyy)
integer, dimension(:), allocatable t_ik2
real, dimension(:), allocatable cple_loc
character(len=17) lastquadfile
real, dimension(:), allocatable nk1d