attenuation.f90 1.34 KB

!____________________________________________________________________
!
!DESCRIPTION: In this routine, the attenuation coefficient is
!             calculated according Kohout and Meylan (2008) and the
!             spectrum attenuation is computed.
!
!____________________________________________________________________

!INTERFACE:
subroutine attenuation

!USES:
use parameters

implicit none

double precision :: q11,q12,q13,q14,q15,q21,q22,q23,q24,q25
double precision :: q31,q32,q33,q34,q35
double precision, allocatable :: p1(:),p2(:),p3(:)
double precision, allocatable :: att(:)

allocate(p1 (nf))
allocate(p2 (nf))
allocate(p3 (nf))
allocate(att(nf))


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

alpha(i,:)=-1*(p1*h(i)**2 + p2*h(i) + p3)

do ii=1,nf
   if ( alpha(i,ii).lt.0 )then
      alpha(i,ii)=0d0
   end if
end do

att=C_ice(i)*alpha(i,:)/(Dave(i)+3e-14)
!S_ice(n,i,1:nf)=-att*E(n,i,1:nf)
S_ice(i,:,n)=-att*E(i,:,n)
E(i,:,n)=E(i,:,n)*exp(-att*Cg*dt)

end subroutine attenuation