attenuation.f90 1.45 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