XBeach
C:/repositories/XBeach/trunk/src/xbeachlibrary/ship.F90
Go to the documentation of this file.
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
 All Classes Files Functions Variables Defines