module mod_matrix use mod_allocation use mod_variables use mod_parameters use mod_matrix_coeffs use mod_LU implicit none contains subroutine create_matrix(j) implicit none integer :: j real(8) :: jj, d,d2 jj = dble(j) call matrix_coeffs(jj) S_U=zero; T_U=zero; S_L=zero; T_L=zero; S_I=int(zero); T_I=int(zero) if(SWITCH_compressible==0)then call generate_matrix_incompressible call internal_boundaries call top_bottom_boundaries elseif(SWITCH_compressible==1)then call generate_matrix_compressible call internal_boundaries call top_bottom_boundaries else stop "For incompressible: 0 or Compressible: 1..." end if call bandec(S_U, smtot, sph_lo_diag, sph_hi_diag, smtot, sph_band, S_L, sph_lo_diag, S_I, d) call bandec(T_U, tmtot, tor_lo_diag, tor_hi_diag, tmtot, tor_band, T_L, tor_lo_diag, T_I, d2) end subroutine create_matrix subroutine generate_matrix_incompressible implicit none integer :: ir, il, smm, nn, tmm real(8) :: hl, hup, hlo, hup_r, hlo_r real(8) :: hi, hup2, hlo2, hup2_r, hlo2_r smm = 0; tmm=0 do il = 1, numoflayers nn = layer(il)%n do ir = 1, nn hl = layer(il)%r_l(ir+1) - layer(il)%r_l(ir) hup = (layer(il)%r_l(ir+1) - layer(il)%r_i(ir)) / hl; hup_r = hup/layer(il)%r_i(ir) hlo = (layer(il)%r_i(ir) - layer(il)%r_l(ir)) / hl; hlo_r = hlo/layer(il)%r_i(ir) !! Balance of mass S_U(smm+7*(ir-1)+1, sph_diag+0) = layer(il)%rho_r_i(ir)*sc%mass_bal_vjm1*hup ! !V_{JM}^{J-1}, LOWER S_U(smm+7*(ir-1)+1, sph_diag+1) = layer(il)%rho_r_i(ir)*sc%mass_bal_vjp1*hup ! !V_{JM}^{J+1}, LOWER S_U(smm+7*(ir-1)+1, sph_diag+2) = 1.0d0/time_step ! !RHO_t S_U(smm+7*(ir-1)+1, sph_diag+7) = layer(il)%rho_r_i(ir)*sc%mass_bal_vjm1*hlo ! !V_{JM}^{J-1}, UPPER S_U(smm+7*(ir-1)+1, sph_diag+8) = layer(il)%rho_r_i(ir)*sc%mass_bal_vjp1*hlo ! !V_{JM}^{J+1}, UPPER !! Constitutive - pressure div(v) = 0 S_U(smm+7*(ir-1)+2, sph_diag-1) = -sc%const_p_vjm1(1)/hl + sc%const_p_vjm1(2)*hup_r S_U(smm+7*(ir-1)+2, sph_diag+0) = -sc%const_p_vjp1(1)/hl + sc%const_p_vjp1(2)*hup_r S_U(smm+7*(ir-1)+2, sph_diag+6) = sc%const_p_vjm1(1)/hl + sc%const_p_vjm1(2)*hlo_r S_U(smm+7*(ir-1)+2, sph_diag+7) = sc%const_p_vjp1(1)/hl + sc%const_p_vjp1(2)*hlo_r !! CONSTITUTIVE EQUATION Y_{JM}^{J-2,2} S_U(smm+7*(ir-1)+5, sph_diag-4) = (-sc%Const_jm2_vjm1(1)/hl + sc%Const_jm2_vjm1(2)*hup_r) S_U(smm+7*(ir-1)+5, sph_diag+0) = 1.0d0 / (2.0d0 * layer(il)%mu_i(ir)*time_step) + 1.0d0/(2.0d0*layer(il)%eta_i(ir)) S_U(smm+7*(ir-1)+5, sph_diag+3) = ( sc%Const_jm2_vjm1(1)/hl + sc%Const_jm2_vjm1(2)*hlo_r) !! CONSTITUTIVE EQUATION Y_{JM}^{J ,2} S_U(smm+7*(ir-1)+6, sph_diag-5) = (-sc%Const_j2__vjm1(1)/hl + sc%Const_j2__vjm1(2)*hup_r) S_U(smm+7*(ir-1)+6, sph_diag-4) = (-sc%Const_j2__vjp1(1)/hl + sc%Const_j2__vjp1(2)*hup_r) S_U(smm+7*(ir-1)+6, sph_diag+0) = 1.0d0 / (2.0d0 * layer(il)%mu_i(ir)*time_step) + 1.0d0/(2.0d0*layer(il)%eta_i(ir)) S_U(smm+7*(ir-1)+6, sph_diag+2) = ( sc%Const_j2__vjm1(1)/hl + sc%Const_j2__vjm1(2)*hlo_r) S_U(smm+7*(ir-1)+6, sph_diag+3) = ( sc%Const_j2__vjp1(1)/hl + sc%Const_j2__vjp1(2)*hlo_r) !! CONSTITUTIVE EQUATION Y_{JM}^{J+2,2} S_U(smm+7*(ir-1)+7, sph_diag-5) = (-sc%Const_jp2_vjp1(1)/hl + sc%Const_jp2_vjp1(2)*hup_r) S_U(smm+7*(ir-1)+7, sph_diag+0) = 1.0d0 / (2.0d0 * layer(il)%mu_i(ir)*time_step) + 1.0d0/(2.0d0*layer(il)%eta_i(ir)) S_U(smm+7*(ir-1)+7, sph_diag+2) = ( sc%Const_jp2_vjp1(1)/hl + sc%Const_jp2_vjp1(2)*hlo_r) if(SWITCH_ocean_status == 'v')then ! !CONSTITUTIVE EQUATION Y_{JM}^{J-1,2} T_U(tmm+3*(ir-1) + 2, tor_diag-1) = (-tc%Const_jm1_vj(1)/hl + tc%Const_jm1_vj(2)*hup_r) T_U(tmm+3*(ir-1) + 2, tor_diag+0) = 1.0d0 / (2.0d0 * layer(il)%mu_i(ir)*time_step) + 1.0d0/(2.0d0*layer(il)%eta_i(ir)) T_U(tmm+3*(ir-1) + 2, tor_diag+2) = ( tc%Const_jm1_vj(1)/hl + tc%Const_jm1_vj(2)*hlo_r) ! !CONSTITUTIVE EQUATION Y_{JM}^{J+1,2} T_U(tmm+3*(ir-1) + 3, tor_diag-2) = (-tc%Const_jp1_vj(1)/hl + tc%Const_jp1_vj(2)*hup_r) T_U(tmm+3*(ir-1) + 3, tor_diag+0) = 1.0d0 / (2.0d0 * layer(il)%mu_i(ir)*time_step) + 1.0d0/(2.0d0*layer(il)%eta_i(ir)) T_U(tmm+3*(ir-1) + 3, tor_diag+1) = ( tc%Const_jp1_vj(1)/hl + tc%Const_jp1_vj(2)*hlo_r) end if end do smm = sum(layer(1:il)%sm); tmm=sum(layer(1:il)%tm) end do smm = 0; tmm=0 do il = 1, numoflayers nn = layer(il)%n do ir = 2, nn hi = layer(il)%r_i(ir) - layer(il)%r_i(ir-1) hup2 = (layer(il)%r_i(ir) - layer(il)%r_l( ir)) / hi; hup2_r = hup2/layer(il)%r_l(ir) hlo2 = (layer(il)%r_l(ir) - layer(il)%r_i(ir-1)) / hi; hlo2_r = hlo2/layer(il)%r_l(ir) !!Y_{JM}^{J-1} COMPONENT S_U(smm+7*(ir-1)+3, sph_diag-7) = (CN_theta)*(layer(il)%g_l(ir)) * sc%Mom__jm1_rho * hup2 ! !RHO_{JM} S_U(smm+7*(ir-1)+3, sph_diag-6) = (CN_theta)*(-sc%Mom__jm1_sj0(1) /hi + sc%Mom__jm1_sj0(2) * hup2_r) ! !S_{JM}^{J ,0} S_U(smm+7*(ir-1)+3, sph_diag-5) = (CN_theta)*(-sc%Mom__jm1_sjm2(1)/hi + sc%Mom__jm1_sjm2(2) * hup2_r) ! !S_{JM}^{J-2,2} S_U(smm+7*(ir-1)+3, sph_diag-4) = (CN_theta)*(-sc%Mom__jm1_sj2(1) /hi + sc%Mom__jm1_sj2(2) * hup2_r) ! !S_{JM}^{J ,2} S_U(smm+7*(ir-1)+3, sph_diag+0) = (CN_theta)*(layer(il)%g_l(ir)) * sc%Mom__jm1_rho * hlo2 ! !RHO_{JM} S_U(smm+7*(ir-1)+3, sph_diag+1) = (CN_theta)*( sc%Mom__jm1_sj0(1) /hi + sc%Mom__jm1_sj0(2) * hlo2_r) ! !S_{JM}^{J ,0} S_U(smm+7*(ir-1)+3, sph_diag+2) = (CN_theta)*( sc%Mom__jm1_sjm2(1)/hi + sc%Mom__jm1_sjm2(2) * hlo2_r) ! !S_{JM}^{J-2,2} S_U(smm+7*(ir-1)+3, sph_diag+3) = (CN_theta)*( sc%Mom__jm1_sj2(1) /hi + sc%Mom__jm1_sj2(2) * hlo2_r) ! !S_{JM}^{J ,2} S_U(smm+7*(ir-1)+3, sph_diag-2) = -SWITCH_inertia*layer(il)%rho_l(ir)/time_step !!Y_{JM}^{J+1} COMPONENT S_U(smm+7*(ir-1)+4, sph_diag-8) = (CN_theta)*(layer(il)%g_l(ir)) * sc%Mom__jp1_rho*hup2 ! !RHO_{JM} S_U(smm+7*(ir-1)+4, sph_diag-7) = (CN_theta)*(-sc%Mom__jp1_sj0(1) /hi + sc%Mom__jp1_sj0(2) * hup2_r) ! !S_{JM}^{J ,0} S_U(smm+7*(ir-1)+4, sph_diag-5) = (CN_theta)*(-sc%Mom__jp1_sj2(1) /hi + sc%Mom__jp1_sj2(2) * hup2_r) ! !S_{JM}^{J ,2} S_U(smm+7*(ir-1)+4, sph_diag-4) = (CN_theta)*(-sc%Mom__jp1_sjp2(1)/hi + sc%Mom__jp1_sjp2(2) * hup2_r) ! !S_{JM}^{J+2,2} S_U(smm+7*(ir-1)+4, sph_diag-1) = (CN_theta)*(layer(il)%g_l(ir)) * sc%Mom__jp1_rho*hlo2 ! !RHO_{JM} S_U(smm+7*(ir-1)+4, sph_diag+0) = (CN_theta)*( sc%Mom__jp1_sj0(1) /hi + sc%Mom__jp1_sj0(2) * hlo2_r) ! !S_{JM}^{J ,0} S_U(smm+7*(ir-1)+4, sph_diag+2) = (CN_theta)*( sc%Mom__jp1_sj2(1) /hi + sc%Mom__jp1_sj2(2) * hlo2_r) ! !S_{JM}^{J ,2} S_U(smm+7*(ir-1)+4, sph_diag+3) = (CN_theta)*( sc%Mom__jp1_sjp2(1)/hi + sc%Mom__jp1_sjp2(2) * hlo2_r) ! !S_{JM}^{J+2,2} S_U(smm+7*(ir-1)+4, sph_diag-2 ) = -SWITCH_inertia*layer(il)%rho_l(ir)/time_step if(SWITCH_ocean_status == 'v')then T_U(tmm+3*(ir-1)+1, tor_diag-2) = (CN_theta)*(-tc%Mom__j_s_jm1(1) /hi + tc%Mom__j_s_jm1(2) * hup2_r) T_U(tmm+3*(ir-1)+1, tor_diag-1) = (CN_theta)*(-tc%Mom__j_s_jp1(1) /hi + tc%Mom__j_s_jp1(2) * hup2_r) T_U(tmm+3*(ir-1)+1, tor_diag+0) = -SWITCH_inertia*layer(il)%rho_l(ir)/time_step T_U(tmm+3*(ir-1)+1, tor_diag+1) = (CN_theta)*( tc%Mom__j_s_jm1(1) /hi + tc%Mom__j_s_jm1(2) * hlo2_r) T_U(tmm+3*(ir-1)+1, tor_diag+2) = (CN_theta)*( tc%Mom__j_s_jp1(1) /hi + tc%Mom__j_s_jp1(2) * hlo2_r) end if end do smm = sum(layer(1:il)%sm); tmm=sum(layer(1:il)%tm) end do end subroutine generate_matrix_incompressible subroutine generate_matrix_compressible implicit none integer :: ir, il, smm, tmm, nn real(8) :: hl, hup, hlo, hup_r, hlo_r real(8) :: hi, hup2, hlo2, hup2_r, hlo2_r real(8) :: anelas smm = 0; tmm=0 do il = 1, numoflayers nn = layer(il)%n if(il==numoflayers-1)then anelas=0.0d0 else anelas=1.0d0 end if do ir = 1, nn hl = layer(il)%r_l(ir+1) - layer(il)%r_l(ir) hup = (layer(il)%r_l(ir+1) - layer(il)%r_i(ir)) / hl; hup_r = hup/layer(il)%r_i(ir) hlo = (layer(il)%r_i(ir) - layer(il)%r_l(ir)) / hl; hlo_r = hlo/layer(il)%r_i(ir) !!Balance of mass S_U(smm+7*(ir-1)+1, sph_diag+0) = layer(il)%rho_r_i(ir)*sc%mass_bal_vjm1*hup + layer(il)%rho_i(ir)*(-sc%const_p_vjm1(1)/hl + sc%const_p_vjm1(2)*hup_r) ! !V_{JM}^{J-1}, LOWER S_U(smm+7*(ir-1)+1, sph_diag+1) = layer(il)%rho_r_i(ir)*sc%mass_bal_vjp1*hup + layer(il)%rho_i(ir)*(-sc%const_p_vjp1(1)/hl + sc%const_p_vjp1(2)*hup_r) ! !V_{JM}^{J+1}, LOWER S_U(smm+7*(ir-1)+1, sph_diag+2) = anelas/time_step ! !RHO_t S_U(smm+7*(ir-1)+1, sph_diag+7) = layer(il)%rho_r_i(ir)*sc%mass_bal_vjm1*hlo + layer(il)%rho_i(ir)*( sc%const_p_vjm1(1)/hl + sc%const_p_vjm1(2)*hlo_r) ! !V_{JM}^{J-1}, UPPER S_U(smm+7*(ir-1)+1, sph_diag+8) = layer(il)%rho_r_i(ir)*sc%mass_bal_vjp1*hlo + layer(il)%rho_i(ir)*( sc%const_p_vjp1(1)/hl + sc%const_p_vjp1(2)*hlo_r) ! !V_{JM}^{J+1}, UPPER if(il==numoflayers-1)then !! Constitutive - pressure p - K rho/rho0 + zeta div(v) = 0 S_U(smm+7*(ir-1)+2, sph_diag+2) = pjmstr/(layer(il)%K_i(ir)) !! pressure S_U(smm+7*(ir-1)+2, sph_diag+1) = -1.0d0 /(layer(il)%rho_i(ir)) !! density S_U(smm+7*(ir-1)+2, sph_diag-1) = -0.0d0*(-sc%const_p_vjm1(1)/hl + sc%const_p_vjm1(2)*hup_r) S_U(smm+7*(ir-1)+2, sph_diag+0) = -0.0d0*(-sc%const_p_vjp1(1)/hl + sc%const_p_vjp1(2)*hup_r) S_U(smm+7*(ir-1)+2, sph_diag+6) = -0.0d0*( sc%const_p_vjm1(1)/hl + sc%const_p_vjm1(2)*hlo_r) S_U(smm+7*(ir-1)+2, sph_diag+7) = -0.0d0*( sc%const_p_vjp1(1)/hl + sc%const_p_vjp1(2)*hlo_r) else !! Constitutive - pressure p_t - vr rho0 g0 + K div(v) = 0 S_U(smm+7*(ir-1)+2, sph_diag+2) = pjmstr/(layer(il)%K_i(ir)*time_step) S_U(smm+7*(ir-1)+2, sph_diag-1) = (-sc%const_p_vjm1(1)/hl + sc%const_p_vjm1(2)*hup_r) - (layer(il)%rho_i(ir)*layer(il)%g_i(ir)*sc%mass_bal_vjm1*hup) * (1.0d0/layer(il)%K_i(ir)) S_U(smm+7*(ir-1)+2, sph_diag+0) = (-sc%const_p_vjp1(1)/hl + sc%const_p_vjp1(2)*hup_r) - (layer(il)%rho_i(ir)*layer(il)%g_i(ir)*sc%mass_bal_vjp1*hup) * (1.0d0/layer(il)%K_i(ir)) S_U(smm+7*(ir-1)+2, sph_diag+6) = ( sc%const_p_vjm1(1)/hl + sc%const_p_vjm1(2)*hlo_r) - (layer(il)%rho_i(ir)*layer(il)%g_i(ir)*sc%mass_bal_vjm1*hlo) * (1.0d0/layer(il)%K_i(ir)) S_U(smm+7*(ir-1)+2, sph_diag+7) = ( sc%const_p_vjp1(1)/hl + sc%const_p_vjp1(2)*hlo_r) - (layer(il)%rho_i(ir)*layer(il)%g_i(ir)*sc%mass_bal_vjp1*hlo) * (1.0d0/layer(il)%K_i(ir)) end if ! !CONSTITUTIVE EQUATION Y_{JM}^{J-2,2} S_U(smm+7*(ir-1)+5, sph_diag-4) = (-sc%Const_jm2_vjm1(1)/hl + sc%Const_jm2_vjm1(2)*hup_r) S_U(smm+7*(ir-1)+5, sph_diag+0) = 1.0d0 / (2.0d0 * layer(il)%mu_i(ir)*time_step) + 1.0d0/(2.0d0*layer(il)%eta_i(ir)) S_U(smm+7*(ir-1)+5, sph_diag+3) = ( sc%Const_jm2_vjm1(1)/hl + sc%Const_jm2_vjm1(2)*hlo_r) ! !CONSTITUTIVE EQUATION Y_{JM}^{J ,2} S_U(smm+7*(ir-1)+6, sph_diag-5) = (-sc%Const_j2__vjm1(1)/hl + sc%Const_j2__vjm1(2)*hup_r) S_U(smm+7*(ir-1)+6, sph_diag-4) = (-sc%Const_j2__vjp1(1)/hl + sc%Const_j2__vjp1(2)*hup_r) S_U(smm+7*(ir-1)+6, sph_diag+0) = 1.0d0 / (2.0d0 * layer(il)%mu_i(ir)*time_step) + 1.0d0/(2.0d0*layer(il)%eta_i(ir)) S_U(smm+7*(ir-1)+6, sph_diag+2) = ( sc%Const_j2__vjm1(1)/hl + sc%Const_j2__vjm1(2)*hlo_r) S_U(smm+7*(ir-1)+6, sph_diag+3) = ( sc%Const_j2__vjp1(1)/hl + sc%Const_j2__vjp1(2)*hlo_r) ! !CONSTITUTIVE EQUATION Y_{JM}^{J+2,2} S_U(smm+7*(ir-1)+7, sph_diag-5) = (-sc%Const_jp2_vjp1(1)/hl + sc%Const_jp2_vjp1(2)*hup_r) S_U(smm+7*(ir-1)+7, sph_diag+0) = 1.0d0 / (2.0d0 * layer(il)%mu_i(ir)*time_step) + 1.0d0/(2.0d0*layer(il)%eta_i(ir)) S_U(smm+7*(ir-1)+7, sph_diag+2) = ( sc%Const_jp2_vjp1(1)/hl + sc%Const_jp2_vjp1(2)*hlo_r) if(SWITCH_ocean_status == 'v')then ! !CONSTITUTIVE EQUATION Y_{JM}^{J-1,2} T_U(tmm+3*(ir-1) + 2, tor_diag-1) = (-tc%Const_jm1_vj(1)/hl + tc%Const_jm1_vj(2)*hup_r) T_U(tmm+3*(ir-1) + 2, tor_diag+0) = 1.0d0 / (2.0d0 * layer(il)%mu_i(ir)*time_step) + 1.0d0/(2.0d0*layer(il)%eta_i(ir)) T_U(tmm+3*(ir-1) + 2, tor_diag+2) = ( tc%Const_jm1_vj(1)/hl + tc%Const_jm1_vj(2)*hlo_r) ! !CONSTITUTIVE EQUATION Y_{JM}^{J+1,2} T_U(tmm+3*(ir-1) + 3, tor_diag-2) = (-tc%Const_jp1_vj(1)/hl + tc%Const_jp1_vj(2)*hup_r) T_U(tmm+3*(ir-1) + 3, tor_diag+0) = 1.0d0 / (2.0d0 * layer(il)%mu_i(ir)*time_step) + 1.0d0/(2.0d0*layer(il)%eta_i(ir)) T_U(tmm+3*(ir-1) + 3, tor_diag+1) = ( tc%Const_jp1_vj(1)/hl + tc%Const_jp1_vj(2)*hlo_r) end if end do smm=sum(layer(1:il)%sm); tmm=sum(layer(1:il)%tm) end do smm = 0; tmm=0 do il = 1, numoflayers nn = layer(il)%n do ir = 2, nn hi = layer(il)%r_i(ir) - layer(il)%r_i(ir-1) hup2 = (layer(il)%r_i(ir) - layer(il)%r_l( ir)) / hi; hup2_r = hup2/layer(il)%r_l(ir) hlo2 = (layer(il)%r_l(ir) - layer(il)%r_i(ir-1)) / hi; hlo2_r = hlo2/layer(il)%r_l(ir) !!Y_{JM}^{J-1} COMPONENT S_U(smm+7*(ir-1)+3, sph_diag-7) = (layer(il)%g_l(ir)) * sc%Mom__jm1_rho * hup2 ! !RHO_{JM} S_U(smm+7*(ir-1)+3, sph_diag-6) = (-sc%Mom__jm1_sj0(1) /hi + sc%Mom__jm1_sj0(2) * hup2_r) ! !S_{JM}^{J ,0} S_U(smm+7*(ir-1)+3, sph_diag-5) = (-sc%Mom__jm1_sjm2(1)/hi + sc%Mom__jm1_sjm2(2) * hup2_r) ! !S_{JM}^{J-2,2} S_U(smm+7*(ir-1)+3, sph_diag-4) = (-sc%Mom__jm1_sj2(1) /hi + sc%Mom__jm1_sj2(2) * hup2_r) ! !S_{JM}^{J ,2} S_U(smm+7*(ir-1)+3, sph_diag+0) = (layer(il)%g_l(ir)) * sc%Mom__jm1_rho * hlo2 ! !RHO_{JM} S_U(smm+7*(ir-1)+3, sph_diag+1) = ( sc%Mom__jm1_sj0(1) /hi + sc%Mom__jm1_sj0(2) * hlo2_r) ! !S_{JM}^{J ,0} S_U(smm+7*(ir-1)+3, sph_diag+2) = ( sc%Mom__jm1_sjm2(1)/hi + sc%Mom__jm1_sjm2(2) * hlo2_r) ! !S_{JM}^{J-2,2} S_U(smm+7*(ir-1)+3, sph_diag+3) = ( sc%Mom__jm1_sj2(1) /hi + sc%Mom__jm1_sj2(2) * hlo2_r) ! !S_{JM}^{J ,2} S_U(smm+7*(ir-1)+3, sph_diag-2) = -SWITCH_inertia*layer(il)%rho_l(ir)/time_step !!Y_{JM}^{J+1} COMPONENT S_U(smm+7*(ir-1)+4, sph_diag-8) = (layer(il)%g_l(ir)) * sc%Mom__jp1_rho*hup2 ! !RHO_{JM} S_U(smm+7*(ir-1)+4, sph_diag-7) = (-sc%Mom__jp1_sj0(1) /hi + sc%Mom__jp1_sj0(2) * hup2_r) ! !S_{JM}^{J ,0} S_U(smm+7*(ir-1)+4, sph_diag-5) = (-sc%Mom__jp1_sj2(1) /hi + sc%Mom__jp1_sj2(2) * hup2_r) ! !S_{JM}^{J ,2} S_U(smm+7*(ir-1)+4, sph_diag-4) = (-sc%Mom__jp1_sjp2(1)/hi + sc%Mom__jp1_sjp2(2) * hup2_r) ! !S_{JM}^{J+2,2} S_U(smm+7*(ir-1)+4, sph_diag-1) = (layer(il)%g_l(ir)) * sc%Mom__jp1_rho*hlo2 ! !RHO_{JM} S_U(smm+7*(ir-1)+4, sph_diag+0) = ( sc%Mom__jp1_sj0(1) /hi + sc%Mom__jp1_sj0(2) * hlo2_r) ! !S_{JM}^{J ,0} S_U(smm+7*(ir-1)+4, sph_diag+2) = ( sc%Mom__jp1_sj2(1) /hi + sc%Mom__jp1_sj2(2) * hlo2_r) ! !S_{JM}^{J ,2} S_U(smm+7*(ir-1)+4, sph_diag+3) = ( sc%Mom__jp1_sjp2(1)/hi + sc%Mom__jp1_sjp2(2) * hlo2_r) ! !S_{JM}^{J+2,2} S_U(smm+7*(ir-1)+4, sph_diag-2 ) = -SWITCH_inertia*layer(il)%rho_l(ir)/time_step if(SWITCH_ocean_status == 'v')then T_U(tmm+3*(ir-1)+1, tor_diag-2) = (-tc%Mom__j_s_jm1(1) /hi + tc%Mom__j_s_jm1(2) * hup2_r) T_U(tmm+3*(ir-1)+1, tor_diag-1) = (-tc%Mom__j_s_jp1(1) /hi + tc%Mom__j_s_jp1(2) * hup2_r) T_U(tmm+3*(ir-1)+1, tor_diag+0) = -SWITCH_inertia*layer(il)%rho_l(ir)/time_step T_U(tmm+3*(ir-1)+1, tor_diag+1) = ( tc%Mom__j_s_jm1(1) /hi + tc%Mom__j_s_jm1(2) * hlo2_r) T_U(tmm+3*(ir-1)+1, tor_diag+2) = ( tc%Mom__j_s_jp1(1) /hi + tc%Mom__j_s_jp1(2) * hlo2_r) end if end do smm = sum(layer(1:il)%sm); tmm=sum(layer(1:il)%tm) end do end subroutine generate_matrix_compressible subroutine internal_boundaries implicit none integer :: il, nn, smm, tmm real(8) :: hl, hup, hlo, hup_r, hlo_r, drho, g do il = 2, numoflayers nn = layer(il-1)%n; smm = sum(layer(1:il-1)%sm); tmm = sum(layer(1:il-1)%tm) drho = (layer(il-1)%rho_i(layer(il-1)%n) - layer(il)%rho_i(1)) g = layer(il-1)% g_i(layer(il-1)%n) hl = layer(il-1)%r_l(nn+1) - layer(il-1)%r_l(nn) hup = (layer(il-1)%r_l(nn+1) - layer(il-1)%r_i(nn))/(hl); hup_r = hup/layer(il-1)%r_i(nn) hlo = (layer(il-1)%r_i(nn) - layer(il-1)%r_l(nn))/(hl); hlo_r = hlo/layer(il-1)%r_i(nn) ! !CONTIUNITY OF TRACTION_{JM}^{J-1} S_U(smm+3, sph_diag-11) = -hup * sc_b%Bound_jm1_vjm1 * drho * g * 0.5d0 * time_step ! !V_{JM}^{J-1}, LOWER S_U(smm+3, sph_diag-10) = -hup * sc_b%Bound_jm1_vjp1 * drho * g * 0.5d0 * time_step ! !V_{JM}^{J-1}, LOWER S_U(smm+3, sph_diag-8 ) = -sc_b%Bound_jm1_s_0 ! !S_{JM}^{J ,0}, LOWER S_U(smm+3, sph_diag-7 ) = -sc_b%Bound_jm1_s_jm2 ! !S_{JM}^{J-2,2}, LOWER S_U(smm+3, sph_diag-6 ) = -sc_b%Bound_jm1_s_j2 ! !S_{JM}^{J ,2}, LOWER S_U(smm+3, sph_diag-4 ) = -hlo * sc_b%Bound_jm1_vjm1 * drho * g * 0.5d0 * time_step ! !V_{JM}^{J-1}, LOWER S_U(smm+3, sph_diag-3 ) = -hlo * sc_b%Bound_jm1_vjp1 * drho * g * 0.5d0 * time_step ! !V_{JM}^{J-1}, LOWER S_U(smm+3, sph_diag+1 ) = sc_b%Bound_jm1_s_0 ! !S_{JM}^{J ,0}, UPPER S_U(smm+3, sph_diag+2 ) = sc_b%Bound_jm1_s_jm2 ! !S_{JM}^{J-2,2}, UPPER S_U(smm+3, sph_diag+3 ) = sc_b%Bound_jm1_s_j2 ! !S_{JM}^{J ,2}, UPPER ! ?========================================================== ! ?========================================================== ! !CONTIUNITY OF TRACTION_{JM}^{J-1} S_U(smm+4, sph_diag-12) = -hup * sc_b%Bound_jp1_vjm1 * drho * g * 0.5d0 * time_step ! !V_{JM}^{J-1}, LOWER S_U(smm+4, sph_diag-11) = -hup * sc_b%Bound_jp1_vjp1 * drho * g * 0.5d0 * time_step ! !V_{JM}^{J+1}, LOWER S_U(smm+4, sph_diag-9 ) = -sc_b%Bound_jp1_s_0 ! !S_{JM}^{J ,0}, LOWER S_U(smm+4, sph_diag-7 ) = -sc_b%Bound_jp1_s_j2 ! !S_{JM}^{J-2,2}, LOWER S_U(smm+4, sph_diag-6 ) = -sc_b%Bound_jp1_s_jp2 ! !S_{JM}^{J ,2}, LOWER S_U(smm+4, sph_diag-5 ) = -hlo * sc_b%Bound_jp1_vjm1 * drho * g * 0.5d0 * time_step ! !V_{JM}^{J-1}, LOWER S_U(smm+4, sph_diag-4 ) = -hlo * sc_b%Bound_jp1_vjp1 * drho * g * 0.5d0 * time_step ! !V_{JM}^{J+1}, LOWER S_U(smm+4, sph_diag+0 ) = sc_b%Bound_jp1_s_0 ! !S_{JM}^{J ,0}, UPPER S_U(smm+4, sph_diag+2 ) = sc_b%Bound_jp1_s_j2 ! !S_{JM}^{J ,2}, UPPER S_U(smm+4, sph_diag+3 ) = sc_b%Bound_jp1_s_jp2 ! !S_{JM}^{J+2,2}, UPPER T_U(tmm+1, tor_diag-3) = -tc_b%Bound_j_s_jm1 !core T_U(tmm+1, tor_diag-2) = -tc_b%Bound_j_s_jp1 !core T_U(tmm+1, tor_diag+1) = tc_b%Bound_j_s_jm1!ocean T_U(tmm+1, tor_diag+2) = tc_b%Bound_j_s_jp1!ocean end do do il = 1, numoflayers-1 smm = sum(layer(1:il)%sm); tmm = sum(layer(1:il)%tm) ! !CONTIUNITY OF V_{JM}^{J-1} S_U(smm-1, sph_diag-7) = -0.5d0 ! !BELOW LAYER S_U(smm-1, sph_diag+0) = -0.5d0 ! !BELOW LAYER S_U(smm-1, sph_diag+2) = 0.5d0 ! !UPPER LAYER S_U(smm-1, sph_diag+9) = 0.5d0 ! !UPPER LAYER ! !CONTIUNITY OF V_{JM}^{J+1} S_U(smm , sph_diag-7) = -0.5d0 ! !BELOW LAYER S_U(smm , sph_diag+0) = -0.5d0 ! !BELOW LAYER S_U(smm , sph_diag+2) = 0.5d0 ! !UPPER LAYER S_U(smm , sph_diag+9) = 0.5d0 ! !UPPER LAYER if(SWITCH_ocean_status=='v')then T_U(tmm, tor_diag-3) = -0.5d0 T_U(tmm, tor_diag+0) = -0.5d0 T_U(tmm, tor_diag+1) = 0.5d0 T_U(tmm, tor_diag+4) = 0.5d0 end if end do end subroutine internal_boundaries subroutine top_bottom_boundaries implicit none integer :: il, nn real(8) :: hl, hup, hlo, rho, g T_U(1, tor_diag+0) = 1.0d0 S_U(3, sph_diag-2) = 1.0d0 S_U(4, sph_diag-2) = 1.0d0 il = numoflayers nn = layer(il)%n rho = layer(il)%rho_i(nn) g = layer(il)%g_i(nn) hl = layer(il)%r_l(nn+1) - layer(il)%r_l(nn) hup = (layer(il)%r_l(nn+1) - layer(il)%r_i(nn))/(hl) hlo = (layer(il)%r_i(nn) - layer(il)%r_l(nn))/(hl) !Free surface: Te_r = - u_r * rho * g = - (u_i + (v_{i+1} + v_{i})*dt/2) * rho * g !Traction vector j-1 S_U(smtot-1, sph_diag-7) = sc_b%Bound_jm1_vjm1 * hup * rho * g * 0.5d0 * time_step ! S_U(smtot-1, sph_diag-6) = sc_b%Bound_jm1_vjp1 * hup * rho * g * 0.5d0 * time_step ! S_U(smtot-1, sph_diag-4) = sc_b%Bound_jm1_s_0 ! !S_{JM}^{J ,0} S_U(smtot-1, sph_diag-3) = sc_b%Bound_jm1_s_jm2 ! !S_{JM}^{J-2,2} S_U(smtot-1, sph_diag-2) = sc_b%Bound_jm1_s_j2 ! !S_{JM}^{J ,2} S_U(smtot-1, sph_diag+0) = sc_b%Bound_jm1_vjm1* hlo * rho * g * 0.5d0 * time_step ! S_U(smtot-1, sph_diag+1) = sc_b%Bound_jm1_vjp1* hlo * rho * g * 0.5d0 * time_step ! !Traction vector j+1 S_U(smtot , sph_diag-8) = sc_b%Bound_jp1_vjm1* hup * rho * g * 0.5d0 * time_step ! S_U(smtot , sph_diag-7) = sc_b%Bound_jp1_vjp1* hup * rho * g * 0.5d0 * time_step ! S_U(smtot , sph_diag-5) = sc_b%Bound_jp1_s_0 ! !S_{JM}^{J ,0} S_U(smtot , sph_diag-3) = sc_b%Bound_jp1_s_j2 ! !S_{JM}^{J ,2} S_U(smtot , sph_diag-2) = sc_b%Bound_jp1_s_jp2 ! !S_{JM}^{J+2,2} S_U(smtot , sph_diag-1) = sc_b%Bound_jp1_vjm1 * hlo * rho * g * 0.5d0 * time_step ! S_U(smtot , sph_diag+0) = sc_b%Bound_jp1_vjp1 * hlo * rho * g * 0.5d0 * time_step ! T_U(tmtot, tor_diag-2) = tc_b%Bound_j_s_jm1 T_U(tmtot, tor_diag-1) = tc_b%Bound_j_s_jp1 end subroutine top_bottom_boundaries end module mod_matrix