! Calculations of forces from stresses subroutine fl_node include 'precision.inc' include 'params.inc' include 'arrays.inc' ! 1 - 3 ! | | ! 2 - 4 ! ! diagonal / : ! ! A: B: ! ! 1---3 1 ! | / / | ! 2 2---3 ! ! diagonal \ : ! ! C: D: ! ! 1 1---3 ! | \ \ | ! 2---3 2 ! ! assemblage of forces is COUNTRE CLOCK-WISE ! ! boff = 0 drat = dt / dt_elastic if (drat .lt. 1.) drat = 1. !!!$DIR LOOP_PARALLEL,LOOP_PRIVATE(j,fx,fy,p_est,rosubg,press_norm_l,dlx_l,dly_l,press_norm_r,dlx_r,dly_r),REDUCTION(boff) !$DIR PREFER_PARALLEL do i = 1,nx do j = 1,nz if(ynstressbc.eq.0.) then force(j,i,1) = 0 force(j,i,2) = 0 balance(j,i,1) = 0 balance(j,i,2) = 0 endif ! REGULAR PART - forces from stresses ! Element (j-1,i-1). Triangles B,C,D if ( j.ne.1 .and. i.ne.1 ) then ! triangle B ! side 2-3 fx = stress0(1,2,j-1,i-1) * (cord(j ,i ,2)-cord(j ,i-1,2)) - & stress0(3,2,j-1,i-1) * (cord(j ,i ,1)-cord(j ,i-1,1)) fy = stress0(3,2,j-1,i-1) * (cord(j ,i ,2)-cord(j ,i-1,2)) - & stress0(2,2,j-1,i-1) * (cord(j ,i ,1)-cord(j ,i-1,1)) force(j,i,1) = force(j,i,1) - 0.25*fx force(j,i,2) = force(j,i,2) - 0.25*fy balance(j,i,1) = balance(j,i,1) + 0.25*abs(fx) balance(j,i,2) = balance(j,i,2) + 0.25*abs(fy) ! side 3-1 fx = stress0(1,2,j-1,i-1) * (cord(j-1,i ,2)-cord(j ,i ,2)) - & stress0(3,2,j-1,i-1) * (cord(j-1,i ,1)-cord(j ,i ,1)) fy = stress0(3,2,j-1,i-1) * (cord(j-1,i ,2)-cord(j ,i ,2)) - & stress0(2,2,j-1,i-1) * (cord(j-1,i ,1)-cord(j ,i ,1)) force(j,i,1) = force(j,i,1) - 0.25*fx force(j,i,2) = force(j,i,2) - 0.25*fy balance(j,i,1) = balance(j,i,1) + 0.25*abs(fx) balance(j,i,2) = balance(j,i,2) + 0.25*abs(fy) ! triangle C ! side 2-3 fx = stress0(1,3,j-1,i-1) * (cord(j ,i ,2)-cord(j ,i-1,2)) - & stress0(3,3,j-1,i-1) * (cord(j ,i ,1)-cord(j ,i-1,1)) fy = stress0(3,3,j-1,i-1) * (cord(j ,i ,2)-cord(j ,i-1,2)) - & stress0(2,3,j-1,i-1) * (cord(j ,i ,1)-cord(j ,i-1,1)) force(j,i,1) = force(j,i,1) - 0.25*fx force(j,i,2) = force(j,i,2) - 0.25*fy balance(j,i,1) = balance(j,i,1) + 0.25*abs(fx) balance(j,i,2) = balance(j,i,2) + 0.25*abs(fy) ! side 3-1 fx = stress0(1,3,j-1,i-1) * (cord(j-1,i-1,2)-cord(j ,i ,2)) - & stress0(3,3,j-1,i-1) * (cord(j-1,i-1,1)-cord(j ,i ,1)) fy = stress0(3,3,j-1,i-1) * (cord(j-1,i-1,2)-cord(j ,i ,2)) - & stress0(2,3,j-1,i-1) * (cord(j-1,i-1,1)-cord(j ,i ,1)) force(j,i,1) = force(j,i,1) - 0.25*fx force(j,i,2) = force(j,i,2) - 0.25*fy balance(j,i,1) = balance(j,i,1) + 0.25*abs(fx) balance(j,i,2) = balance(j,i,2) + 0.25*abs(fy) ! triangle D ! side 1-2 fx = stress0(1,4,j-1,i-1) * (cord(j ,i ,2)-cord(j-1,i-1,2)) - & stress0(3,4,j-1,i-1) * (cord(j ,i ,1)-cord(j-1,i-1,1)) fy = stress0(3,4,j-1,i-1) * (cord(j ,i ,2)-cord(j-1,i-1,2)) - & stress0(2,4,j-1,i-1) * (cord(j ,i ,1)-cord(j-1,i-1,1)) force(j,i,1) = force(j,i,1) - 0.25*fx force(j,i,2) = force(j,i,2) - 0.25*fy balance(j,i,1) = balance(j,i,1) + 0.25*abs(fx) balance(j,i,2) = balance(j,i,2) + 0.25*abs(fy) ! side 2-3 fx = stress0(1,4,j-1,i-1) * (cord(j-1,i ,2)-cord(j ,i ,2)) - & stress0(3,4,j-1,i-1) * (cord(j-1,i ,1)-cord(j ,i ,1)) fy = stress0(3,4,j-1,i-1) * (cord(j-1,i ,2)-cord(j ,i ,2)) - & stress0(2,4,j-1,i-1) * (cord(j-1,i ,1)-cord(j ,i ,1)) force(j,i,1) = force(j,i,1) - 0.25*fx force(j,i,2) = force(j,i,2) - 0.25*fy balance(j,i,1) = balance(j,i,1) + 0.25*abs(fx) balance(j,i,2) = balance(j,i,2) + 0.25*abs(fy) endif ! Element (j-1,i). Triangles A,B,C. if ( j.ne.1 .and. i.ne.nx ) then ! triangle A ! side 1-2 fx = stress0(1,1,j-1,i ) * (cord(j ,i ,2)-cord(j-1,i ,2)) - & stress0(3,1,j-1,i ) * (cord(j ,i ,1)-cord(j-1,i ,1)) fy = stress0(3,1,j-1,i ) * (cord(j ,i ,2)-cord(j-1,i ,2)) - & stress0(2,1,j-1,i ) * (cord(j ,i ,1)-cord(j-1,i ,1)) force(j,i,1) = force(j,i,1) - 0.25*fx force(j,i,2) = force(j,i,2) - 0.25*fy balance(j,i,1) = balance(j,i,1) + 0.25*abs(fx) balance(j,i,2) = balance(j,i,2) + 0.25*abs(fy) ! side 2-3 fx = stress0(1,1,j-1,i ) * (cord(j-1,i+1,2)-cord(j ,i ,2)) - & stress0(3,1,j-1,i ) * (cord(j-1,i+1,1)-cord(j ,i ,1)) fy = stress0(3,1,j-1,i ) * (cord(j-1,i+1,2)-cord(j ,i ,2)) - & stress0(2,1,j-1,i ) * (cord(j-1,i+1,1)-cord(j ,i ,1)) force(j,i,1) = force(j,i,1) - 0.25*fx force(j,i,2) = force(j,i,2) - 0.25*fy balance(j,i,1) = balance(j,i,1) + 0.25*abs(fx) balance(j,i,2) = balance(j,i,2) + 0.25*abs(fy) ! triangle B ! side 1-2 fx = stress0(1,2,j-1,i ) * (cord(j ,i ,2)-cord(j-1,i+1,2)) - & stress0(3,2,j-1,i ) * (cord(j ,i ,1)-cord(j-1,i+1,1)) fy = stress0(3,2,j-1,i ) * (cord(j ,i ,2)-cord(j-1,i+1,2)) - & stress0(2,2,j-1,i ) * (cord(j ,i ,1)-cord(j-1,i+1,1)) force(j,i,1) = force(j,i,1) - 0.25*fx force(j,i,2) = force(j,i,2) - 0.25*fy balance(j,i,1) = balance(j,i,1) + 0.25*abs(fx) balance(j,i,2) = balance(j,i,2) + 0.25*abs(fy) ! side 2-3 fx = stress0(1,2,j-1,i ) * (cord(j ,i+1,2)-cord(j ,i ,2)) - & stress0(3,2,j-1,i ) * (cord(j ,i+1,1)-cord(j ,i ,1)) fy = stress0(3,2,j-1,i ) * (cord(j ,i+1,2)-cord(j ,i ,2)) - & stress0(2,2,j-1,i ) * (cord(j ,i+1,1)-cord(j ,i ,1)) force(j,i,1) = force(j,i,1) - 0.25*fx force(j,i,2) = force(j,i,2) - 0.25*fy balance(j,i,1) = balance(j,i,1) + 0.25*abs(fx) balance(j,i,2) = balance(j,i,2) + 0.25*abs(fy) ! triangle C ! side 1-2 fx = stress0(1,3,j-1,i ) * (cord(j ,i ,2)-cord(j-1,i ,2)) - & stress0(3,3,j-1,i ) * (cord(j ,i ,1)-cord(j-1,i ,1)) fy = stress0(3,3,j-1,i ) * (cord(j ,i ,2)-cord(j-1,i ,2)) - & stress0(2,3,j-1,i ) * (cord(j ,i ,1)-cord(j-1,i ,1)) force(j,i,1) = force(j,i,1) - 0.25*fx force(j,i,2) = force(j,i,2) - 0.25*fy balance(j,i,1) = balance(j,i,1) + 0.25*abs(fx) balance(j,i,2) = balance(j,i,2) + 0.25*abs(fy) ! side 2-3 fx = stress0(1,3,j-1,i ) * (cord(j ,i+1,2)-cord(j ,i ,2)) - & stress0(3,3,j-1,i ) * (cord(j ,i+1,1)-cord(j ,i ,1)) fy = stress0(3,3,j-1,i ) * (cord(j ,i+1,2)-cord(j ,i ,2)) - & stress0(2,3,j-1,i ) * (cord(j ,i+1,1)-cord(j ,i ,1)) force(j,i,1) = force(j,i,1) - 0.25*fx force(j,i,2) = force(j,i,2) - 0.25*fy balance(j,i,1) = balance(j,i,1) + 0.25*abs(fx) balance(j,i,2) = balance(j,i,2) + 0.25*abs(fy) endif ! Element (j,i-1). Triangles A,B,D if ( j.ne.nz .and. i.ne.1 ) then ! triangle A ! side 2-3 fx = stress0(1,1,j ,i-1) * (cord(j ,i ,2)-cord(j+1,i-1,2)) - & stress0(3,1,j ,i-1) * (cord(j ,i ,1)-cord(j+1,i-1,1)) fy = stress0(3,1,j ,i-1) * (cord(j ,i ,2)-cord(j+1,i-1,2)) - & stress0(2,1,j ,i-1) * (cord(j ,i ,1)-cord(j+1,i-1,1)) force(j,i,1) = force(j,i,1) - 0.25*fx force(j,i,2) = force(j,i,2) - 0.25*fy balance(j,i,1) = balance(j,i,1) + 0.25*abs(fx) balance(j,i,2) = balance(j,i,2) + 0.25*abs(fy) ! side 3-1 fx = stress0(1,1,j ,i-1) * (cord(j ,i-1,2)-cord(j ,i ,2)) - & stress0(3,1,j ,i-1) * (cord(j ,i-1,1)-cord(j ,i ,1)) fy = stress0(3,1,j ,i-1) * (cord(j ,i-1,2)-cord(j ,i ,2)) - & stress0(2,1,j ,i-1) * (cord(j ,i-1,1)-cord(j ,i ,1)) force(j,i,1) = force(j,i,1) - 0.25*fx force(j,i,2) = force(j,i,2) - 0.25*fy balance(j,i,1) = balance(j,i,1) + 0.25*abs(fx) balance(j,i,2) = balance(j,i,2) + 0.25*abs(fy) ! triangle B ! side 1-2 fx = stress0(1,2,j ,i-1) * (cord(j+1,i-1,2)-cord(j ,i ,2)) - & stress0(3,2,j ,i-1) * (cord(j+1,i-1,1)-cord(j ,i ,1)) fy = stress0(3,2,j ,i-1) * (cord(j+1,i-1,2)-cord(j ,i ,2)) - & stress0(2,2,j ,i-1) * (cord(j+1,i-1,1)-cord(j ,i ,1)) force(j,i,1) = force(j,i,1) - 0.25*fx force(j,i,2) = force(j,i,2) - 0.25*fy balance(j,i,1) = balance(j,i,1) + 0.25*abs(fx) balance(j,i,2) = balance(j,i,2) + 0.25*abs(fy) ! side 3-1 fx = stress0(1,2,j ,i-1) * (cord(j ,i ,2)-cord(j+1,i ,2)) - & stress0(3,2,j ,i-1) * (cord(j ,i ,1)-cord(j+1,i ,1)) fy = stress0(3,2,j ,i-1) * (cord(j ,i ,2)-cord(j+1,i ,2)) - & stress0(2,2,j ,i-1) * (cord(j ,i ,1)-cord(j+1,i ,1)) force(j,i,1) = force(j,i,1) - 0.25*fx force(j,i,2) = force(j,i,2) - 0.25*fy balance(j,i,1) = balance(j,i,1) + 0.25*abs(fx) balance(j,i,2) = balance(j,i,2) + 0.25*abs(fy) ! triangle D ! side 2-3 fx = stress0(1,4,j ,i-1) * (cord(j ,i ,2)-cord(j+1,i ,2)) - & stress0(3,4,j ,i-1) * (cord(j ,i ,1)-cord(j+1,i ,1)) fy = stress0(3,4,j ,i-1) * (cord(j ,i ,2)-cord(j+1,i ,2)) - & stress0(2,4,j ,i-1) * (cord(j ,i ,1)-cord(j+1,i ,1)) force(j,i,1) = force(j,i,1) - 0.25*fx force(j,i,2) = force(j,i,2) - 0.25*fy balance(j,i,1) = balance(j,i,1) + 0.25*abs(fx) balance(j,i,2) = balance(j,i,2) + 0.25*abs(fy) ! side 3-1 fx = stress0(1,4,j ,i-1) * (cord(j ,i-1,2)-cord(j ,i ,2)) - & stress0(3,4,j ,i-1) * (cord(j ,i-1,1)-cord(j ,i ,1)) fy = stress0(3,4,j ,i-1) * (cord(j ,i-1,2)-cord(j ,i ,2)) - & stress0(2,4,j ,i-1) * (cord(j ,i-1,1)-cord(j ,i ,1)) force(j,i,1) = force(j,i,1) - 0.25*fx force(j,i,2) = force(j,i,2) - 0.25*fy balance(j,i,1) = balance(j,i,1) + 0.25*abs(fx) balance(j,i,2) = balance(j,i,2) + 0.25*abs(fy) endif ! Element (j,i). Triangles A,C,D if ( j.ne.nz .and. i.ne.nx ) then ! triangle A ! side 1-2 fx = stress0(1,1,j ,i ) * (cord(j+1,i ,2)-cord(j ,i ,2)) - & stress0(3,1,j ,i ) * (cord(j+1,i ,1)-cord(j ,i ,1)) fy = stress0(3,1,j ,i ) * (cord(j+1,i ,2)-cord(j ,i ,2)) - & stress0(2,1,j ,i ) * (cord(j+1,i ,1)-cord(j ,i ,1)) force(j,i,1) = force(j,i,1) - 0.25*fx force(j,i,2) = force(j,i,2) - 0.25*fy balance(j,i,1) = balance(j,i,1) + 0.25*abs(fx) balance(j,i,2) = balance(j,i,2) + 0.25*abs(fy) ! side 3-1 fx = stress0(1,1,j ,i ) * (cord(j ,i ,2)-cord(j ,i+1,2)) - & stress0(3,1,j ,i ) * (cord(j ,i ,1)-cord(j ,i+1,1)) fy = stress0(3,1,j ,i ) * (cord(j ,i ,2)-cord(j ,i+1,2)) - & stress0(2,1,j ,i ) * (cord(j ,i ,1)-cord(j ,i+1,1)) force(j,i,1) = force(j,i,1) - 0.25*fx force(j,i,2) = force(j,i,2) - 0.25*fy balance(j,i,1) = balance(j,i,1) + 0.25*abs(fx) balance(j,i,2) = balance(j,i,2) + 0.25*abs(fy) ! triangle C ! side 1-2 fx = stress0(1,3,j ,i ) * (cord(j+1,i ,2)-cord(j ,i ,2)) - & stress0(3,3,j ,i ) * (cord(j+1,i ,1)-cord(j ,i ,1)) fy = stress0(3,3,j ,i ) * (cord(j+1,i ,2)-cord(j ,i ,2)) - & stress0(2,3,j ,i ) * (cord(j+1,i ,1)-cord(j ,i ,1)) force(j,i,1) = force(j,i,1) - 0.25*fx force(j,i,2) = force(j,i,2) - 0.25*fy balance(j,i,1) = balance(j,i,1) + 0.25*abs(fx) balance(j,i,2) = balance(j,i,2) + 0.25*abs(fy) ! side 3-1 fx = stress0(1,3,j ,i ) * (cord(j ,i ,2)-cord(j+1,i+1,2)) - & stress0(3,3,j ,i ) * (cord(j ,i ,1)-cord(j+1,i+1,1)) fy = stress0(3,3,j ,i ) * (cord(j ,i ,2)-cord(j+1,i+1,2)) - & stress0(2,3,j ,i ) * (cord(j ,i ,1)-cord(j+1,i+1,1)) force(j,i,1) = force(j,i,1) - 0.25*fx force(j,i,2) = force(j,i,2) - 0.25*fy balance(j,i,1) = balance(j,i,1) + 0.25*abs(fx) balance(j,i,2) = balance(j,i,2) + 0.25*abs(fy) ! triangle D ! side 1-2 fx = stress0(1,4,j ,i ) * (cord(j+1,i+1,2)-cord(j ,i ,2)) - & stress0(3,4,j ,i ) * (cord(j+1,i+1,1)-cord(j ,i ,1)) fy = stress0(3,4,j ,i ) * (cord(j+1,i+1,2)-cord(j ,i ,2)) - & stress0(2,4,j ,i ) * (cord(j+1,i+1,1)-cord(j ,i ,1)) force(j,i,1) = force(j,i,1) - 0.25*fx force(j,i,2) = force(j,i,2) - 0.25*fy balance(j,i,1) = balance(j,i,1) + 0.25*abs(fx) balance(j,i,2) = balance(j,i,2) + 0.25*abs(fy) ! side 3-1 fx = stress0(1,4,j ,i ) * (cord(j ,i ,2)-cord(j ,i+1,2)) - & stress0(3,4,j ,i ) * (cord(j ,i ,1)-cord(j ,i+1,1)) fy = stress0(3,4,j ,i ) * (cord(j ,i ,2)-cord(j ,i+1,2)) - & stress0(2,4,j ,i ) * (cord(j ,i ,1)-cord(j ,i+1,1)) force(j,i,1) = force(j,i,1) - 0.25*fx force(j,i,2) = force(j,i,2) - 0.25*fy balance(j,i,1) = balance(j,i,1) + 0.25*abs(fx) balance(j,i,2) = balance(j,i,2) + 0.25*abs(fy) endif ! GRAVITY FORCE force(j,i,2) = force(j,i,2) - rmass(j,i)*g balance(j,i,2) = balance(j,i,2) + abs(rmass(j,i)*g) ! BOUNDARY CONDITIONS ! bottom support - Archimed force (normal to the surface, shear component = 0) if(nyhydro.gt.0 .and. j.eq.nz) then p_est = pisos + 0.5*(den(iphsub)+drosub)*g*(cord(nz,i,2)-rzbo) rosubg = g * (den(iphsub)+drosub) * (1-alfa(iphsub)*temp(nz,i)+beta(iphsub)*p_est) if(i.eq.1) then press_norm_l = 0 dlx_l = 0 dly_l = 0 press_norm_r = pisos-rosubg*((cord(nz,i+1,2)+cord(nz,i,2))/2-rzbo) dlx_r = cord(nz,i+1,1)-cord(nz,i ,1) dly_r = cord(nz,i+1,2)-cord(nz,i ,2) elseif(i.eq.nx) then press_norm_l = pisos-rosubg*((cord(nz,i-1,2)+cord(nz,i,2))/2-rzbo) dlx_l = cord(nz,i ,1)-cord(nz,i-1,1) dly_l = cord(nz,i ,2)-cord(nz,i-1,2) press_norm_r = 0 dlx_r = 0 dly_r = 0 else press_norm_l = pisos-rosubg*((cord(nz,i-1,2)+cord(nz,i,2))/2-rzbo) dlx_l = cord(nz,i ,1)-cord(nz,i-1,1) dly_l = cord(nz,i ,2)-cord(nz,i-1,2) press_norm_r = pisos-rosubg*((cord(nz,i+1,2)+cord(nz,i,2))/2-rzbo) dlx_r = cord(nz,i+1,1)-cord(nz,i ,1) dly_r = cord(nz,i+1,2)-cord(nz,i ,2) endif force(nz,i,1) = force(nz,i,1)-0.5*press_norm_l*dly_l-0.5*press_norm_r*dly_r force(nz,i,2) = force(nz,i,2)+0.5*press_norm_l*dlx_l+0.5*press_norm_r*dlx_r balance(nz,i,1) = 1.e+17 endif ! BALANCE-OFF if( iand(ncod(j,i,1),1).eq.1 .or. j.le.n_boff_cutoff ) then balance(j,i,1) = 0 else balance(j,i,1) = abs(force(j,i,1)) / (balance(j,i,1) + 1.e-9) endif if( iand(ncod(j,i,2),2).eq.2 .or. j.le.n_boff_cutoff ) then balance(j,i,2) = 0 else balance(j,i,2) = abs(force(j,i,2)) / (balance(j,i,2) + 1.e-9) endif ! DAMPING if( iand(ncod(j,i,1),1).ne.1 .and. abs(vel(j,i,1)).gt.1.e-13 ) then force(j,i,1) = force(j,i,1) - demf*sign(force(j,i,1),vel(j,i,1)) endif if( iand(ncod(j,i,2),2).ne.2 .and. abs(vel(j,i,2)).gt.1.e-13 ) then force(j,i,2) = force(j,i,2) - demf*sign(force(j,i,2),vel(j,i,2)) endif ! VELOCITIES FROM FORCES if( ncod(j,i,1) .eq. 1 ) then vel(j,i,1) = bc(j,i,1) ! if(j.eq.30) vel(j,i,1) =0. ! write(*,*) i,j,vel(j,i,1) else vel(j,i,1) = vel(j,i,1) + dt*force(j,i,1)/(amass(j,i)*drat*drat) endif if( ncod(j,i,2) .eq. 1 ) then vel(j,i,2) = bc(j,i,2) ! write(*,*) i,j,vel(j,i,2) else vel(j,i,2) = vel(j,i,2) + dt*force(j,i,2)/(amass(j,i)*drat*drat) endif ! MAX balance-off boff = max(boff,balance(j,i,1)) boff = max(boff,balance(j,i,2)) end do end do return end