module mod_attenuation_ice use mod_xc c --- provides attenuation coefficient, damping, ice wavelength and c --- group velocity; c --- uses 2D (thickness, period) chebyshev interpolation from pre-computed c --- coefficients to get them; c --- TODO this slows down program may need to speed up c --- WAVES_ROLLOVER flag -> try and incorporate effects in attenuation c --- makes no difference though - probably proper wind generation c --- effects are the best way to go; implicit none c !real :: atten_coeff,damping,wavlen,groupvel,t,h !integer :: ACoption contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine attenuation_ice(atten_coeff,damping,groupvel,wavlen, & t,h,ACoption,AGoption) implicit none !! real, intent(in) :: t,h integer, intent(in) :: ACoption integer, intent(in) :: AGoption !! real, intent(out) :: atten_coeff,damping,wavlen,groupvel real, parameter :: pi = 3.1415927 damping = 0.0 wavlen = 0.0 groupvel = 0.0 if (AGoption.eq.2) then wavlen = fn_Kcheb(2*pi/t,h) groupvel = fn_AGcheb(2*pi/t,h) else if (AGoption.eq.1) then wavlen = fn_Kcheb(2*pi/t,h) end if if (ACoption.eq.0) then !Dany's original polynomial approx for ac of Kohout & Meylan; !no draft, no damping; ! (2008); atten_coeff = alpKM08(t,h) !print*, 'hi - using alpKM08' elseif (ACoption.eq.1) then !model 'A1', alpha=-2*log(1-|R|^2); !draft, R=semi-infinite reflection coefficient; !no damping; atten_coeff = alpA1(2*pi/t,h) elseif (ACoption.eq.2) then !model 'AD1', alpha=-2*log(1-|R|^2); !draft, R=semi-infinite reflection coefficient; !extra damping; atten_coeff = alpA1(2*pi/t,h) !!use hard coded damping - visc_RP ~ 10 damping = dampingAD1(2*pi/t,h) elseif (ACoption.eq.3) then !model 'A0', alpha=-2*log(1-|R|^2); !no draft, R=semi-infinite reflection coefficient; !no damping; atten_coeff = alpA0(2*pi/t,h) elseif (ACoption.eq.4) then !model 'E1', ensemble, Toyota distribution; !draft, no damping; atten_coeff = alpE1(2*pi/t,h) !print*, 'hi - using alpE1' elseif (ACoption.eq.5) then !model 'E0', ensemble, long floe, Rayleigh; !draft, no damping; atten_coeff = alpE0(2*pi/t,h) !print*, 'hi - using alpE0' !print*, 't,h=',t,h end if #if defined (WAVES_ROLLOVER) roll_rat = rollover_ratio(2*pi/t) atten_coeff = roll_rat*atten_coeff #endif end subroutine attenuation_ice !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine OP_chebinterp2d(z,x,y,chebys,Ncx,Limits) real, intent(out) :: z real, intent(in) :: x,y integer, intent(in) :: Ncx !! real, intent(in), dimension(:) :: chebys real, intent(in), dimension(4) :: Limits ! real, intent(in) :: chebys, Limits !! real :: Amn, am real :: Tn0_x, Tn1_x, Tn_x, tx real :: Tm0_x, Tm1_x, Tm_x real :: Tn0_y, Tn1_y, Tn_y, ty real :: xmin, xmax, ymin, ymax integer :: s,nx,ny,Ncy !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Ncy = size(chebys)/(Ncx+1)-1! aint(size(chebys)/(Ncx+1))-1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !print*, 'hi - using chebys' !!print*,'TWtest0, Ncy = ',Ncy,Ncx,size(chebys),size(chebys)/(Ncx+1) xmin = Limits(1) xmax = Limits(2) tx = -1.0+2*(x-xmin)/(xmax-xmin) !! ymin = Limits(3) ymax = Limits(4) ty = -1.0+2*(y-ymin)/(ymax-ymin) !! z = 0.0 Tm0_x = 1.0 Tm1_x = tx !print*,'TW OP: x,y,tx,ty,z start',x,y,tx,ty,z ! ! do s=1,size(chebys) ! print*, 'chebys(s),s',chebys(s),s ! end do !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! Need to calculate z=\sum_{m=0}^Ncx a_m*T_m(tx), !! where a_m=\sum_{n=0}^Ncy A_{mn}*T_n(ty); !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! m=0 term;!!TWch nx = 0 Tm_x = Tm0_x am = 0.0 !! ny = 0 Tn0_y = 1.0 !s = (ny+1) + nx*(Ncx+1) s = (nx+1) + ny*(Ncx+1) Amn = chebys(s) am = am+Amn*Tn0_y !print*,'Amn,nx,ny,am,s:',Amn,nx,ny,am,s !! ny = 1 Tn1_y = ty !s = (ny+1) + nx*(Ncx+1) s = (nx+1) + ny*(Ncx+1) Amn = chebys(s) am = am+Amn*Tn1_y !print*,'Amn,nx,ny,am,s:',Amn,nx,ny,am,s !! !print*,'Ncx,Ncy,size(chebys)=',Ncx,Ncy,size(chebys) do ny=2,Ncy Tn_y = 2*ty*Tn1_y-Tn0_y Tn0_y = Tn1_y Tn1_y = Tn_y !! s = (nx+1) + ny*(Ncx+1) Amn = chebys(s) am = am+Amn*Tn1_y !print*,'Amn,nx,ny,am,s:',Amn,nx,ny,am,s end do !! z = z+am*Tm_x !print*,'TW OP: m=0 end (z,am)',z,am !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! m=1 term; nx = 1 Tm_x = Tm1_x am = 0.0 !! ny = 0 Tn0_y = 1.0 s = (nx+1) + ny*(Ncx+1) Amn = chebys(s) am = am+Amn*Tn0_y !! ny = 1 Tn1_y = ty s = (nx+1) + ny*(Ncx+1) Amn = chebys(s) am = am+Amn*Tn1_y !! do ny=2,Ncy Tn_y = 2*ty*Tn1_y-Tn0_y Tn0_y = Tn1_y Tn1_y = Tn_y !! s = (nx+1) + ny*(Ncx+1) Amn = chebys(s) am = am+Amn*Tn1_y end do !! z = z+am*Tm_x !print*,'TW OP: m=1 end (z,am)',z,am !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! Now sum over rest of terms (m=2 to m=Ncx); do nx = 2,Ncx Tm_x = 2*tx*Tm1_x-Tm0_x Tm0_x = Tm1_x Tm1_x = Tm_x am = 0.0 !! ny = 0 Tn0_y = 1.0 s = (nx+1) + ny*(Ncx+1) Amn = chebys(s) am = am+Amn*Tn0_y !! ny = 1 Tn1_y = ty s = (nx+1) + ny*(Ncx+1) Amn = chebys(s) am = am+Amn*Tn1_y !! do ny=2,Ncy Tn_y = 2*ty*Tn1_y-Tn0_y Tn0_y = Tn1_y Tn1_y = Tn_y !! s = (nx+1) + ny*(Ncx+1) Amn = chebys(s) am = am+Amn*Tn1_y end do !! z = z+am*Tm_x !print*,'TW OP: m=',nx,' end (z,am)',z,am end do !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! end subroutine OP_chebinterp2d !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! real function alpKM08(t,h) real, intent(in) :: t,h real :: p1,p2,p3 real :: q11,q12,q13,q14,q15 real :: q21,q22,q23,q24,q25 real :: q31,q32,q33,q34,q35 ! A polynomial fit of the attenuation coefficient ! Kohout and Meylan (2008) q11 = -0.00000777 q12 = 0.00032080 q13 = -0.00437542 q14 = 0.02047559 q15 = -0.01356537 q21 = 0.00003635 q22 = -0.00153484 q23 = 0.02121709 q24 = -0.09289399 q25 = -0.03693082 q31 = -0.00004509 q32 = 0.00214484 q33 = -0.03663425 q34 = 0.26065369 q35 = -0.62474085 p1 = q11*t**4 + q12*t**3 + q13*t**2 + q14*t + q15 p2 = q21*t**4 + q22*t**3 + q23*t**2 + q24*t + q25 p3 = q31*t**4 + q32*t**3 + q33*t**2 + q34*t + q35 if ( h .lt. 0.4 ) then alpKM08 = p1*0.4**2+p2*0.4+p3 else if ( h .gt. 5.0 ) then alpKM08 = p1*5.0**2+p2*5.0+p3 else alpKM08 = p1*h**2+p2*h+p3 end if if ( alpKM08 .gt. 0.0 ) alpKM08 = 0.0 end function alpKM08 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! real function alpE0(om,h) !! different fit to Kohout and Meylan (2008) ac's; real, intent(in) :: om,h integer, parameter :: Ncheb_f = 10;! deg of poly in period/freq; integer, parameter :: Ncheb_h = 10;! deg of poly in thickness; !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!SET INPUT TYPES: !X_INPUT = 1 for period; !X_INPUT = 2 for freq; !X_INPUT = 3 for omega; integer, parameter :: X_INPUT = 1 !AC_INPUT = 1 if input AC is energy AC; !AC_INPUT = 2 if input AC is amplitude AC; integer, parameter :: AC_INPUT = 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !RANGE OF VALIDITY VARIES FROM RUN TO RUN, ! SO NEED TO CHECK THIS EACH TIME; real, parameter :: hmax = 3.7 real, parameter :: hmin = 0.4 real, parameter :: xmin = 6.0 real, parameter :: xmax = 15.9 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! finally define vector of chebyshev coefficients from !! !! /home/nersc/timill/validation/miz1d/ALPfxns_23May2011/attenKM08_cheb2F.dat !! !! NB rest of function should be the same for all !! attenuation models; real, dimension((Ncheb_f+1)*(Ncheb_h+1)), parameter :: & chebys = (/ & -5.1718041e+00, & -3.0630368e+00, & 3.4726520e-02, & 7.0718546e-02, & -3.9074912e-02, & 1.5855768e-03, & 9.0309985e-03, & -3.4686153e-03, & -4.2571231e-03, & 1.5937636e-04, & 4.5062664e-04, & 2.8165857e+00, & 1.2963481e+00, & -6.4203762e-01, & -1.5116760e-01, & 9.8121338e-02, & -8.8303648e-04, & -2.0931019e-02, & 6.4713514e-03, & 5.0155575e-03, & -2.0781515e-03, & -1.6425678e-03, & -7.8534639e-01, & 1.3654805e-01, & 5.4160136e-01, & -3.4941646e-02, & -1.1651259e-01, & 1.8363695e-02, & 2.0405273e-02, & -8.1491733e-03, & -6.5704112e-03, & 3.2722509e-03, & -8.3666633e-04, & 1.9697720e-01, & -2.8087459e-01, & -2.0228810e-01, & 1.4647672e-01, & 7.6894240e-02, & -3.5861425e-02, & -1.8312850e-02, & 7.5376454e-03, & 7.7566678e-04, & -3.5595637e-03, & -2.9015505e-03, & -7.1638712e-02, & 1.2246831e-01, & -8.9356920e-03, & -1.3126669e-01, & -1.9097635e-02, & 4.5828588e-02, & 7.0505742e-03, & -1.2623172e-02, & -3.1646086e-03, & 4.1729696e-03, & 1.5960447e-03, & 4.2343603e-02, & -1.8246647e-02, & 5.0920843e-02, & 6.8929330e-02, & -2.1566054e-02, & -3.4364256e-02, & 2.7802594e-03, & 1.0290463e-02, & -1.1270861e-03, & -4.9161222e-03, & 1.0156776e-04, & -2.5390146e-02, & -1.6067386e-03, & -3.6717573e-02, & -1.7028931e-02, & 3.1235878e-02, & 1.7293652e-02, & -1.0982658e-02, & -1.0916749e-02, & 1.8212220e-03, & 3.4067499e-03, & -9.2195159e-04, & 9.1526141e-03, & -2.4790889e-03, & 1.0221778e-02, & -6.1965240e-03, & -2.1160829e-02, & -3.9377485e-03, & 1.2706661e-02, & 2.9905997e-03, & -4.0412531e-03, & -3.6382220e-03, & 9.0460708e-04, & -2.1065702e-03, & 4.7903427e-03, & -5.0119180e-04, & 7.3132784e-03, & 1.2368661e-02, & -5.6285727e-03, & -5.9370646e-03, & -1.7908780e-03, & 7.4853299e-03, & -6.3878541e-04, & -3.4269901e-03, & 5.3871224e-05, & -3.3029307e-03, & -2.8465866e-03, & -3.0070479e-03, & -3.0906007e-03, & 6.2453948e-03, & 4.4413054e-03, & -2.0603195e-03, & -3.3106400e-03, & -3.6574527e-05, & 3.6338643e-03, & 5.5765073e-04, & 3.2074228e-03, & 5.6822777e-06, & 8.9152339e-04, & -8.7498480e-04, & -4.3412112e-03, & -6.9399559e-04, & 2.7676035e-03, & 5.1223259e-03, & -2.8656280e-03, & 7.3840654e-04 & /) !end of 'chebys' matrix (E0) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! rest should be same for each program; real, parameter :: pi = 3.1415927 real :: x,y,z,zfac real, dimension(4) :: Limits zfac = 1.0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! print*,'in: om (t),h,Limits(1:4)',om,'(',2*pi/om,')',h,Limits(1), ! & Limits(2),Limits(3),Limits(4) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! check thickness is OK y = h if (h.gt.hmax) then !print*, 'TW (E0): h > ',hmax, '; h = ',h y = hmax ! call xcstop('attenE0') elseif (h.lt.hmin) then !use linear interpolation for small h zfac = h/hmin y = hmin end if !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! check 'x' is OK if (X_INPUT.eq.1) then x = (2*pi)/om !xstr = 'period'; elseif (X_INPUT.eq.2) then x = om/(2*pi) !xstr = 'freq'; else x = om !xstr = 'omega'; end if if (x.gt.xmax) then print*, 'x > ',xmax, '; x = ',x call xcstop('attenE0') elseif (x.lt.xmin) then print*, 'x < ',xmin, '; x = ',x call xcstop('attenE0') end if !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! Do interpolation Limits(1) = xmin Limits(2) = xmax Limits(3) = hmin Limits(4) = hmax call OP_chebinterp2d(z,x,y,chebys,Ncheb_f,Limits) !print*,'out from OP_chebinterp2d',exp(z) !! want to OUTPUT ENERGY attenuation coefficients: alpE0 = zfac*AC_INPUT*exp(z) end function!! alpE0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! real function alpE1(om,h) !! different fit to Kohout and Meylan (2008) ac's; real, intent(in) :: om,h integer, parameter :: Ncheb_f = 25;! deg of poly in period/freq; integer, parameter :: Ncheb_h = 25;! deg of poly in thickness; !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!SET INPUT TYPES: !X_INPUT = 1 for period; !X_INPUT = 2 for freq; !X_INPUT = 3 for omega; integer, parameter :: X_INPUT = 1 !AC_INPUT = 1 if input AC is energy AC; !AC_INPUT = 2 if input AC is amplitude AC; integer, parameter :: AC_INPUT = 2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !RANGE OF VALIDITY VARIES FROM RUN TO RUN, ! SO NEED TO CHECK THIS EACH TIME; real, parameter :: hmax = 5.0 real, parameter :: hmin = 0.1 real, parameter :: xmin = 2.5 real, parameter :: xmax = 25.0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! finally define vector of chebyshev coefficients from !! /home/nersc/timill/validation/miz1d/ALPfxns_23May2011/... !! attenTWrslPL_c90_cheb2.dat !! !! NB rest of function should be the same for all !! attenuation models; real, dimension((Ncheb_f+1)*(Ncheb_h+1)), parameter :: & chebys = (/ & -8.3662374e+00, & -8.4283985e+00, & 1.1547019e+00, & -4.1302974e-02, & -4.6868249e-01, & 3.0865993e-01, & 3.5716480e-01, & -1.5830466e-01, & -3.4793682e-03, & 9.0858199e-02, & 9.3333584e-02, & -1.2837279e-03, & -6.1739822e-02, & 4.5754541e-03, & 1.0306449e-02, & -2.0909625e-02, & -8.4736775e-02, & -3.2592847e-02, & 3.9062508e-02, & -2.6843630e-02, & -5.7361742e-02, & -1.8844557e-02, & -7.5143369e-03, & -2.4737010e-02, & -6.0097639e-02, & -1.4756351e-02, & 5.0474962e+00, & 2.3397436e+00, & -1.0879456e+00, & 8.1328415e-02, & -4.9729252e-01, & -1.3398148e-01, & 6.7128370e-01, & -2.2984616e-01, & -3.6994226e-01, & 2.4986688e-01, & 7.7288409e-02, & -1.1547666e-02, & -2.0146949e-02, & 8.8262905e-02, & 2.1135622e-01, & 7.3836412e-02, & -1.5189678e-02, & -3.7433776e-02, & 4.2158095e-02, & 7.0595194e-02, & -5.8215704e-02, & -4.8172618e-02, & 1.2159463e-02, & -3.8781631e-02, & -8.0691749e-02, & -7.4961639e-02, & -1.5045246e+00, & 4.8207369e-01, & 1.5365525e+00, & 2.4078323e-01, & -3.4106031e-02, & -3.1111748e-01, & 3.4249780e-01, & 4.4874301e-01, & -2.1539841e-01, & -1.6657296e-02, & 2.8190769e-01, & 3.8020365e-02, & -1.0222629e-01, & -3.8165577e-02, & 1.8874930e-02, & -8.5309789e-03, & -1.4576568e-01, & -9.5605150e-02, & 6.1217192e-02, & -3.5557414e-03, & -9.7045359e-02, & -5.8852987e-02, & -5.1434831e-03, & -3.2811288e-02, & -1.0896159e-01, & -3.8648807e-02, & 4.1107732e-01, & -7.4207130e-01, & -5.9303000e-01, & 5.5677778e-01, & 4.1362917e-01, & 6.1494115e-02, & -3.2242870e-01, & 2.0501736e-01, & 9.4800168e-02, & -2.1413823e-01, & 4.3624007e-02, & 1.5185802e-01, & -1.5549111e-03, & -3.3690306e-03, & 1.3401742e-01, & 1.0926770e-01, & 3.1697681e-02, & -9.3807331e-02, & -2.8012301e-02, & 1.0444106e-01, & -2.6535139e-02, & -1.1692068e-01, & -4.1271554e-03, & -1.6866227e-02, & -7.2037930e-02, & -1.1047555e-01, & -2.9783392e-01, & 3.1195054e-01, & 8.4378397e-02, & -4.3343614e-01, & -5.4787079e-02, & 2.6327352e-01, & 1.8031827e-01, & -8.6567669e-02, & 3.0163554e-01, & 1.1597234e-01, & -1.3791765e-01, & 1.3439102e-01, & 9.7980663e-02, & -9.2247698e-02, & -7.9495904e-02, & 1.9761924e-02, & -4.2512683e-02, & -9.7350236e-02, & -4.0551686e-02, & 1.2688252e-02, & -4.7106294e-03, & -1.1962453e-01, & -4.5120196e-02, & -8.3208269e-03, & -6.3262915e-02, & -7.3963657e-02, & 1.5447957e-01, & -7.9851518e-02, & 1.6119479e-01, & 4.1563374e-01, & -2.9704458e-02, & -1.3606532e-01, & -2.6877353e-02, & 1.7021096e-02, & -2.1697938e-01, & 1.9017393e-01, & 2.8134184e-02, & -1.6412770e-01, & 1.2885736e-01, & 1.1949535e-01, & -3.0452530e-02, & 3.0870338e-02, & 1.1975121e-01, & -1.4382364e-02, & -7.1896075e-02, & 3.6959715e-02, & 9.7484398e-02, & -6.3277663e-02, & -6.1365171e-02, & 8.2169441e-03, & 2.5806329e-02, & -9.4040327e-02, & -9.0657930e-02, & 1.5459918e-01, & 2.1876872e-02, & 2.2522247e-02, & 2.4108091e-01, & 2.2847771e-01, & 2.7557426e-02, & 1.0187506e-01, & 1.2132767e-01, & -1.1036645e-01, & 1.3308192e-01, & 8.0193357e-02, & -1.3676299e-01, & 8.2196723e-02, & 2.3510945e-02, & -1.2141176e-01, & -2.6926939e-02, & 2.2873066e-02, & -2.7567170e-02, & -8.0294057e-02, & 2.5951138e-02, & 3.0401004e-03, & -5.8456963e-02, & -6.2835393e-02, & 3.9446734e-02, & -1.8691166e-02, & -5.8332879e-02, & -2.4971185e-01, & -8.1437311e-02, & -3.5424987e-02, & -1.3025558e-01, & -3.4154958e-02, & 6.9511989e-02, & -4.0699089e-02, & -1.3077449e-02, & 1.0075631e-01, & -1.2673466e-01, & 3.1302889e-02, & 7.3159687e-02, & -6.6350015e-02, & 1.1639126e-01, & 6.1613518e-02, & -4.5316588e-02, & 3.2041604e-02, & 3.5863206e-02, & 4.6350967e-03, & 3.3830327e-03, & 4.1919659e-02, & 9.6396743e-03, & -4.2271025e-02, & -3.1834703e-05, & 3.1804791e-02, & 3.0853972e-02, & 1.3361369e-01, & 1.0485792e-01, & 4.3647094e-02, & 9.7031506e-02, & -3.1110281e-02, & -4.1355715e-02, & 3.6461492e-02, & 3.1349584e-02, & -4.9575844e-03, & 1.3587262e-01, & -5.8650824e-02, & 2.3382289e-02, & 4.4591732e-02, & -1.0976292e-01, & 3.5433224e-02, & -8.7393078e-03, & -6.9981435e-02, & 1.8166705e-02, & 3.1027707e-02, & -6.2915342e-02, & -1.5031376e-02, & 2.9628069e-03, & -1.0250745e-02, & -8.2729756e-02, & 3.1539708e-02, & 2.2979664e-03, & -3.8914987e-02, & 4.2666482e-02, & 6.5528941e-02, & 5.2830220e-02, & 9.6192979e-02, & 6.6730528e-02, & -4.5647645e-02, & -6.5303784e-03, & -7.3114528e-03, & -6.0532942e-02, & 6.1605942e-02, & -4.1783026e-02, & 2.1741581e-02, & 7.4860205e-02, & -7.2091286e-02, & 4.2323299e-02, & 1.0177118e-02, & -7.3686153e-02, & 2.8015712e-02, & 1.7386535e-02, & -5.1492911e-02, & -1.4791293e-02, & 2.6605068e-02, & -3.3709689e-02, & -5.4208702e-02, & -3.6042074e-02, & 1.4726675e-02, & -4.8883506e-02, & -7.6219522e-03, & -2.9197810e-02, & -8.7894514e-03, & -2.9342315e-03, & 9.2255681e-02, & 8.1981077e-03, & 3.7655556e-02, & 2.2177272e-02, & -4.5125967e-02, & 4.4951508e-02, & -3.5161952e-02, & -1.9904095e-02, & 4.7882545e-02, & -6.4732014e-02, & 2.4286707e-02, & 3.5892766e-02, & -6.5570349e-02, & 3.5969425e-02, & -5.3922386e-03, & -3.3691910e-02, & -2.8713714e-02, & 3.9813703e-02, & -3.9327678e-02, & 2.2868849e-02, & -1.6979304e-02, & -7.9415373e-03, & -4.5469071e-02, & -4.3484922e-02, & -3.4397431e-02, & -3.0870933e-02, & -6.7404448e-02, & 2.1750796e-03, & -2.1716505e-02, & 1.5163756e-03, & 3.6463258e-02, & -4.7465468e-02, & 2.9985599e-02, & -1.9342078e-02, & -3.1915157e-02, & 4.8991598e-02, & -4.6984350e-02, & 2.8433052e-02, & 4.0321621e-02, & -3.9066749e-02, & 5.1322374e-02, & 3.5104624e-02, & -4.5878770e-03, & 1.8800197e-02, & 4.7823862e-02, & -1.9226949e-02, & -2.7780977e-02, & -2.6033307e-02, & -2.5285944e-02, & 1.7260343e-03, & -3.5208010e-02, & -2.7749655e-02, & -2.7751790e-02, & -5.9696034e-02, & -3.5348632e-02, & -3.0419618e-02, & -3.3913219e-02, & 3.9637203e-02, & -3.4882566e-02, & 2.5036799e-02, & 2.1176814e-02, & -1.2474542e-02, & 6.8949473e-02, & -4.3619761e-02, & 2.0048293e-02, & 3.5151879e-02, & -2.3778639e-02, & 3.0427746e-02, & 4.1205130e-02, & -1.3245610e-02, & 2.6488226e-02, & 7.1575690e-02, & 7.0548984e-02, & 1.7655842e-02, & -3.2959130e-02, & -2.9965924e-02, & 1.0205483e-02, & 4.6327923e-02, & 4.0129336e-02, & 5.3580779e-02, & -1.4249633e-02, & -4.5859669e-03, & -4.6699255e-03, & -5.4185202e-02, & -4.3143417e-03, & -5.2149602e-02, & -9.5354611e-03, & 2.9551683e-04, & -3.1497195e-02, & 5.4986544e-02, & -4.6388507e-02, & -1.7524007e-02, & 4.8940157e-02, & -3.4625932e-02, & 2.0584076e-02, & 1.7791421e-02, & -1.0100843e-02, & -7.4014811e-02, & -9.8372136e-02, & -1.3190285e-01, & -7.7061390e-02, & -8.1634074e-02, & -7.9788427e-02, & -1.1250340e-01, & -6.6883799e-02, & -7.3448485e-02, & -2.6625503e-04, & -2.4705332e-02, & -1.1484820e-02, & 1.3225935e-02, & -2.4135377e-02, & 2.8118773e-02, & 1.1644702e-02, & 3.9816768e-02, & 4.0471970e-02, & -1.9247900e-02, & 4.9717874e-02, & 1.5621069e-02, & 1.0241904e-03, & 8.2124060e-02, & -1.2152046e-02, & 6.3184950e-02, & 2.0356901e-02, & 1.2269586e-01, & 1.0579580e-01, & 6.0980302e-02, & -6.1197782e-02, & -2.6812461e-02, & -3.3974209e-02, & 4.6040776e-02, & 4.8078166e-03, & 2.3923645e-02, & -6.8352227e-02, & -7.5404329e-03, & -9.0402328e-03, & -3.0293465e-02, & -2.3236227e-02, & -5.6620217e-02, & -1.4309425e-02, & -3.9815389e-02, & 3.0352905e-03, & 3.2757432e-02, & -3.5375148e-02, & 2.3308640e-03, & 4.9533541e-02, & -1.7507751e-02, & 7.3281920e-02, & -3.0582219e-02, & 6.2040212e-02, & -6.8112776e-02, & -8.5671579e-02, & -1.2807058e-01, & -6.1582850e-02, & -6.4803912e-02, & -5.6854078e-02, & -1.0008646e-01, & -5.5177538e-02, & -6.5414937e-02, & -9.8924527e-03, & -4.9796517e-02, & -3.1160333e-02, & -9.9657718e-03, & -1.5630335e-02, & 2.7980976e-02, & 2.8749503e-02, & 6.1648919e-02, & 2.3568303e-02, & 9.6647382e-04, & 3.5334201e-02, & 8.0508379e-04, & 1.9916542e-02, & 5.1074844e-02, & -8.0644913e-03, & 7.8470689e-02, & -2.5225654e-02, & 1.0626573e-01, & 9.4698229e-02, & 1.4802024e-02, & -8.2483809e-02, & -6.6546354e-02, & -3.9574944e-02, & 3.5608091e-02, & 2.7996424e-02, & 3.2371695e-02, & -3.0096504e-02, & 1.6964380e-02, & 1.0482192e-02, & -3.5077555e-02, & -5.8997584e-02, & -7.9031878e-02, & -4.1226364e-02, & -4.1781096e-02, & 1.8255380e-02, & 3.4377537e-02, & -2.3051899e-02, & 1.8271517e-02, & 4.2343973e-02, & 2.8900726e-02, & 4.3526519e-02, & -2.8492140e-03, & 7.6669064e-02, & -5.0719717e-02, & -9.6864483e-02, & -1.0468344e-01, & -7.8881807e-02, & -4.6821010e-02, & -8.2245285e-02, & -1.0519775e-01, & -9.4697104e-02, & -7.3768228e-02, & -4.6653637e-02, & -5.2702559e-02, & -3.0478130e-02, & 1.1951399e-02, & 1.5297283e-02, & 3.6948352e-02, & 4.0434467e-02, & 4.8094992e-02, & 1.6762028e-02, & -4.5035668e-03, & 3.0261526e-02, & -6.4188478e-03, & 3.2016601e-02, & 1.8277364e-02, & 3.4227518e-02, & 4.6431580e-02, & -1.6428566e-03, & 9.8291679e-02, & 1.2006655e-01, & 1.9994809e-02, & -4.4050081e-02, & -5.7657296e-02, & -8.1750688e-03, & 3.5316735e-02, & 4.0881073e-02, & 1.7413434e-02, & -2.2497754e-02, & -1.6221592e-03, & -3.8310122e-03, & -5.3696353e-02, & -7.1427916e-02, & -6.9875429e-02, & -4.0622409e-02, & -2.9862924e-02, & 1.2238246e-02, & 2.2256397e-02, & -4.1708486e-02, & -8.0707704e-03, & 1.1296182e-02, & 2.4697406e-02, & -3.8270115e-04, & 2.1537978e-02, & 4.8280080e-02, & -3.1228346e-02, & -7.6721454e-02, & -6.8062074e-02, & -6.1835623e-02, & -2.3588901e-02, & -6.8876651e-02, & -7.5243327e-02, & -7.6979433e-02, & -4.4841723e-02, & -2.8883329e-02, & -2.4028196e-02, & -1.9570395e-02, & 1.3384785e-02, & 1.3958883e-02, & 2.1380503e-02, & 3.2816559e-02, & 3.9510682e-02, & 2.4550858e-02, & 3.5939357e-03, & 3.1444848e-02, & 6.4498869e-03, & 3.6898350e-02, & 1.4882664e-02, & 4.5062865e-02, & 2.1975905e-02, & 1.5843098e-02, & 7.1100027e-02, & 8.7051847e-02, & -6.0775324e-03, & -5.7077944e-02, & -7.5333716e-02, & -3.0275680e-02, & 2.5516642e-03, & 1.9168757e-02, & 6.7108744e-04, & -2.4971514e-02, & -1.2734005e-02, & -5.0344718e-03, & -4.0562667e-02, & -5.3715513e-02, & -4.6952160e-02, & -3.3673116e-02, & -2.8865063e-02, & 1.2384125e-03, & 1.7144502e-02, & -2.8313472e-02, & -9.3643438e-03, & 1.2573054e-02, & 2.4217606e-02, & -5.3865370e-04, & 2.9797448e-02, & 2.8356832e-02, & -2.0506516e-02, & -4.7731417e-02, & -3.4999854e-02, & -2.5778041e-02, & -2.2455897e-03, & -3.6922315e-02, & -4.8118709e-02, & -5.4602488e-02, & -3.8994207e-02, & -3.3257505e-02, & -2.7955448e-02, & -2.0479142e-02, & 7.6154495e-03, & 1.5340374e-02, & 2.2805989e-02, & 3.7772992e-02, & 4.1243450e-02, & 2.1724053e-02, & -5.4016676e-03, & 3.6773038e-03, & -4.5264207e-03, & 1.2544675e-02, & 5.3858605e-03, & 2.6155512e-02, & 1.0562263e-02, & 1.2628957e-02, & 4.5604837e-02, & 5.6055102e-02, & -6.5303928e-03, & -4.1115891e-02, & -4.9994022e-02, & -1.7144661e-02, & 1.0812350e-02, & 2.5094927e-02, & 1.3688704e-02, & -3.3741309e-03, & -1.9387442e-03, & -4.3749156e-03, & -3.1987984e-02, & -4.4711413e-02, & -4.0193037e-02, & -3.1367765e-02, & -2.3035258e-02, & 4.7467753e-03, & 2.3083008e-02, & -5.9627897e-03, & -4.4952200e-03, & 1.3234452e-02, & 1.7352332e-02, & -3.8695632e-04, & 1.5577596e-02, & 1.2077597e-02, & -1.4231916e-02, & -3.4668590e-02, & -2.7940883e-02, & -1.9369642e-02, & -8.2111733e-03, & -2.9262485e-02, & -3.6617467e-02, & -3.7246184e-02, & -2.5651486e-02, & -2.0867098e-02, & -1.2706216e-02, & -3.7278751e-03, & 1.1656245e-02, & 1.5306222e-02, & 1.5298721e-02, & 1.8333182e-02, & 1.6831384e-02, & 6.9685831e-03, & -6.1861895e-03, & 8.5097910e-04, & 1.3493193e-03, & 6.1964217e-03, & 2.3286579e-03, & 1.3622444e-02, & 7.2949065e-03, & 9.7162175e-03, & 1.8660258e-02, & 2.5053612e-02, & -2.5429952e-03, & -1.7966794e-02, & -2.1766543e-02, & -6.6149383e-03, & 2.8611476e-03, & 6.4771225e-03, & -2.2176426e-03, & -9.0823885e-03, & -8.2409964e-03, & -6.5534596e-03, & -1.4879434e-02, & -1.5975703e-02, & -8.5618582e-03, & -2.1155861e-03, & -5.1498946e-04, & 5.3685867e-03, & 8.0029692e-03, & -9.1526193e-03, & -8.6027026e-03, & 1.6159462e-03, & 5.8102390e-03, & 1.8590681e-03, & 1.1555711e-02, & 4.3003470e-03 & /) !end of 'chebys' matrix (E1) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! rest should be same for each program; real, parameter :: pi = 3.1415927 real :: x,y,z,zfac real, dimension(4) :: Limits zfac = 1.0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! print*,'in: om (t),h,Limits(1:4)',om,'(',2*pi/om,')',h,Limits(1), ! & Limits(2),Limits(3),Limits(4) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! check thickness is OK y = h if (h.gt.hmax) then !print*, 'TW (E1): h > ',hmax, '; h = ',h y = hmax ! call xcstop('attenE1') elseif (h.lt.hmin) then !use linear interpolation for small h zfac = h/hmin y = hmin end if !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! check 'x' is OK if (X_INPUT.eq.1) then x = (2*pi)/om !xstr = 'period'; elseif (X_INPUT.eq.2) then x = om/(2*pi) !xstr = 'freq'; else x = om !xstr = 'omega'; end if if (x.gt.xmax) then print*, 'x > ',xmax, '; x = ',x call xcstop('attenE1') elseif (x.lt.xmin) then print*, 'x < ',xmin, '; x = ',x call xcstop('attenE1') end if !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! Do interpolation Limits(1) = xmin Limits(2) = xmax Limits(3) = hmin Limits(4) = hmax call OP_chebinterp2d(z,x,y,chebys,Ncheb_f,Limits) !print*,'out from OP_chebinterp2d',exp(z) !! want to OUTPUT ENERGY attenuation coefficients: alpE1 = zfac*AC_INPUT*exp(z) end function!! alpE1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! real function alpA0(om,h) !! -2*log(1-|R|^2), R is semi-infinite ref coeff; !! no draft real, intent(in) :: om,h integer, parameter :: Ncheb_f = 10;! deg of poly in period/freq; integer, parameter :: Ncheb_h = 10;! deg of poly in thickness; !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!SET INPUT TYPES: !X_INPUT = 1 for period; !X_INPUT = 2 for freq; !X_INPUT = 3 for omega; integer, parameter :: X_INPUT = 1 !AC_INPUT = 1 if input AC is energy AC; !AC_INPUT = 2 if input AC is amplitude AC; integer, parameter :: AC_INPUT = 2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !RANGE OF VALIDITY VARIES FROM RUN TO RUN, ! SO NEED TO CHECK THIS EACH TIME; real, parameter :: hmax = 5.0 real, parameter :: hmin = 0.1 real, parameter :: xmin = 2.5 real, parameter :: xmax = 25.0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! finally define vector of chebyshev coefficients from !! /home/nersc/timill/validation/miz1d/ALPfxns_23May2011/... !! attenTWrslPL_c90_cheb2.dat !! !! NB rest of function should be the same for all !! attenuation models; real, dimension((Ncheb_f+1)*(Ncheb_h+1)), parameter :: & chebys = (/ & -6.3648136e+00, & -5.9408354e+00, & 6.0823789e-01, & 1.5169419e-01, & 2.4826164e-02, & 2.0454896e-02, & -3.6035906e-02, & 9.0007967e-03, & -6.4465368e-04, & -3.4958570e-03, & 4.4870990e-03, & 3.4685912e+00, & 9.6252301e-01, & -1.6341767e+00, & 3.6501777e-01, & 2.5319557e-01, & -7.9817426e-02, & 2.8837743e-02, & -3.5071614e-02, & 2.3583588e-03, & 6.7852646e-03, & -5.9944778e-03, & -1.3218570e+00, & 8.2277092e-02, & 4.9041759e-01, & -6.4724337e-01, & 6.4169237e-02, & 2.0080448e-01, & -6.8595064e-02, & 9.1966953e-03, & -7.2502222e-03, & -6.6158350e-03, & 9.4112277e-03, & 6.9367788e-01, & -8.6700675e-02, & -8.8653441e-02, & 3.4359908e-01, & -2.7490483e-01, & -6.7983152e-02, & 1.4720605e-01, & -3.2016416e-02, & -9.4066399e-03, & 2.4269262e-03, & -6.7678627e-03, & -3.9655814e-01, & 7.8473574e-02, & 1.3738607e-02, & -1.3608019e-01, & 2.2028195e-01, & -7.8832728e-02, & -1.0894200e-01, & 9.1473493e-02, & 7.2344011e-04, & -1.6653275e-02, & 6.6437295e-03, & 2.4312181e-01, & -6.0533708e-02, & 9.4380333e-03, & 5.6954281e-02, & -1.2661095e-01, & 1.1121865e-01, & 2.0709289e-02, & -9.7523323e-02, & 3.9686739e-02, & 1.7773563e-02, & -1.5782650e-02, & -1.5433188e-01, & 4.3040435e-02, & -1.5611582e-02, & -2.3966901e-02, & 6.9329410e-02, & -8.6387080e-02, & 2.8620852e-02, & 5.9793960e-02, & -6.3151360e-02, & 5.7950282e-03, & 2.5014762e-02, & 9.9977637e-02, & -3.0395925e-02, & 1.4665657e-02, & 8.6048269e-03, & -3.8311806e-02, & 5.7398803e-02, & -3.9220324e-02, & -2.2656700e-02, & 5.6567836e-02, & -3.2533368e-02, & -1.7155952e-02, & -6.5743136e-02, & 2.1360701e-02, & -1.1885306e-02, & -1.8381123e-03, & 2.1113921e-02, & -3.6572033e-02, & 3.4064767e-02, & 1.6692253e-03, & -3.7783197e-02, & 4.2096749e-02, & -1.3615713e-03, & 3.9370482e-02, & -1.3368301e-02, & 7.9838551e-03, & -4.8427528e-04, & -1.0790830e-02, & 2.1044687e-02, & -2.2942918e-02, & 5.3101289e-03, & 2.0191020e-02, & -3.2755924e-02, & 1.1804324e-02, & -2.6078045e-02, & 9.2128632e-03, & -5.8379292e-03, & 1.2977307e-03, & 5.8160425e-03, & -1.3294198e-02, & 1.6230862e-02, & -8.2406917e-03, & -9.8938125e-03, & 2.2804706e-02, & -1.7486848e-02 & /) !end of 'chebys' matrix (A0) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! rest should be same for each program; real, parameter :: pi = 3.1415927 real :: x,y,z,zfac real, dimension(4) :: Limits zfac = 1.0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! print*,'in: om (t),h,Limits(1:4)',om,'(',2*pi/om,')',h,Limits(1), ! & Limits(2),Limits(3),Limits(4) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! check thickness is OK y = h if (h.gt.hmax) then !print*, 'TW (A0): h > ',hmax, '; h = ',h y = hmax ! call xcstop('attenA0') elseif (h.lt.hmin) then !use linear interpolation for small h zfac = h/hmin y = hmin end if !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! check 'x' is OK if (X_INPUT.eq.1) then x = (2*pi)/om !xstr = 'period'; elseif (X_INPUT.eq.2) then x = om/(2*pi) !xstr = 'freq'; else x = om !xstr = 'omega'; end if if (x.gt.xmax) then print*, 'x > ',xmax, '; x = ',x call xcstop('attenA0') elseif (x.lt.xmin) then print*, 'x < ',xmin, '; x = ',x call xcstop('attenA0') end if !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! Do interpolation Limits(1) = xmin Limits(2) = xmax Limits(3) = hmin Limits(4) = hmax call OP_chebinterp2d(z,x,y,chebys,Ncheb_f,Limits) !print*,'out from OP_chebinterp2d',exp(z) !! want to OUTPUT ENERGY attenuation coefficients: alpA0 = zfac*AC_INPUT*exp(z) end function!! alpA0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! real function alpA1(om,h) !! -2*log(1-|R|^2), R is semi-infinite ref coeff; !! has draft real, intent(in) :: om,h integer, parameter :: Ncheb_f = 10;! deg of poly in period/freq; integer, parameter :: Ncheb_h = 10;! deg of poly in thickness; !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!SET INPUT TYPES: !X_INPUT = 1 for period; !X_INPUT = 2 for freq; !X_INPUT = 3 for omega; integer, parameter :: X_INPUT = 2 !AC_INPUT = 1 if input AC is energy AC; !AC_INPUT = 2 if input AC is amplitude AC; integer, parameter :: AC_INPUT = 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !RANGE OF VALIDITY VARIES FROM RUN TO RUN, ! SO NEED TO CHECK THIS EACH TIME; real, parameter :: hmax = 5.0 real, parameter :: hmin = 0.1 real, parameter :: xmin = 0.042 real, parameter :: xmax = 0.4 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! finally define vector of chebyshev coefficients from !! /home/nersc/timill/validation/miz1d/ALPfxns_23May2011/... !! ALPfxn_SAdamp_SUB_cheb2.m !! !! NB rest of function should be the same for all !! attenuation models; real, dimension((Ncheb_f+1)*(Ncheb_h+1)), parameter :: & chebys = (/ & -1.8501579e+00, & 4.9101877e+00, & -1.9944048e+00, & 7.3571385e-01, & -2.2388319e-01, & 4.3700266e-02, & 1.4351128e-02, & -2.5879087e-02, & 1.9640453e-02, & -1.0363481e-02, & 3.2663770e-03, & 2.7695669e+00, & -1.2063711e+00, & -5.0784397e-02, & 3.7480698e-01, & -2.6425166e-01, & 9.1308132e-02, & 6.9547408e-03, & -4.3420410e-02, & 4.6092399e-02, & -3.1540644e-02, & 1.5837220e-02, & -1.2036662e+00, & 3.3680013e-01, & 1.9808329e-01, & -2.3814800e-01, & 8.0319397e-02, & 3.7840528e-02, & -6.0775037e-02, & 3.5031731e-02, & -6.0416843e-03, & -8.9097036e-03, & 1.5726785e-02, & 6.4898367e-01, & -1.2109648e-01, & -1.5806140e-01, & 1.3688927e-01, & -3.3480367e-03, & -6.1277527e-02, & 4.1973575e-02, & -5.6549417e-03, & -1.1632590e-02, & 1.2089542e-02, & -6.4837898e-03, & -3.7936273e-01, & 5.3998050e-02, & 1.1371674e-01, & -7.5373235e-02, & -2.6600942e-02, & 5.5102220e-02, & -2.1563646e-02, & -8.0610915e-03, & 1.2680156e-02, & -5.9534129e-03, & -1.1215486e-03, & 2.3485146e-01, & -2.8291759e-02, & -8.0908767e-02, & 3.8793499e-02, & 3.4616020e-02, & -3.9923118e-02, & 4.7538859e-03, & 1.4474075e-02, & -1.0138135e-02, & 1.0790529e-03, & 3.2584851e-03, & -1.5053499e-01, & 1.2525344e-02, & 5.7500633e-02, & -1.8200037e-02, & -3.1466412e-02, & 2.5868808e-02, & 4.0116291e-03, & -1.4768318e-02, & 6.5896646e-03, & 1.8758681e-03, & -3.7027599e-03, & 1.0190258e-01, & -2.2429425e-04, & -4.0608910e-02, & 6.9273313e-03, & 2.5656968e-02, & -1.5656594e-02, & -7.7986690e-03, & 1.2440697e-02, & -3.0118289e-03, & -3.7434000e-03, & 3.4392031e-03, & -6.7257320e-02, & -4.5018646e-03, & 2.6230343e-02, & -1.4532063e-03, & -1.8832335e-02, & 8.5065359e-03, & 8.4602431e-03, & -9.2579927e-03, & 4.6304848e-04, & 4.2940931e-03, & -2.7030121e-03, & 3.8848120e-02, & 3.8380170e-03, & -1.4649663e-02, & -2.4250313e-04, & 1.1576818e-02, & -3.9017295e-03, & -6.5477149e-03, & 5.7063587e-03, & 7.0878167e-04, & -3.4892451e-03, & 1.7178745e-03, & -2.3909477e-02, & -2.8499561e-03, & 8.8697690e-03, & 5.1143991e-04, & -7.5240865e-03, & 1.6398133e-03, & 5.0596886e-03, & -3.4272682e-03, & -1.3717067e-03, & 2.8594988e-03, & -9.5172114e-04 & /) !end of 'chebys' matrix (A0) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! rest should be same for each program; real, parameter :: pi = 3.1415927 real :: x,y,z,zfac real, dimension(4) :: Limits zfac = 1.0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! print*,'in: om (t),h,Limits(1:4)',om,'(',2*pi/om,')',h,Limits(1), ! & Limits(2),Limits(3),Limits(4) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! check thickness is OK y = h if (h.gt.hmax) then !print*, 'TW (A1): h > ',hmax, '; h = ',h y = hmax ! call xcstop('attenA1') elseif (h.lt.hmin) then !use linear interpolation for small h zfac = h/hmin y = hmin end if !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! check 'x' is OK if (X_INPUT.eq.1) then x = (2*pi)/om !xstr = 'period'; elseif (X_INPUT.eq.2) then x = om/(2*pi) !xstr = 'freq'; else x = om !xstr = 'omega'; end if if (x.gt.xmax) then print*, 'x > ',xmax, '; x = ',x call xcstop('attenA1') elseif (x.lt.xmin) then print*, 'x < ',xmin, '; x = ',x call xcstop('attenA1') end if !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! Do interpolation Limits(1) = xmin Limits(2) = xmax Limits(3) = hmin Limits(4) = hmax call OP_chebinterp2d(z,x,y,chebys,Ncheb_f,Limits) !print*,'out from OP_chebinterp2d',exp(z) !! want to OUTPUT ENERGY attenuation coefficients: !! interpolated log(alpha) to make sure alpha>0 alpA1 = zfac*AC_INPUT*exp(z) end function!! alpA1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! real function dampingAD1(om,h) !! -2*log(1-|R|^2), R is semi-infinite ref coeff; !! no draft real, intent(in) :: om,h integer, parameter :: Ncheb_f = 10;! deg of poly in period/freq; integer, parameter :: Ncheb_h = 10;! deg of poly in thickness; !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!SET INPUT TYPES: !X_INPUT = 1 for period; !X_INPUT = 2 for freq; !X_INPUT = 3 for omega; integer, parameter :: X_INPUT = 2 !AC_INPUT = 1 if input AC is energy AC; !AC_INPUT = 2 if input AC is amplitude AC; integer, parameter :: AC_INPUT = 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !RANGE OF VALIDITY VARIES FROM RUN TO RUN, ! SO NEED TO CHECK THIS EACH TIME; real, parameter :: hmax = 5.0 real, parameter :: hmin = 0.1 real, parameter :: xmin = 0.042 real, parameter :: xmax = 0.4 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! finally define vector of chebyshev coefficients from !! /home/nersc/timill/validation/miz1d/ALPfxns_23May2011/... !! ALPfxn_SAdamp_SUB_cheb2.m !! !! NB rest of function should be the same for all !! attenuation models; real, dimension((Ncheb_f+1)*(Ncheb_h+1)), parameter :: & chebys = (/ & -1.2233019e+01, & 1.4519969e-01, & -3.0148111e-01, & 2.3742724e-01, & -1.5423088e-01, & 9.6480359e-02, & -5.9199747e-02, & 3.4991675e-02, & -1.9572012e-02, & 9.8009109e-03, & -4.2839433e-03, & -1.5073001e+00, & -7.3553740e-01, & 5.2598861e-01, & -2.5163648e-01, & 7.1908038e-02, & 6.2572098e-03, & -2.8200328e-02, & 2.8164817e-02, & -2.0418112e-02, & 1.1440691e-02, & -6.1785019e-03, & 4.9484433e-01, & 4.0946407e-01, & -2.2304337e-01, & 3.0487214e-02, & 5.7320175e-02, & -5.9445032e-02, & 3.4480403e-02, & -1.1908350e-02, & -1.3558438e-03, & 5.6946671e-03, & -6.1196024e-03, & -2.3318404e-01, & -2.3101814e-01, & 8.7637783e-02, & 3.3088361e-02, & -5.9667711e-02, & 3.1089041e-02, & -1.3657105e-03, & -1.1877202e-02, & 1.2965055e-02, & -8.8384312e-03, & 4.1215163e-03, & 1.2414779e-01, & 1.3891949e-01, & -3.5294341e-02, & -3.8149307e-02, & 3.6841709e-02, & -5.7353116e-03, & -1.3700209e-02, & 1.4190742e-02, & -6.9953807e-03, & 1.2273073e-03, & 2.4705757e-03, & -7.1257123e-02, & -8.8753470e-02, & 1.3123439e-02, & 3.0938423e-02, & -1.8761548e-02, & -6.3752909e-03, & 1.4583399e-02, & -7.6341361e-03, & -7.1560475e-04, & 3.7549782e-03, & -3.8705399e-03, & 4.2180028e-02, & 5.7105361e-02, & -3.5010876e-03, & -2.2144227e-02, & 8.3109209e-03, & 9.4134848e-03, & -1.0648099e-02, & 1.9406586e-03, & 4.2310236e-03, & -4.4944972e-03, & 2.2127542e-03, & -2.5778203e-02, & -3.7181215e-02, & -7.3173365e-04, & 1.5259532e-02, & -2.8733395e-03, & -8.8122776e-03, & 6.4546850e-03, & 1.3593037e-03, & -4.7253396e-03, & 3.1414945e-03, & -8.4790829e-05, & 1.5835451e-02, & 2.3867144e-02, & 2.2131525e-03, & -1.0122262e-02, & 3.6525172e-04, & 6.8649923e-03, & -3.3519888e-03, & -2.5467776e-03, & 3.7813389e-03, & -1.5179202e-03, & -1.1687535e-03, & -8.8847762e-03, & -1.3757294e-02, & -2.0636176e-03, & 5.9459058e-03, & 4.5093754e-04, & -4.3833673e-03, & 1.4296488e-03, & 2.2285197e-03, & -2.3517064e-03, & 4.3926964e-04, & 1.3171873e-03, & 5.3308274e-03, & 8.4306204e-03, & 1.7888482e-03, & -3.6963192e-03, & -7.1646906e-04, & 2.8757952e-03, & -4.4809105e-04, & -1.8383124e-03, & 1.4040431e-03, & 1.7212564e-04, & -1.2364047e-03 & /) !end of 'chebys' matrix (AD1) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! rest should be same for each program; real, parameter :: pi = 3.1415927 real :: x,y,z,zfac real, dimension(4) :: Limits zfac = 1.0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! print*,'in: om (t),h,Limits(1:4)',om,'(',2*pi/om,')',h,Limits(1), ! & Limits(2),Limits(3),Limits(4) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! check thickness is OK y = h if (h.gt.hmax) then !print*, 'TW (AD1): h > ',hmax, '; h = ',h y = hmax ! call xcstop('attenAD1') elseif (h.lt.hmin) then !use linear interpolation for small h zfac = h/hmin y = hmin end if !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! check 'x' is OK if (X_INPUT.eq.1) then x = (2*pi)/om !xstr = 'period'; elseif (X_INPUT.eq.2) then x = om/(2*pi) !xstr = 'freq'; else x = om !xstr = 'omega'; end if if (x.gt.xmax) then print*, 'x > ',xmax, '; x = ',x call xcstop('attenAD1') elseif (x.lt.xmin) then print*, 'x < ',xmin, '; x = ',x call xcstop('attenAD1') end if !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! Do interpolation Limits(1) = xmin Limits(2) = xmax Limits(3) = hmin Limits(4) = hmax call OP_chebinterp2d(z,x,y,chebys,Ncheb_f,Limits) !print*,'out from OP_chebinterp2d',exp(z) !! want to OUTPUT ENERGY attenuation coefficients: dampingAD1 = zfac*AC_INPUT*exp(z) end function!! dampingAD1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! real function fn_Kcheb(om,h) !! OUTPUTS WAVELENGTH, but uses chebyshev interp for k; !! NB infinite depth; real, intent(in) :: om,h integer, parameter :: Ncheb_f = 15! deg of poly in period/freq; integer, parameter :: Ncheb_h = 15! deg of poly in thickness; real, parameter :: g = 9.81! accel due to gravity; !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!SET INPUT TYPES: !X_INPUT = 1 for period; !X_INPUT = 2 for freq; !X_INPUT = 3 for omega; integer, parameter :: X_INPUT = 2 !AC_INPUT = 1 if input AC is energy AC; !AC_INPUT = 2 if input AC is amplitude AC; integer, parameter :: AC_INPUT = 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !RANGE OF VALIDITY VARIES FROM RUN TO RUN, ! SO NEED TO CHECK THIS EACH TIME; real, parameter :: hmax = 5.0 real, parameter :: hmin = 0.1 real, parameter :: xmin = 0.04!1/25 real, parameter :: xmax = 0.4!1/25 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! finally define vector of chebyshev coefficients from !! /home/nersc/timill/validation/miz1d/ALPfxns_23May2011/... !! ALP_kag_cheb2.m !! !! NB rest of function should be the same for all !! attenuation models; real, dimension((Ncheb_f+1)*(Ncheb_h+1)), parameter :: & chebys = (/ & 6.1178545e-02, & 4.7630410e-02, & -7.4259923e-03, & 3.2266473e-04, & 3.7901261e-04, & -1.1588517e-04, & 4.1972568e-05, & -5.9144780e-05, & 5.3706333e-05, & -3.4037453e-05, & 1.8068989e-05, & -7.3541636e-06, & 1.8531325e-06, & -1.0334165e-06, & 1.4510847e-06, & -8.7636342e-07, & -4.9508643e-02, & -5.1391714e-02, & 4.1553351e-03, & 3.3137764e-03, & -2.1700125e-03, & 5.9917414e-04, & -5.9214049e-05, & -1.6344748e-05, & 2.9132114e-05, & -3.8335656e-05, & 3.5439518e-05, & -2.6645713e-05, & 1.6178753e-05, & -5.2472743e-06, & -9.7492797e-07, & 1.8679828e-06, & 2.7291373e-02, & 3.1942152e-02, & 6.7997197e-05, & -3.7187010e-03, & 1.5280376e-03, & 5.9308700e-05, & -3.5186595e-04, & 2.0708785e-04, & -8.3156904e-05, & 2.4348386e-05, & 2.5108329e-06, & -1.0974554e-05, & 1.0475810e-05, & -1.0891658e-05, & 8.5787319e-06, & -5.7588071e-06, & -1.6405280e-02, & -2.0706880e-02, & -1.5359721e-03, & 3.0076357e-03, & -6.7668199e-04, & -5.0153322e-04, & 4.3988276e-04, & -1.2543626e-04, & -2.7321917e-05, & 5.4155528e-05, & -4.0730265e-05, & 2.0718965e-05, & -7.0046994e-06, & 1.1241325e-06, & -7.2378169e-07, & -1.0832924e-06, & 1.0309569e-02, & 1.3772135e-02, & 1.9192888e-03, & -2.1733218e-03, & 9.9952032e-05, & 6.0805845e-04, & -3.0666292e-04, & -2.5090785e-05, & 1.0730898e-04, & -6.8652688e-05, & 2.2360961e-05, & 5.0421273e-06, & -1.4366296e-05, & 1.1940474e-05, & -5.2569252e-06, & 1.9388250e-06, & -6.6251107e-03, & -9.2578293e-03, & -1.8206651e-03, & 1.4805456e-03, & 1.9539768e-04, & -5.4166234e-04, & 1.4488556e-04, & 1.2157015e-04, & -1.1791786e-04, & 3.4152369e-05, & 1.4792602e-05, & -2.4661119e-05, & 1.8083553e-05, & -6.3716163e-06, & -3.2669161e-07, & 3.1715416e-06, & 4.3400058e-03, & 6.3072789e-03, & 1.5779574e-03, & -9.6644444e-04, & -3.2246307e-04, & 4.1681554e-04, & -1.7925436e-05, & -1.5234846e-04, & 8.2856441e-05, & 8.7206238e-06, & -3.4928034e-05, & 2.2284403e-05, & -6.5409445e-06, & -4.6291050e-06, & 6.4741230e-06, & -6.0677142e-06, & -2.8602240e-03, & -4.2973115e-03, & -1.2789934e-03, & 6.0959398e-04, & 3.4108001e-04, & -2.9230442e-04, & -5.4241175e-05, & 1.4189498e-04, & -4.0765835e-05, & -3.6645251e-05, & 3.7064289e-05, & -9.7677898e-06, & -7.1088469e-06, & 1.0753666e-05, & -7.3248753e-06, & 3.3923980e-06, & 1.9036805e-03, & 2.9454096e-03, & 1.0100220e-03, & -3.7161811e-04, & -3.1243978e-04, & 1.8875755e-04, & 8.8937617e-05, & -1.1352913e-04, & 5.1433170e-06, & 4.7756235e-05, & -2.7230714e-05, & -3.8829904e-06, & 1.4215705e-05, & -8.8818281e-06, & 3.0079238e-06, & 1.0584637e-06, & -1.2654101e-03, & -2.0069310e-03, & -7.7142413e-04, & 2.1828396e-04, & 2.5987874e-04, & -1.1259011e-04, & -9.5754447e-05, & 8.1879830e-05, & 1.6552815e-05, & -4.6162343e-05, & 1.4879830e-05, & 1.2585437e-05, & -1.4755180e-05, & 3.5671619e-06, & 2.2386246e-06, & -4.2272123e-06, & 8.3681012e-04, & 1.3554121e-03, & 5.7297708e-04, & -1.2219149e-04, & -2.0281130e-04, & 6.1109289e-05, & 8.7352221e-05, & -5.4094958e-05, & -2.6424121e-05, & 3.8080691e-05, & -4.4310688e-06, & -1.6007225e-05, & 1.1348632e-05, & 1.6633267e-06, & -5.8844342e-06, & 4.9508328e-06, & -5.5292136e-04, & -9.1265624e-04, & -4.1923767e-04, & 6.3653708e-05, & 1.5255007e-04, & -2.8695646e-05, & -7.3278058e-05, & 3.2698729e-05, & 2.8804529e-05, & -2.8309388e-05, & -2.8275334e-06, & 1.5860746e-05, & -6.6788286e-06, & -5.2141738e-06, & 7.3338471e-06, & -3.7589059e-06, & 3.4713231e-04, & 5.8153359e-04, & 2.8437757e-04, & -3.0103527e-05, & -1.0488350e-04, & 1.0944782e-05, & 5.4031580e-05, & -1.7858610e-05, & -2.4371164e-05, & 1.8707729e-05, & 5.5597217e-06, & -1.2735839e-05, & 3.0134030e-06, & 5.9234022e-06, & -6.5319620e-06, & 2.1507619e-06, & -2.2359942e-04, & -3.8059985e-04, & -1.9802957e-04, & 1.1338773e-05, & 7.3417980e-05, & -1.2576875e-06, & -4.0079341e-05, & 8.4143130e-06, & 2.0250160e-05, & -1.1746141e-05, & -6.9352272e-06, & 9.8841441e-06, & -2.3549971e-07, & -5.8548701e-06, & 5.2265108e-06, & -4.1831097e-07, & 1.1591372e-04, & 1.9928658e-04, & 1.0755040e-04, & -3.0803884e-06, & -3.9919928e-05, & -1.3596812e-06, & 2.2469746e-05, & -3.1873999e-06, & -1.2015005e-05, & 5.8642093e-06, & 4.7920773e-06, & -5.6437551e-06, & -5.5823936e-07, & 3.7334566e-06, & -3.0039024e-06, & -2.1520886e-07, & -7.0667294e-05, & -1.2330442e-04, & -6.9775011e-05, & -8.6094602e-07, & 2.5823015e-05, & 2.7318635e-06, & -1.4990901e-05, & 6.0350409e-07, & 8.6016601e-06, & -3.0858163e-06, & -4.0859101e-06, & 3.7383264e-06, & 1.1460049e-06, & -2.8091025e-06, & 1.8601025e-06, & 7.7113350e-07 & /) !end of 'chebys' matrix (k) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! rest should be same for each program; real, parameter :: pi = 3.1415927 real :: x,y,z,zfac,kwtr real, dimension(4) :: Limits zfac = 1.0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! print*,'in: om (t),h,Limits(1:4)',om,'(',2*pi/om,')',h,Limits(1), ! & Limits(2),Limits(3),Limits(4) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! check thickness is OK y = h if (h.gt.hmax) then !print*, 'h > ',hmax, '; h = ',h y = hmax ! call xcstop('wave number') elseif (h.lt.hmin) then !use linear interpolation for small h zfac = h/hmin y = hmin end if !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! check 'x' is OK if (X_INPUT.eq.1) then x = (2*pi)/om !xstr = 'period'; elseif (X_INPUT.eq.2) then x = om/(2*pi) !xstr = 'freq'; else x = om !xstr = 'omega'; end if if (x.gt.xmax) then print*, 'x > ',xmax, '; x = ',x call xcstop('attenE0') elseif (x.lt.xmin) then print*, 'x < ',xmin, '; x = ',x call xcstop('attenE0') end if !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! Do interpolation Limits(1) = xmin Limits(2) = xmax Limits(3) = hmin Limits(4) = hmax call OP_chebinterp2d(z,x,y,chebys,Ncheb_f,Limits) !print*,'out from OP_chebinterp2d',exp(z) !! want to OUTPUT WAVELENGTH; kwtr = om**2/g if (zfac.lt.1.0) then fn_Kcheb = 2*pi/(kwtr+zfac*(z-kwtr)) else fn_Kcheb = 2*pi/z end if end function!! fn_Kcheb !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! real function fn_AGcheb(om,h) !! OUTPUTS WAVELENGTH, but uses chebyshev interp for k; !! NB infinite depth; real, intent(in) :: om,h integer, parameter :: Ncheb_f = 15! deg of poly in period/freq; integer, parameter :: Ncheb_h = 15! deg of poly in thickness; real, parameter :: g = 9.81! accel due to gravity; !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!SET INPUT TYPES: !X_INPUT = 1 for period; !X_INPUT = 2 for freq; !X_INPUT = 3 for omega; integer, parameter :: X_INPUT = 2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !RANGE OF VALIDITY VARIES FROM RUN TO RUN, ! SO NEED TO CHECK THIS EACH TIME; real, parameter :: hmax = 5.0 real, parameter :: hmin = 0.1 real, parameter :: xmin = 0.04!1/25 real, parameter :: xmax = 0.4!1/2.5 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! finally define vector of chebyshev coefficients from !! /home/nersc/timill/validation/miz1d/ALPfxns_23May2011/... !! ALP_kag_cheb2.m !! !! NB rest of function should be the same for all !! attenuation models; real, dimension((Ncheb_f+1)*(Ncheb_h+1)), parameter :: & chebys = (/ & 5.1476875e+01, & 3.5851287e+01, & -2.0161371e+00, & -1.4837287e+00, & 1.6490688e+00, & -1.1975536e+00, & 7.3526817e-01, & -3.9658347e-01, & 1.8234219e-01, & -6.4078222e-02, & 7.7216554e-03, & 1.4029248e-02, & -1.7316756e-02, & 1.5338240e-02, & -9.1825204e-03, & 5.9979845e-03, & 3.6502424e+01, & 2.9531673e+01, & -6.6663840e+00, & 2.0355269e+00, & -1.2501887e-01, & -5.0881454e-01, & 5.5884175e-01, & -4.1665714e-01, & 2.4768632e-01, & -1.1558954e-01, & 3.1503723e-02, & 1.2749817e-02, & -2.6217869e-02, & 2.7173690e-02, & -1.6964660e-02, & 1.1741352e-02, & -4.2415753e+00, & -5.3811072e+00, & 7.3541933e-01, & 5.6388848e-01, & -7.8082743e-01, & 5.0879122e-01, & -1.8714576e-01, & -1.8442216e-02, & 9.9250439e-02, & -1.0455480e-01, & 7.7373012e-02, & -4.3844903e-02, & 1.8859651e-02, & -1.3356739e-03, & -3.7038181e-03, & 6.5473543e-03, & 1.2052275e+00, & 1.8451298e+00, & 5.9175094e-02, & -5.0278925e-01, & 3.5066498e-01, & -5.2990562e-02, & -1.2191401e-01, & 1.4797341e-01, & -9.5831221e-02, & 3.3622633e-02, & 8.4185472e-03, & -2.6773926e-02, & 2.7437315e-02, & -2.1633818e-02, & 1.2176758e-02, & -6.3237138e-03, & -4.4886591e-01, & -7.9652916e-01, & -1.5211732e-01, & 2.9115610e-01, & -9.9368017e-02, & -8.3342481e-02, & 1.2677035e-01, & -7.1590978e-02, & 4.8332277e-03, & 3.2258123e-02, & -3.7220326e-02, & 2.4276868e-02, & -1.0009270e-02, & -1.5183783e-03, & 3.9915427e-03, & -5.5654531e-03, & 1.9232676e-01, & 3.8241241e-01, & 1.3631155e-01, & -1.5994025e-01, & 6.8458477e-03, & 8.6239124e-02, & -6.8130606e-02, & 6.0763650e-03, & 3.2535772e-02, & -3.4647419e-02, & 1.7075322e-02, & 1.4367773e-03, & -1.0696378e-02, & 1.2384502e-02, & -7.9740973e-03, & 4.5022135e-03, & -8.8290652e-02, & -1.9820132e-01, & -1.0260593e-01, & 8.3361365e-02, & 2.3420927e-02, & -6.2964818e-02, & 2.5353989e-02, & 1.9112081e-02, & -3.1237503e-02, & 1.6571002e-02, & 2.2279228e-03, & -1.1574224e-02, & 1.1755152e-02, & -5.9888024e-03, & 1.7303169e-03, & 1.7394416e-03, & 4.1308215e-02, & 1.0851794e-01, & 7.1171670e-02, & -4.0769685e-02, & -2.9688390e-02, & 4.0782334e-02, & -3.7558719e-03, & -2.2414734e-02, & 1.9336283e-02, & -2.0063403e-03, & -9.5812274e-03, & 9.5951096e-03, & -4.4404231e-03, & -1.9206747e-03, & 3.3149056e-03, & -3.9226502e-03, & -1.9216063e-02, & -6.0975981e-02, & -4.8254762e-02, & 1.8731780e-02, & 2.6853063e-02, & -2.3930689e-02, & -5.4905477e-03, & 1.8199186e-02, & -8.5579193e-03, & -5.1686027e-03, & 9.2125307e-03, & -4.4419757e-03, & -1.4439083e-03, & 4.6439660e-03, & -3.8257305e-03, & 2.2626609e-03, & 8.8137611e-03, & 3.4219580e-02, & 3.2502492e-02, & -8.0943701e-03, & -2.0947420e-02, & 1.2600898e-02, & 8.2558699e-03, & -1.2737563e-02, & 1.8772807e-03, & 7.0407214e-03, & -6.2698834e-03, & 5.0930571e-04, & 3.7712139e-03, & -3.7441257e-03, & 2.0374799e-03, & 9.3991658e-05, & -3.9681945e-03, & -1.8936057e-02, & -2.1906068e-02, & 3.1344733e-03, & 1.4979860e-02, & -5.7233010e-03, & -8.0588383e-03, & 8.1048388e-03, & 1.3997287e-03, & -6.3050107e-03, & 3.3289193e-03, & 1.5164958e-03, & -3.6650127e-03, & 1.7821583e-03, & -1.0537692e-04, & -1.4484547e-03, & 1.7318914e-03, & 1.0348464e-02, & 1.4856272e-02, & -8.2448193e-04, & -1.0173182e-02, & 1.9432952e-03, & 6.7816294e-03, & -4.7375501e-03, & -2.6037648e-03, & 4.7384130e-03, & -1.2100509e-03, & -2.1819679e-03, & 2.5788828e-03, & -1.4663512e-04, & -1.0963861e-03, & 1.6882532e-03, & -7.1724650e-04, & -5.4635319e-03, & -9.5655341e-03, & -9.6147013e-05, & 6.3875418e-03, & -2.4359363e-04, & -4.9346444e-03, & 2.5171319e-03, & 2.4727245e-03, & -3.0913108e-03, & 1.0659441e-04, & 1.9460660e-03, & -1.4640464e-03, & -6.1646589e-04, & 1.3771327e-03, & -1.3085087e-03, & 2.5239071e-04, & 2.9409511e-03, & 6.3268055e-03, & 5.2889100e-04, & -4.0303006e-03, & -4.7745904e-04, & 3.5446164e-03, & -1.1745115e-03, & -2.1076058e-03, & 1.8926369e-03, & 4.6777170e-04, & -1.5674525e-03, & 6.2010347e-04, & 9.2288824e-04, & -1.2913578e-03, & 7.6688307e-04, & -6.7863261e-05, & -1.3543689e-03, & -3.3338923e-03, & -4.3428975e-04, & 2.0511680e-03, & 4.4845883e-04, & -1.9439579e-03, & 4.4513624e-04, & 1.2523133e-03, & -9.2518088e-04, & -4.2893467e-04, & 9.0391982e-04, & -1.8640283e-04, & -6.5408807e-04, & 7.9079259e-04, & -3.3622583e-04, & -1.1412539e-05, & 7.2495716e-04, & 2.0658833e-03, & 4.2496864e-04, & -1.2121821e-03, & -4.2114508e-04, & 1.2426069e-03, & -1.0286042e-04, & -8.7654352e-04, & 4.7497350e-04, & 4.2596918e-04, & -5.9685538e-04, & -3.9986196e-05, & 5.1960506e-04, & -5.2769769e-04, & 7.7011440e-05 & /) !end of 'chebys' matrix (ag) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! rest should be same for each program; real, parameter :: pi = 3.1415927 real :: x,y,z,zfac,agw real, dimension(4) :: Limits zfac = 1.0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! print*,'in: om (t),h,Limits(1:4)',om,'(',2*pi/om,')',h,Limits(1), ! & Limits(2),Limits(3),Limits(4) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! check thickness is OK y = h if (h.gt.hmax) then !print*, 'h > ',hmax, '; h = ',h y = hmax ! call xcstop('group vel') elseif (h.lt.hmin) then !use linear interpolation for small h zfac = h/hmin y = hmin end if !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! check 'x' is OK if (X_INPUT.eq.1) then x = (2*pi)/om !xstr = 'period'; elseif (X_INPUT.eq.2) then x = om/(2*pi) !xstr = 'freq'; else x = om !xstr = 'omega'; end if if (x.gt.xmax) then print*, 'x > ',xmax, '; x = ',x call xcstop('attenE0') elseif (x.lt.xmin) then print*, 'x < ',xmin, '; x = ',x call xcstop('attenE0') end if !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! Do interpolation Limits(1) = xmin Limits(2) = xmax Limits(3) = hmin Limits(4) = hmax call OP_chebinterp2d(z,x,y,chebys,Ncheb_f,Limits) !print*,'out from OP_chebinterp2d',exp(z) agw = g/2.0/om if (zfac.lt.1.0) then fn_AGcheb = agw+zfac*(z-agw) else fn_AGcheb = z end if end function!! fn_AGcheb !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #if defined (WAVES_ROLLOVER) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine OP_chebinterp1d(z,x,chebys,Ncx,Limits) real, intent(out) :: z real, intent(in) :: x integer, intent(in) :: Ncx !! real, intent(in), dimension(:) :: chebys real, intent(in), dimension(2) :: Limits ! real, intent(in) :: chebys, Limits !! real :: Tm0_x, Tm1_x, Tm_x real :: tx, xmin, xmax integer :: nx !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !print*, 'hi - using chebys' !!print*,'TWtest0, Ncy = ',Ncy,Ncx,size(chebys),size(chebys)/(Ncx+1) xmin = Limits(1) xmax = Limits(2) tx = -1.0+2*(x-xmin)/(xmax-xmin) !! z = chebys(1) Tm0_x = 1.0 Tm1_x = tx !! if (Ncx.ge.1) then z = z+chebys(2)*Tm1_x end if do nx=2,Ncx Tm_x = 2*tx*Tm1_x-Tm0_x z = z+chebys(nx+1)*Tm_x Tm0_x = Tm1_x Tm1_x = Tm_x end do end subroutine OP_chebinterp1d !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #endif end module mod_attenuation_ice