XBeach
|
00001 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00002 ! Copyright (C) 2011 UNESCO-IHE, WL|Delft Hydraulics and Delft University ! 00003 ! Dano Roelvink, Ap van Dongeren, Ad Reniers, Jamie Lescinski, ! 00004 ! Jaap van Thiel de Vries, Robert McCall ! 00005 ! ! 00006 ! d.roelvink@unesco-ihe.org ! 00007 ! UNESCO-IHE Institute for Water Education ! 00008 ! P.O. Box 3015 ! 00009 ! 2601 DA Delft ! 00010 ! The Netherlands ! 00011 ! ! 00012 ! This library is free software; you can redistribute it and/or ! 00013 ! modify it under the terms of the GNU Lesser General Public ! 00014 ! License as published by the Free Software Foundation; either ! 00015 ! version 2.1 of the License, or (at your option) any later version. ! 00016 ! ! 00017 ! This library is distributed in the hope that it will be useful, ! 00018 ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! 00019 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ! 00020 ! Lesser General Public License for more details. ! 00021 ! ! 00022 ! You should have received a copy of the GNU Lesser General Public ! 00023 ! License along with this library; if not, write to the Free Software ! 00024 ! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 ! 00025 ! USA ! 00026 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00027 module ship_module 00028 use typesandkinds, only: slen 00029 implicit none 00030 save 00031 type ship 00032 character(slen) :: name 00033 real*8 :: dx 00034 real*8 :: dy 00035 integer :: nx 00036 integer :: ny 00037 integer :: compute_force ! option (0/1) to compute forces on ship (can be switched on/off per ship) 00038 integer :: compute_motion ! option (0/1) to compute ship motion due to waves 00039 integer :: flying 00040 integer :: read_heading 00041 character(slen) :: shipgeom 00042 real*8 :: xCG ! x location of center of gravity w.r.t. ship grid 00043 real*8 :: yCG ! y location of center of gravity w.r.t. ship grid 00044 real*8 :: zCG ! z location of center of gravity w.r.t. ship grid 00045 real*8, dimension(:,:), pointer :: x 00046 real*8, dimension(:,:), pointer :: y 00047 real*8, dimension(:,:), pointer :: depth ! ship depth defined on local grid relative to plane fixed to ship 00048 real*8, dimension(:,:), pointer :: zhull ! actual z-coordinate of ship hull relative to horizontal reference plane 00049 real*8, dimension(:,:), pointer :: zs ! water level defined on ship grid, interpolated from XBeach 00050 real*8, dimension(:,:), pointer :: ph ! pressure head at ship hull = zs-zhull 00051 character(slen) :: shiptrack 00052 integer :: track_nt 00053 real*8 , dimension(:) , pointer :: track_t 00054 real*8 , dimension(:) , pointer :: track_x 00055 real*8 , dimension(:) , pointer :: track_y 00056 real*8 , dimension(:) , pointer :: track_z 00057 real*8 , dimension(:) , pointer :: track_dir 00058 real*8 :: mass 00059 real*8 :: Jx 00060 real*8 :: Jy 00061 real*8 :: Jz 00062 real*8 , dimension(:) , pointer :: xs 00063 real*8 , dimension(:) , pointer :: ys 00064 integer, dimension(:) , pointer :: nrx 00065 integer, dimension(:) , pointer :: nry 00066 integer, dimension(:) , pointer :: nrin 00067 integer, dimension(:) , pointer :: iflag 00068 integer, dimension(:,:), pointer :: iref 00069 real*8 , dimension(:,:), pointer :: w 00070 integer, dimension(:) , pointer :: covered 00071 end type ship 00072 00073 contains 00074 00075 subroutine ship_init(s,par,sh) 00076 use params, only: parameters 00077 use xmpi_module 00078 use spaceparams 00079 use readkey_module, only: count_lines, readkey_dbl, readkey_int, readkey_intvec, readkey_name, read_v 00080 use filefunctions, only: create_new_fid 00081 use interp 00082 00083 implicit none 00084 00085 type(parameters) :: par 00086 type(spacepars),target :: s 00087 type(ship), dimension(:), pointer :: sh 00088 00089 integer :: i,fid,iy,it 00090 integer :: n2 00091 logical :: toall = .true. 00092 00093 !include 's.ind' 00094 !include 's.inp' 00095 if(par%ships==1) then 00096 ! Read ship names (== filenames with ship geometry and track data) 00097 00098 par%nship = count_lines(par%shipfile) 00099 00100 allocate(sh(par%nship)) 00101 if (xmaster) then 00102 fid=create_new_fid() 00103 open(fid,file=par%shipfile) 00104 do i=1,par%nship 00105 read(fid,'(a)') sh(i)%name 00106 enddo 00107 close(fid) 00108 endif 00109 #ifdef USEMPI 00110 do i=1,par%nship 00111 call xmpi_bcast(sh(i)%name,toall) 00112 enddo 00113 #endif 00114 do i=1,par%nship 00115 ! Read ship geometry 00116 sh(i)%dx = readkey_dbl(sh(i)%name,'dx', 5.d0, 0.d0, 100.d0) 00117 sh(i)%dy = readkey_dbl(sh(i)%name,'dy', 5.d0, 0.d0, 100.d0) 00118 sh(i)%nx = readkey_int(sh(i)%name,'nx', 20, 1, 1000 ) 00119 sh(i)%ny = readkey_int(sh(i)%name,'ny', 20, 1, 1000 ) 00120 sh(i)%shipgeom = readkey_name(sh(i)%name,'shipgeom',required=.true.) 00121 sh(i)%xCG = readkey_dbl(sh(i)%name,'xCG', 0.d0, -1000.d0, 1000.d0) 00122 sh(i)%yCG = readkey_dbl(sh(i)%name,'yCG', 0.d0, -1000.d0, 1000.d0) 00123 sh(i)%zCG = readkey_dbl(sh(i)%name,'zCG', 0.d0, -1000.d0, 1000.d0) 00124 sh(i)%shiptrack = readkey_name(sh(i)%name,'shiptrack',required=.true.) 00125 sh(i)%compute_force = readkey_int(sh(i)%name,'compute_force' , 0, 0, 1) 00126 sh(i)%compute_motion = readkey_int(sh(i)%name,'compute_motion', 0, 0, 1) 00127 sh(i)%flying = readkey_int(sh(i)%name,'flying', 0, 0, 1) 00128 sh(i)%read_heading = readkey_int(sh(i)%name,'read_heading', 0, 0, 1) 00129 00130 allocate (sh(i)%depth(sh(i)%nx+1,sh(i)%ny+1)) 00131 allocate (sh(i)%zhull(sh(i)%nx+1,sh(i)%ny+1)) 00132 allocate (sh(i)%ph(sh(i)%nx+1,sh(i)%ny+1)) 00133 allocate (sh(i)%zs(sh(i)%nx+1,sh(i)%ny+1)) 00134 allocate (sh(i)%x(sh(i)%nx+1,sh(i)%ny+1)) 00135 allocate (sh(i)%y(sh(i)%nx+1,sh(i)%ny+1)) 00136 n2=(sh(i)%nx+1)*(sh(i)%ny+1) 00137 allocate (sh(i)%xs(n2)) 00138 allocate (sh(i)%ys(n2)) 00139 allocate (sh(i)%nrx(n2)) 00140 allocate (sh(i)%nry(n2)) 00141 allocate (sh(i)%nrin(n2)) 00142 allocate (sh(i)%iflag(n2)) 00143 allocate (sh(i)%iref(4,n2)) 00144 allocate (sh(i)%w(4,n2)) 00145 allocate (sh(i)%covered(n2)) 00146 sh(i)%ph = 0.d0 00147 sh(i)%zs = 0.d0 00148 00149 fid=create_new_fid() 00150 !open(fid,file=sh(i)%shipgeom) 00151 if(xmaster) open(fid,file=sh(i)%shipgeom) 00152 00153 do iy=1,sh(i)%ny+1 00154 call read_v(fid,sh(i)%depth(:,iy)) 00155 enddo 00156 if(xmaster) close(fid) 00157 00158 if (sh(i)%compute_motion==0) then 00159 sh(i)%zhull=-sh(i)%depth 00160 endif 00161 if (sh(i)%compute_force==0) then 00162 sh(i)%ph=sh(i)%depth 00163 endif 00164 00165 ! Read t,x,y of ship position 00166 00167 sh(i)%track_nt=count_lines(sh(i)%shiptrack) 00168 00169 fid=create_new_fid() 00170 if(xmaster) open(fid,file=sh(i)%shiptrack) 00171 00172 allocate(sh(i)%track_t(sh(i)%track_nt)) 00173 allocate(sh(i)%track_x(sh(i)%track_nt)) 00174 allocate(sh(i)%track_y(sh(i)%track_nt)) 00175 allocate(sh(i)%track_z(sh(i)%track_nt)) 00176 allocate(sh(i)%track_dir(sh(i)%track_nt)) 00177 00178 if (sh(i)%flying==0) then 00179 if (sh(i)%read_heading==0) then 00180 do it=1,sh(i)%track_nt 00181 call read_v(fid,sh(i)%track_t(it),sh(i)%track_x(it),sh(i)%track_y(it)) 00182 enddo 00183 else 00184 do it=1,sh(i)%track_nt 00185 ! also read in track_dir, in nautical deg 00186 call read_v(fid,sh(i)%track_t(it),sh(i)%track_x(it),sh(i)%track_y(it),sh(i)%track_dir(it)) 00187 ! convert to cartesian direction 00188 sh(i)%track_dir(it)=par%px/180d0*(270d0-sh(i)%track_dir(it)) 00189 enddo 00190 endif 00191 sh(i)%track_z=0.d0 00192 else 00193 if (sh(i)%read_heading==0) then 00194 do it=1,sh(i)%track_nt 00195 call read_v(fid,sh(i)%track_t(it),sh(i)%track_x(it),sh(i)%track_y(it),sh(i)%track_z(it)) 00196 enddo 00197 else 00198 do it=1,sh(i)%track_nt 00199 ! also read in track_dir, in nautical deg 00200 call read_v(fid,sh(i)%track_t(it),sh(i)%track_x(it),sh(i)%track_y(it),sh(i)%track_dir(it),sh(i)%track_z(it)) 00201 ! convert to cartesian direction 00202 sh(i)%track_dir(it)=par%px/180d0*(270d0-sh(i)%track_dir(it)) 00203 enddo 00204 endif 00205 endif 00206 if(xmaster)close(fid) 00207 00208 ! Compute ship direction 00209 00210 if (sh(i)%read_heading==0) then 00211 sh(i)%track_dir(1)=atan2(sh(i)%track_y(2)-sh(i)%track_y(1),sh(i)%track_x(2)-sh(i)%track_x(1)) 00212 do it=2,sh(i)%track_nt-1 00213 sh(i)%track_dir(it)=atan2(sh(i)%track_y(it+1)-sh(i)%track_y(it-1),sh(i)%track_x(it+1)-sh(i)%track_x(it-1)) 00214 enddo 00215 it=sh(i)%track_nt 00216 sh(i)%track_dir(it)=atan2(sh(i)%track_y(it)-sh(i)%track_y(it-1),sh(i)%track_x(it)-sh(i)%track_x(it-1)) 00217 endif 00218 enddo ! loop over ships 00219 00220 ! Initialize ship-related global variables 00221 if (xmaster) then 00222 ! only on xmaster, rest is done automatically by call from libxbeach 00223 allocate(s%shipxCG (par%nship)) 00224 allocate(s%shipyCG (par%nship)) 00225 allocate(s%shipzCG (par%nship)) 00226 allocate(s%shipFx (par%nship)) 00227 allocate(s%shipFy (par%nship)) 00228 allocate(s%shipFz (par%nship)) 00229 allocate(s%shipMx (par%nship)) 00230 allocate(s%shipMy (par%nship)) 00231 allocate(s%shipMz (par%nship)) 00232 allocate(s%shipphi (par%nship)) 00233 allocate(s%shipchi (par%nship)) 00234 allocate(s%shippsi (par%nship)) 00235 s%shipxCG=0.d0 00236 s%shipyCG=0.d0 00237 s%shipzCG=0.d0 00238 s%shipFx=0.d0 00239 s%shipFy=0.d0 00240 s%shipFz=0.d0 00241 endif 00242 else ! (par%ships==0) 00243 if (xmaster) then 00244 ! just allocate address for memory, only on xmaster, rest is 00245 ! done automatically by call from libxbeach 00246 allocate(s%shipxCG (par%nship)) 00247 allocate(s%shipyCG (par%nship)) 00248 allocate(s%shipzCG (par%nship)) 00249 allocate(s%shipFx (par%nship)) 00250 allocate(s%shipFy (par%nship)) 00251 allocate(s%shipFz (par%nship)) 00252 allocate(s%shipMx (par%nship)) 00253 allocate(s%shipMy (par%nship)) 00254 allocate(s%shipMz (par%nship)) 00255 allocate(s%shipphi (par%nship)) 00256 allocate(s%shipchi (par%nship)) 00257 allocate(s%shippsi (par%nship)) 00258 endif 00259 endif 00260 00261 end subroutine ship_init 00262 00263 subroutine shipwave(s,par,sh) 00264 use params, only: parameters 00265 use xmpi_module 00266 use spaceparams 00267 use readkey_module 00268 use filefunctions 00269 use interp, only: grmap, grmap2, mkmap, linear_interp 00270 00271 implicit none 00272 00273 type(parameters) :: par 00274 type(spacepars),target :: s 00275 type(ship), dimension(:), pointer :: sh 00276 00277 integer :: i,ix,iy,shp_indx 00278 logical, save :: firstship=.true. 00279 real*8 :: shipx_old,shipy_old,dirship,radius,cosdir,sindir 00280 integer :: n1,n2,iprint=0 00281 real :: xymiss=-999 00282 00283 real*8, dimension(:,:),allocatable :: zsvirt 00284 !include 's.ind' 00285 !include 's.inp' 00286 00287 if (.not. allocated(zsvirt)) allocate(zsvirt(s%nx+1,s%ny+1)) 00288 zsvirt=s%zs+s%ph 00289 s%ph=0.d0 00290 00291 do i=1,par%nship 00292 00293 ! Find actual position and orientation of ship 00294 00295 shipx_old = s%shipxCG(i) 00296 shipy_old = s%shipyCG(i) 00297 call linear_interp(sh(i)%track_t,sh(i)%track_x,sh(i)%track_nt,par%t,s%shipxCG(i),shp_indx) 00298 call linear_interp(sh(i)%track_t,sh(i)%track_y,sh(i)%track_nt,par%t,s%shipyCG(i),shp_indx) 00299 if (sh(i)%flying==1) then 00300 call linear_interp(sh(i)%track_t,sh(i)%track_z,sh(i)%track_nt,par%t,s%shipzCG(i),shp_indx) 00301 else 00302 s%shipzCG(i) = 0.d0 00303 endif 00304 call linear_interp(sh(i)%track_t,sh(i)%track_dir,sh(i)%track_nt,par%t,dirship,shp_indx) 00305 radius=max(sh(i)%nx*sh(i)%dx,sh(i)%ny*sh(i)%dy)/2 00306 cosdir=cos(dirship) 00307 sindir=sin(dirship) 00308 00309 ! Compute pressures on ship based on water levels from XBeach 00310 !------------------------------------------------------------ 00311 ! Update locations of x ship grid points 00312 do iy=1,sh(i)%ny+1 00313 do ix=1,sh(i)%nx+1 00314 sh(i)%x(ix,iy)=s%shipxCG(i)+(ix-sh(i)%nx/2-1)*sh(i)%dx*cosdir - (iy-sh(i)%ny/2-1)*sh(i)%dy*sindir 00315 sh(i)%y(ix,iy)=s%shipyCG(i)+(ix-sh(i)%nx/2-1)*sh(i)%dx*sindir + (iy-sh(i)%ny/2-1)*sh(i)%dy*cosdir 00316 end do 00317 end do 00318 ! Interpolate XBeach water levels to ship grid when required 00319 n1=(s%nx+1)*(s%ny+1) 00320 n2=(sh(i)%nx+1)*(sh(i)%ny+1) 00321 ! Only carry out (costly) MKMAP once when ship is not moving 00322 !! If XBeach grid is regular rectangular a more efficient mapper can be used; TBD 00323 call MKMAP (s%wetz ,s%xz ,s%yz ,s%nx+1 ,s%ny+1 , & 00324 & sh(i)%x ,sh(i)%y ,n2 ,sh(i)%xs ,sh(i)%ys, & 00325 & sh(i)%nrx ,sh(i)%nry ,sh(i)%iflag ,sh(i)%nrin ,sh(i)%w , & 00326 & sh(i)%iref,iprint ,sh(i)%covered ,xymiss) 00327 call GRMAP (zsvirt ,n1 ,sh(i)%zs ,n2 ,sh(i)%iref, & 00328 & sh(i)%w ,4 ,iprint ) 00329 00330 if (sh(i)%compute_force==1) then 00331 ! Compute pressure head (m) on ship grid including small-scale motions 00332 !----------------------------------------------------- 00333 sh(i)%ph = sh(i)%zs-sh(i)%zhull 00334 00335 ! Compute forces on ship when required 00336 !------------------------------------- 00337 call ship_force(i,sh(i),s,par) 00338 00339 if (sh(i)%compute_motion==1) then 00340 ! Update vertical position and rotations when required 00341 !----------------------------------------------------- 00342 ! TBD 00343 00344 ! Compute actual z position of ship hull on ship grid when required 00345 !------------------------------------------------------------------ 00346 do iy=1,sh(i)%ny+1 00347 do ix=1,sh(i)%nx+1 00348 sh(i)%zhull(ix,iy)=s%shipzCG(i)-sh(i)%zCG-sh(i)%depth(ix,iy) & 00349 & -(sh(i)%x(ix,iy)-sh(i)%xCG)*sin(s%shipchi(i)) & 00350 & +(sh(i)%y(ix,iy)-sh(i)%yCG)*sin(s%shipphi(i)) 00351 enddo 00352 enddo 00353 endif ! compute_motion 00354 ! Compute pressure head (m) on ship grid ti feed back into XBeach 00355 !----------------------------------------------------- 00356 ! Next line to be implemented in combination with vertical motion only 00357 ! sh(i)%ph = sum(sh(i)%zs)/((sh(i)%nx+1)*(sh(i)%ny+1))-sh(i)%zhull 00358 sh(i)%ph = -sh(i)%zhull 00359 do iy=1,sh(i)%ny+1 00360 do ix=1,sh(i)%nx+1 00361 if (sh(i)%depth(ix,iy)==0) sh(i)%ph(ix,iy)=0.d0 00362 enddo 00363 enddo 00364 else 00365 sh(i)%ph = -sh(i)%zhull-s%shipzCG(i) 00366 sh(i)%ph = max(sh(i)%ph,0.d0) 00367 endif ! compute_forces 00368 00369 ! Compute pressure head (m) on XBeach grid 00370 !----------------------------------------- 00371 call grmap2(s%ph, s%dsdnzi ,n1 ,sh(i)%ph, sh(i)%dx*sh(i)%dy ,n2 ,sh(i)%iref, & 00372 & sh(i)%w ,4 ) 00373 ! do iy=1,s%ny+1 00374 ! do ix=1,s%nx+1 00375 ! ! Convert XBeach cell center coordinates to coordinates in local ship grid 00376 ! x1 = (xz(ix,iy)-shipxCG(i))*cosdir + (yz(ix,iy)-shipyCG(i))*sindir! 00377 ! y1 = -(xz(ix,iy)-shipxCG(i))*sindir + (yz(ix,iy)-shipyCG(i))*cosdir 00378 ! ! Convert to (float) grid number 00379 ! xrel=x1/sh(i)%dx+sh(i)%nx/2 00380 ! yrel=y1/sh(i)%dy+sh(i)%ny/2 00381 ! i1=floor(xrel) 00382 ! j1=floor(yrel) 00383 ! ! Carry out bilinear interpolation 00384 ! if (i1>=0 .and. i1<SH(i)%nx .and. j1>=0 .and. j1<sh(i)%ny) then 00385 ! s%ph(ix,iy)=s%ph(ix,iy)+(1.d0-(xrel-float(i1)))*(1.d0-(yrel-float(j1)))*sh(i)%depth(i1+1,j1+1) & 00386 ! +( xrel-float(i1) )*(1.d0-(yrel-float(j1)))*sh(i)%depth(i1+2,j1+1 ) & 00387 ! +(1.d0-(xrel-float(i1)))*( yrel-float(j1) )*sh(i)%depth(i1+1,j1+2) & 00388 ! +( xrel-float(i1) )*( yrel-float(j1) )*sh(i)%depth(i1+2,j1+2) 00389 ! s%ph(ix,iy)=s%ph(ix,iy)+(1.d0-(xrel-float(i1)))*(1.d0-(yrel-float(j1)))*sh(i)%ph(i1+1,j1+1) & 00390 ! +( xrel-float(i1) )*(1.d0-(yrel-float(j1)))*sh(i)%ph(i1+2,j1+1 ) & 00391 ! +(1.d0-(xrel-float(i1)))*( yrel-float(j1) )*sh(i)%ph(i1+1,j1+2) & 00392 ! +( xrel-float(i1) )*( yrel-float(j1) )*sh(i)%ph(i1+2,j1+2) 00393 00394 ! endif 00395 ! enddo 00396 ! enddo 00397 enddo 00398 00399 00400 ! Compute 00401 00402 if (firstship) then 00403 00404 ! apply initial condition 00405 00406 s%zs=s%zs-s%ph 00407 endif 00408 00409 firstship=.false. 00410 00411 end subroutine shipwave 00412 00413 subroutine ship_force(i,sh,s,par) 00414 use params, only: parameters 00415 use spaceparams 00416 type (ship) :: sh 00417 type (parameters) :: par 00418 type (spacepars),target :: s 00419 integer :: ix,iy,i 00420 real*8 :: dFx,dFy,dFz,hdx,hdy 00421 00422 !include 's.ind' 00423 !include 's.inp' 00424 00425 s%shipFx(i) = 0.d0 00426 s%shipFy(i) = 0.d0 00427 s%shipFz(i) = 0.d0 00428 s%shipMx(i) = 0.d0 00429 s%shipMy(i) = 0.d0 00430 s%shipMz(i) = 0.d0 00431 hdx=.5d0*sh%dx 00432 hdy=.5d0*sh%dy 00433 do iy=1,sh%ny 00434 do ix=1,sh%nx 00435 dFx=.5*(sh%ph(ix,iy)+sh%ph(ix+1,iy))*(sh%depth(ix+1,iy)-sh%depth(ix,iy))*sh%dy 00436 dFy=.5*(sh%ph(ix,iy)+sh%ph(ix,iy+1))*(sh%depth(ix,iy+1)-sh%depth(ix,iy))*sh%dx 00437 dFz=sh%ph(ix,iy)*sh%dx*sh%dy 00438 s%shipFx(i)=s%shipFx(i)+dFx 00439 s%shipFy(i)=s%shipFy(i)+dFy 00440 s%shipFz(i)=s%shipFz(i)+dFz 00441 s%shipMx(i)=s%shipMx(i)+((sh%y(ix,iy)+hdy)-sh%yCG)*dFz-(.5d0*(sh%zhull(ix,iy)+sh%zhull(ix,iy+1))-sh%zCG)*dFy 00442 s%shipMy(i)=s%shipMy(i)-((sh%x(ix,iy)+hdx)-sh%xCG)*dFz+(.5d0*(sh%zhull(ix,iy)+sh%zhull(ix+1,iy))-sh%zCG)*dFx 00443 !shipMz(i)=shipMz(i)-((sh%y(ix,iy)+hdy)-sh%yCG)*dFx+((sh%x(ix,iy)+hdx)-sh%xCG)*dFy 00444 s%shipMz(i)=s%shipMz(i)-(.5*(((iy-sh%ny/2-1)*sh%dy-sh%yCG)+((iy+1-sh%ny/2-1)*sh%dy-sh%yCG)))*dFx & 00445 & +(.5*(((ix-sh%nx/2-1)*sh%dx-sh%xCG)+((ix+1-sh%nx/2-1)*sh%dx-sh%xCG)))*dFy 00446 enddo 00447 enddo 00448 s%shipFx(i)=s%shipFx(i) *par%rho*par%g 00449 s%shipFy(i)=s%shipFy(i) *par%rho*par%g 00450 s%shipFz(i)=s%shipFz(i) *par%rho*par%g 00451 s%shipMx(i)=s%shipMx(i) *par%rho*par%g 00452 s%shipMy(i)=s%shipMy(i) *par%rho*par%g 00453 s%shipMz(i)=s%shipMz(i) *par%rho*par%g 00454 00455 00456 end subroutine ship_force 00457 00458 end module ship_module