module mod_math use mod_parameters use mod_variables use mod_matrix_coeffs implicit none type :: div_s_c real(8) :: jm1_dsj0, jm1_sj0, jm1_dsj2, jm1_sj2, jm1_dsjm2, jm1_sjm2 real(8) :: jp1_dsj0, jp1_sj0, jp1_dsj2, jp1_sj2, jp1_dsjp2, jp1_sjp2 real(8) :: jm0_dsjm1, jm0_sjm1, jm0_dsjp1, jm0_sjp1 real(8) :: jm1_rho, jp1_rho end type type(div_s_c) :: dsc contains subroutine peaks(tt, it2) implicit none real(8) :: tt integer :: it2 real(8) :: Delta_20_1 , Delta_20_2 real(8) :: Delta_r22_1, Delta_r22_2 real(8) :: Delta_i22_1, Delta_i22_2 real(8) :: Delta_d_1 , Delta_d_2 if(it2<3) go to 10 gp20_2 = k20; rgp22_2 = real_k22; igp22_2 = imag_k22; dissipation_2 = layer(layer_ocean)%dissipation*W_to_GW Delta_20_1 = gp20_2 - gp20_1 Delta_20_2 = gp20_1 - gp20_0 Delta_r22_1 = rgp22_2 - rgp22_1 Delta_r22_2 = rgp22_1 - rgp22_0 Delta_i22_1 = igp22_2 - igp22_1 Delta_i22_2 = igp22_1 - igp22_0 Delta_d_1 = dissipation_2 - dissipation_1 Delta_d_2 = dissipation_1 - dissipation_0 ! if((Delta_20_1<0) .and.(Delta_20_2>0 )) print*, '20 coefficient peak is at...', tt, gp20_2 if((Delta_r22_1<0).and.(Delta_r22_2>0)) then ! print*, '22 real coefficient peak is at...', (tt-time_step)/rotation_period, rgp22_1 peak_time_1 = peak_time_2 peak_time_2 = tt ! print*, 'Peak time difference...', (peak_time_2-peak_time_1)/rotation_period end if if((Delta_d_1<0).and.(Delta_d_2>0)) then ! print*, 'Dissipaiton peak at...', (tt-time_step)/rotation_period, dissipation_1 pt_d_1 = pt_d_2 pt_d_2 = tt ! print*, 'Peak time difference (dissipation)...', (pt_d_2-pt_d_1)/rotation_period end if gp20_0 = gp20_1 gp20_1 = gp20_2 rgp22_0 = rgp22_1 rgp22_1 = rgp22_2 igp22_0 = igp22_1 igp22_1 = igp22_2 dissipation_0 = dissipation_1 dissipation_1 = dissipation_2 10 continue end subroutine peaks subroutine determine_dt implicit none integer :: j, m, ir, il, nn integer :: jm1, jm0, jp1 real(8) :: vr, vt, dtmin_r, dtmin_t, erjm1, erjp1, jj, etjm1, etjp1 real(8) :: cr, ct, elas_wave, a, ur dtmin_r = 1d30; dtmin_t = 1d30; elas_wave=1.0d30 do j = jstart, jfinal jj = dble(j) erjm1 = dsqrt( jj /(2.0d0*jj+1.0d0)) erjp1 = -dsqrt((jj+1.0d0)/(2.0d0*jj+1.0d0)) etjm1 = (jj+1.0d0) / (2.0d0*jj+1.0d0) etjp1 = dsqrt(jj*(jj+1.0d0)) / (2.0d0*jj+1.0d0) do m = 0, min(j,2), 2 jm1 = jml2(j,m,j-1) jm0 = jml2(j,m,j ) jp1 = jml2(j,m,j+1) do il = numoflayers, numoflayers nn = layer(il)%n do ir = 1, nn vr = abs(erjm1 * layer(il)%vel(ir, jm1) + erjp1 * layer(il)%vel(ir, jp1)) ur = abs(erjm1 * layer(il)%dis(ir, jm1) + erjp1 * layer(il)%dis(ir, jp1)) vt = abs(etjm1 * layer(il)%vel(ir, jm1) + etjp1 * layer(il)%vel(ir, jp1)) * dsqrt(jj*(jj+1.0d0)) a = abs(layer(il)%r_l(ir+1)-layer(il)%r_l(ir))/(dsqrt(layer(il)%mu_i(ir)/layer(il)%rho_i(ir)) * (ur)) cr = abs(layer(il)%r_l(ir+1)-layer(il)%r_l(ir)) / vr ct = layer(il)%r_l(ir) / vt ! print*, vr, cr, rotation_period/100.0d0 if(cr < dtmin_r) dtmin_r = cr if(ct < dtmin_t) dtmin_t = ct if(a < elas_wave) elas_wave = a end do end do end do end do ! time_step = CFL * min(dtmin_r, dtmin_t, elas_wave) ! print*, 'aaa', dtmin_r, dtmin_t, elas_wave end subroutine determine_dt subroutine Divergence_tensor implicit none integer :: j, m, il, ir, nn integer :: jm1, jm0, jp1, jms complex(8) :: dsj0, asj0, dsj2, asj2, dsjm2, asjm2, dsjp2, asjp2, rho_jm1, rho_jp1, dsjm1, asjm1, dsjp1, asjp1 real(8) :: hi, hup2, hlo2, hup2_r, hlo2_r do il = 1, numoflayers nn = layer(il)%n do ir = 1, nn-1 hi = layer(il)%r_i(ir+1) - layer(il)%r_i(ir ) hup2 = (layer(il)%r_i(ir+1) - layer(il)%r_l(ir+1)) / hi; hup2_r = hup2/layer(il)%r_l(ir+1) hlo2 = (layer(il)%r_l(ir+1) - layer(il)%r_i(ir )) / hi; hlo2_r = hlo2/layer(il)%r_l(ir+1) do j = 1, jmax call div_tensor_coeffs(dble(j)) do m = 0, min(2,j), 2 jm1 = jml2(j,m,j-1) jm0 = jml2(j,m,j ) jp1 = jml2(j,m,j+1) jms = (jm2(j,m) - 1)*6 dsj0 = dsc%jm1_dsj0 * (layer(il)%str(ir+1, jms+1) - layer(il)%str(ir, jms+1)) / hi asj0 = dsc%jm1_sj0 * (layer(il)%str(ir+1, jms+1)*hlo2_r + layer(il)%str(ir, jms+1)*hup2_r) dsjm2 = dsc%jm1_dsjm2 * (layer(il)%str(ir+1, jms+2) - layer(il)%str(ir, jms+2)) / hi asjm2 = dsc%jm1_sjm2 * (layer(il)%str(ir+1, jms+2)*hlo2_r + layer(il)%str(ir, jms+2)*hup2_r) dsj2 = dsc%jm1_dsj2 * (layer(il)%str(ir+1, jms+4) - layer(il)%str(ir, jms+4)) / hi asj2 = dsc%jm1_sj2 * (layer(il)%str(ir+1, jms+4)*hlo2_r + layer(il)%str(ir, jms+4)*hup2_r) rho_jm1 = dsc%jm1_rho * (layer(il)%g_l(ir+1) * (layer(il)%rhoe(ir+1, jm2(j,m))*hlo2 + layer(il)%rhoe(ir, jm2(j,m))*hup2)) layer(il)%div_tensor(ir, jm1) = dsj0 + asj0 + dsj2 + asj2 + dsjm2 + asjm2 + rho_jm1 dsj0 = dsc%jp1_dsj0 * (layer(il)%str(ir+1, jms+1) - layer(il)%str(ir, jms+1)) / hi asj0 = dsc%jp1_sj0 * (layer(il)%str(ir+1, jms+1)*hlo2_r + layer(il)%str(ir, jms+1)*hup2_r) dsjp2 = dsc%jp1_dsjp2 * (layer(il)%str(ir+1, jms+6) - layer(il)%str(ir, jms+6)) / hi asjp2 = dsc%jp1_sjp2 * (layer(il)%str(ir+1, jms+6)*hlo2_r + layer(il)%str(ir, jms+6)*hup2_r) dsj2 = dsc%jp1_dsj2 * (layer(il)%str(ir+1, jms+4) - layer(il)%str(ir, jms+4)) / hi asj2 = dsc%jp1_sj2 * (layer(il)%str(ir+1, jms+4)*hlo2_r + layer(il)%str(ir, jms+4)*hup2_r) rho_jp1 = dsc%jp1_rho * (layer(il)%g_l(ir+1) * (layer(il)%rhoe(ir+1, jm2(j,m))*hlo2 + layer(il)%rhoe(ir, jm2(j,m))*hup2)) layer(il)%div_tensor(ir, jp1) = dsj0 + asj0 + dsj2 + asj2 + dsjp2 + asjp2 + rho_jp1 dsjm1 = dsc%jm0_dsjm1 * (layer(il)%str(ir+1, jms+3) - layer(il)%str(ir, jms+3)) / hi asjm1 = dsc%jm0_sjm1 * (layer(il)%str(ir+1, jms+3)*hlo2_r + layer(il)%str(ir, jms+3)*hup2_r) dsjp1 = dsc%jm0_dsjp1 * (layer(il)%str(ir+1, jms+5) - layer(il)%str(ir, jms+5)) / hi asjp1 = dsc%jm0_sjp1 * (layer(il)%str(ir+1, jms+5)*hlo2_r + layer(il)%str(ir, jms+5)*hup2_r) layer(il)%div_tensor(ir, jm0) = dsjm1 + asjm1 + dsjp1 + asjp1 end do end do end do end do end subroutine Divergence_tensor subroutine div_tensor_coeffs(j) implicit none real(8) :: j ! !Y_{JM}^{J-1} dsc%jm1_dsj0 = -dsqrt( j / (3.0d0 * (2.0d0 * j + 1.0d0))) dsc%jm1_dsjm2 = dsqrt((j - 1.0d0) / (2.0d0 * j - 1.0d0)) dsc%jm1_dsj2 = -dsqrt( (j + 1.0d0)*(2.0d0 * j + 3.0d0) /(6.d0*(2.d0*j - 1.d0)*(2.d0*j + 1.d0) )) dsc%jm1_sj0 = dsc%jm1_dsj0 * (j + 1.0d0) dsc%jm1_sjm2 = -dsc%jm1_dsjm2 * (j - 2.0d0) dsc%jm1_sj2 = dsc%jm1_dsj2 * (j + 1.0d0) dsc%jm1_rho = -dsqrt( j /(2.0d0 * j + 1.0d0)) ! !Y_{JM}^{J+1} dsc%jp1_dsj0 = dsqrt((j + 1.0d0)/((3.d0*(2.d0*j + 1.d0)))) dsc%jp1_dsjp2 = -dsqrt((j +2.d0) / (2.d0*j + 3.d0)) dsc%jp1_dsj2 = dsqrt(( j*(2.d0*j - 1.d0) )/(6.d0*(2.d0*j + 3.d0)*(2.d0*j + 1.d0))) dsc%jp1_sj0 = -dsc%jp1_dsj0 * j dsc%jp1_sjp2 = dsc%jp1_dsjp2 * (j + 3.0d0) dsc%jp1_sj2 = -dsc%jp1_dsj2 * j dsc%jp1_rho = -(-dsqrt((j+1.0d0)/(2.0d0 * j + 1.0d0))) dsc%jm0_dsjm1 = dsqrt((j - 1.0d0)/(2.0d0*(2.0d0 * j + 1.0d0))) dsc%jm0_sjm1 = -(j - 1.0d0) * dsc%jm0_dsjm1 dsc%jm0_dsjp1 = -dsqrt((j + 2.0d0)/(2.0d0*(2.0d0 * j + 1.0d0))) dsc%jm0_sjp1 = (j + 2.0d0) * dsc%jm0_dsjp1 end subroutine div_tensor_coeffs end module mod_math