module mod_rhs use mod_parameters use mod_variables use mod_potential use mod_math implicit none contains subroutine generate_rhs(t2,j,m) implicit none integer :: j, m, jm1, jp1, jm0, jms, smm, tmm integer :: il, ir, nn real(8) :: dt_inertia, t2, jj, erjm1, erjp1, dtmu real(8) :: uc1, uc2, vc1, vc2, vc3, rg, dt2 complex(8) :: ur, avjm1, avjp1 jm1 = jml2(j,m,j-1) jm0 = jml2(j,m,j ) jp1 = jml2(j,m,j+1) jms = jm2(j,m)-1 jj = dble(j) erjm1 = dsqrt( jj / (2.0d0 * jj + 1.0d0)) erjp1 = -dsqrt( (jj+1.0d0) / (2.0d0 * jj + 1.0d0)) uc1 = -erjm1; uc2 = -erjp1; vc1 = jj / (2.0d0 * jj + 1.0d0) vc2 = (jj+1.0d0) / (2.0d0 * jj + 1.0d0) vc3 = dsqrt(jj*(jj+1.0d0)) / (2.0d0 * jj + 1.0d0) dt_inertia = -SWITCH_inertia/time_step dt2 = 0.5d0 * time_step smm = 0; tmm =0; RHS%spheriodal(:) = dzero RHS%torodial (:) = dzero if(abs(CN_theta - 1.0d0)>1d-5) call Divergence_tensor do il = 1, numoflayers nn = layer(il)%n do ir = 1, nn ! ?===================================================================================================================================================================== ! !======================================BALANCE OF MASS=================================================================================== if(il/=layer_ocean) RHS%spheriodal(smm+7*(ir-1)+1) = layer(il)%rhoe(ir, jms+1)/time_step ! ?===================================================================================================================================================================== ! ?===================================================================================================================================================================== ! !======================================PRESSURE PART OF THE CONSTITUTIVE EQUATION=================================================================================== if(SWITCH_compressible==1)then if(il == layer_ocean) then RHS%spheriodal(smm+7*(ir-1)+2) = dzero else RHS%spheriodal(smm+7*(ir-1)+2) = pjmstr*layer(il)%str(ir, 6*jms+1)/(layer(il)%K_i(ir)*time_step) end if end if ! ?===================================================================================================================================================================== ! ?===================================================================================================================================================================== ! !======================================DEVIATORIC PART OF THE CONSTITUTIVE EQUATION=================================================================================== dtmu = 1.0d0/(2.0d0*layer(il)%mu_i(ir)*time_step) RHS%spheriodal(smm+7*(ir-1)+5) = layer(il)%str (ir, 6*jms+2)*dtmu ! !DEVIATORIC PART OF THE STRESS TENSOR TIME DERIVATIVE RHS%spheriodal(smm+7*(ir-1)+6) = layer(il)%str (ir, 6*jms+4)*dtmu ! !DEVIATORIC PART OF THE STRESS TENSOR TIME DERIVATIVE RHS%spheriodal(smm+7*(ir-1)+7) = layer(il)%str (ir, 6*jms+6)*dtmu ! !DEVIATORIC PART OF THE STRESS TENSOR TIME DERIVATIVE RHS%torodial(tmm+3*(ir-1)+2) = layer(il)%str (ir, 6*jms+3)*dtmu RHS%torodial(tmm+3*(ir-1)+3) = layer(il)%str (ir, 6*jms+5)*dtmu ! ?===================================================================================================================================================================== if(ir>1)then RHS%spheriodal(smm+7*(ir-1)+3) = dt_inertia * layer(il)%vel(ir,jm1) * layer(il)%rho_l(ir) - (1.0d0-CN_theta)*layer(il)%div_tensor(ir-1, jm1) ! !MOMENTUM EQUATION TIME DERIVATIVE OF VELOCITY V_{JM}^{J-1} RHS%spheriodal(smm+7*(ir-1)+4) = dt_inertia * layer(il)%vel(ir,jp1) * layer(il)%rho_l(ir) - (1.0d0-CN_theta)*layer(il)%div_tensor(ir-1, jp1) ! !MOMENTUM EQUATION TIME DERIVATIVE OF VELOCITY V_{JM}^{J+1} RHS%torodial (tmm+3*(ir-1)+1) = dt_inertia * layer(il)%vel(ir,jm0) * layer(il)%rho_l(ir) - (1.0d0-CN_theta)*layer(il)%div_tensor(ir-1, jm0) ! !MOMENTUM EQUATION TIME DERIVATIVE OF VELOCITY V_{JM}^{J} if(SWITCH_self_grav == 1 )then RHS%spheriodal(smm+7*(ir-1)+3) = RHS%spheriodal(smm+7*(ir-1)+3) - layer(il)%rho_l(ir) * layer(il)%grad_pot(ir-1, jm1) * weight_sg RHS%spheriodal(smm+7*(ir-1)+4) = RHS%spheriodal(smm+7*(ir-1)+4) - layer(il)%rho_l(ir) * layer(il)%grad_pot(ir-1, jp1) * weight_sg end if if((SWITCH_coriolis == 1).and.(il == layer_ocean)) then RHS%spheriodal(smm+7*(ir-1)+3) = RHS%spheriodal (smm+7*(ir-1)+3) - layer(il)%fc(ir-1, jm1) RHS%spheriodal(smm+7*(ir-1)+4) = RHS%spheriodal (smm+7*(ir-1)+4) - layer(il)%fc(ir-1, jp1) RHS%torodial (tmm+3*(ir-1)+1) = RHS%torodial (tmm+3*(ir-1)+1) - layer(il)%fc(ir-1, jm0) end if if(j==2)then RHS%spheriodal(smm+7*(ir-1)+3) = RHS%spheriodal(smm+7*(ir-1)+3) + tides(layer(il)%rho_l(ir), layer(il)%r_l(ir), m, t2) if(SWITCH_mmtides==1) RHS%spheriodal(smm+7*(ir-1)+3) = RHS%spheriodal(smm+7*(ir-1)+3) + mm_tides(layer(il)%rho_l(ir), layer(il)%r_l(ir), m, t2) end if end if end do smm = sum(layer(1:il)%sm); tmm = sum(layer(1:il)%tm) end do do il = 1, numoflayers-1 smm = sum(layer(1:il)%sm); tmm = sum(layer(1:il)%tm) nn = layer(il)%n ur = erjm1 * layer(il)%dis(nn, jm1) + erjp1 * layer(il)%dis(nn, jp1) avjm1 = 0.5d0 * (layer(il)%vel(nn, jm1) + layer(il)%vel(nn+1, jm1)) avjp1 = 0.5d0 * (layer(il)%vel(nn, jp1) + layer(il)%vel(nn+1, jp1)) ! !INTERFACE CONDITIONS (INTEGRATION OF THE RADIAL VELOCITY) rg = (layer(il)%rho_i(nn) - layer(il+1)%rho_i(1)) * layer(il)%g_i(nn) RHS%spheriodal(smm+3) = (-rg * (uc1 * ur - dt2 * vc1 * avjm1 + dt2 * vc3 * avjp1)) RHS%spheriodal(smm+4) = (-rg * (uc2 * ur + dt2 * vc3 * avjm1 - dt2 * vc2 * avjp1)) end do il = numoflayers; nn = layer(il)%n ur = layer(il)%ur_i(nn, jm2(j,m)) avjm1 = 0.5d0 * (layer(il)%vel(nn, jm1) + layer(il)%vel(nn+1, jm1)) avjp1 = 0.5d0 * (layer(il)%vel(nn, jp1) + layer(il)%vel(nn+1, jp1)) rg = (layer(numoflayers)%rho_i(nn)) * layer(numoflayers)%g_i(nn) RHS%spheriodal(smtot-1) = ( rg * (uc1 * ur - dt2 * vc1 * avjm1 + dt2 * vc3 * avjp1)) RHS%spheriodal(smtot ) = ( rg * (uc2 * ur + dt2 * vc3 * avjm1 - dt2 * vc2 * avjp1)) end subroutine generate_rhs end module mod_rhs