!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
module drifters,1
!BOP
! !MODULE: drifters
!
! !DESCRIPTION:
! This module contains routines necessary for moving drifters.
! NOTE: this module currently does not really exist - this version
! has old CM-5 code and will be replaced with some new code when
! Mat is ready to add it in...
!
! !REVISION HISTORY:
! SVN:$Id: drifters.F90 2290 2006-10-25 18:23:10Z njn01 $
!
! !USES
use kinds_mod
implicit none
private
save
!EOP
!BOC
!EOC
!***********************************************************************
contains
!***********************************************************************
!BOP
! !IROUTINE: init_drifters
! !INTERFACE:
subroutine init_drifters
! !DESCRIPTION:
! Initializes all variables for drifter diagnostics.
!
! !REVISION HISTORY:
! same as module
!EOP
!BOC
!-----------------------------------------------------------------------
!EOC
end subroutine init_drifters
!***********************************************************************
end module drifters
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
#ifdef does_not_exist
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
c begin file drifters.H
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
c-----------------------------------------------------------------------
c Some things for drifters. Note especially that there
c is a new 3D array here (WVEL) that can increase memory
c significantly if "drifter_particles != 2"
c-----------------------------------------------------------------------
integer, parameter :: ndrifters_max = 400
integer ndrifters, ndrifters_total, array_size, arrays_deployed
TYPE drifter_ijk(3,ndrifters_max)
common/drifter_stuff/ndrifters, ndrifters_total, array_size
& , arrays_deployed, drifter_ijk
#if drifter_particles != 2
TYPE, dimension(imt,jmt,km) :: WVEL
common/drifter_3d/WVEL
#endif
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
c end file drifters.H
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
c begin file drifter_read.F
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
subroutine drifter_read,6
c-----------------------------------------------------------------------
c Read in initial logical positions of drifters.
c
c If introducing drifters for the first time in a given run, then
c positions are calculated using the pre-processing routine drifter_init.f.
c If continuing a run then use the final dump on the out_drifters
c file from the previous portion of the run:
c tail -n out_drifters > in_drifters
c where n is the number of drifters
c then edit the in_drifters file and put n (# of drifters) on the
c first line
c addition of drifters resulted in alteration of the following files:
c scalar.H/d_scalar.H
c files.H
c in.dat
c in_files.dat
c makefile
c initial.F
c ocean.F
c baroclinic.F
c advu.F
c time_manager.F
c-----------------------------------------------------------------------
use grid
implicit none
integer np,i
real u,v,w,t,s
logical dump_grid_data
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c Use the following to dump grid data for use with pre-processor
c program drifter_init.f
c-----------------------------------------------------------------------
dump_grid_data = .false.
if(dump_grid_data) then
call get_unit(nu)
open(nu,file='out_grid.dat',form='unformatted',status='unknown')
write(nu)ulong
write(nu)ulat
write(nu)kmt
write(nu)angle
close(nu)
call release_unit(nu)
stop 'Wrote grid data to out_grid.dat'
endif
c-----------------------------------------------------------------------
c Read in drifter positions from the in_drifters file. These must
c be logical position (i,j,k), not physical (x,y,z).
c-----------------------------------------------------------------------
call get_unit(nu)
open(nu,file=in_drifters,status='old')
read(nu,*)ndrifters_total,array_size,arrays_deployed
if(ndrifters_total.gt.ndrifters_max) then
write(*,*)' ** WARNING: ndrifters > ndrifters_max'
write(*,*)' Only first ',ndrifters_max,
& ' drifters will be read in'
ndrifters_total = ndrifters_max
endif
do np = 1,ndrifters_total
read(nu,*)i,
& drifter_ijk(1,np),drifter_ijk(2,np),drifter_ijk(3,np)
enddo
write(*,*)' '
write(*,*)ndrifters_total,' drifters read in...'
write(*,*)arrays_deployed,' drifter arrays deployed...'
write(*,*)' Each drifter array has ',array_size,' elements'
write(*,*)' '
close(nu)
call release_unit(nu)
ndrifters = array_size*arrays_deployed
c-----------------------------------------------------------------------
c Open the drifters data file and write # of drifters.
c This file will stay open throughout the calculation.
c-----------------------------------------------------------------------
call get_unit(nu_drift)
open(nu_drift,file=out_drifters,status='unknown')
write(nu_drift,*)ndrifters_total
end subroutine drifter_read
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
c end file drifter_read.F
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
c begin file drifter_dump.F
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
subroutine drifter_dump(filename),3
c-----------------------------------------------------------------------
c Write logical positions of drifters.
c-----------------------------------------------------------------------
use io
implicit none
integer np
integer (kind=int_kind) ::
& nu ! i/o unit attached to output file
character*80 filename
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c Write in drifter positions in logical space for restart.
c-----------------------------------------------------------------------
call get_unit(nu)
open(nu,file=filename,status='unknown')
write(nu,*)ndrifters_total,array_size,arrays_deployed
do np = 1,ndrifters_total
write(nu,*)np,
& drifter_ijk(1,np),drifter_ijk(2,np),drifter_ijk(3,np)
enddo
close(nu)
call release_unit(nu)
write(stdout,*)' '
write(stdout,*)'file written: ',filename
end subroutine drifter_dump
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
c end file drifter_dump.F
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
c begin file drifter_intp.F
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
real function drifter_intp (A1,A2, i, ip, j, jp, wijk, 8
& wijpk, wijkp, wijpkp, wipjk, wipjpk, wipjkp, wipjpkp)
c-----------------------------------------------------------------------
c This function tri-linearly interpolates a field to a point
c given an array at 2 k-levels (A1,A2) and the weights (w*).
c-----------------------------------------------------------------------
implicit none
integer i, ip, j, jp
real wijk, wipjk, wijpk, wipjpk, wijkp,
& wipjkp, wijpkp, wipjpkp
real, dimension(imt,jmt) :: A1, A2 ! CM array
c***********************************************************************
drifter_intp = wijk*A1(i,j) +
& wijpk*A1(i,jp) +
& wijkp*A2(i,j) +
& wijpkp*A2(i,jp) +
& wipjk*A1(ip,j) +
& wipjpk*A1(ip,jp) +
& wipjkp*A2(ip,j) +
& wipjpkp*A2(ip,jp)
end
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
c end file drifter_intp.F
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
c begin file drifter_move.F
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
subroutine drifter_move 1,6
c-----------------------------------------------------------------------
c Move each drifter using a predictor-corrector scheme. Trilinear
c interpolation is used to find velocity at drifter locations.
c Drifters are moved in logical space=> array drifter_ijk holds
c the logical positions of the drifters. Only when drifter properties
c are dumped (in drifter_prop) is the conversion from logical space
c to physical (lat-long-depth) performed. Algorithm provided
c by Doug Kothe.
c-----------------------------------------------------------------------
use grid
use time_management
use prognostic
implicit none
integer i, ip, j, jp, k, kp, np, im, jm, kk
real drifter_intp
real alpho, beto, dalpho_dx, dbeto_dy, dgammo_dz, dtpred,
& gammo, gammw, up, vp, wp, xp, yp, zp
real wijk, wipjk, wijpk, wipjpk, wijkp, wipjkp, wijpkp,
& wipjpkp
logical predictor
real, dimension(imt,jmt) :: U1, U2, TEMP
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c Cycle through all drifters
c-----------------------------------------------------------------------
do 1000 np = 1,ndrifters
c-----------------------------------------------------------------------
c Predictor delta-t
c-----------------------------------------------------------------------
predictor = .true.
dtpred = p5*dtu
if (mod(itt,navg).eq.0) then ! correct for averaging timesteps
dtpred = p5*dtpred
endif
c-----------------------------------------------------------------------
c Initialize the particle positions
c-----------------------------------------------------------------------
xp = drifter_ijk(1,np)
yp = drifter_ijk(2,np)
zp = drifter_ijk(3,np)
c-----------------------------------------------------------------------
c Get the cell indices containing the particle and surrounding cells
c-----------------------------------------------------------------------
10 i = int(xp)
j = int(yp)
k = int(zp)
if(i.gt.imt) i = 1
if(i.lt.1 ) i = imt
j = min(max(j,1),jmt)
k = min(max(k,0),km) ! k can be zero
ip = i + 1
jp = j + 1
kp = k + 1
im = i - 1
jm = j - 1
if(ip.gt.imt) ip = 1
if(im.lt.1 ) im = imt
jp = min(max(jp,1),jmt)
jm = min(max(jm,1),jmt)
kp = min(max(kp,1),km)
c-----------------------------------------------------------------------
c Compute weights for the tri-linear interpolation coefficients
c-----------------------------------------------------------------------
alpho = xp - i
beto = yp - j
gammw = zp - k
kk = int(zp + p5)
kp = kk + 1
kk = max(kk,1)
kp = min(kp,km)
if(gammw.lt.p5) then
gammo = (p5*dz(kk) + gammw*dz(kp))*dzwr(kk)
else
gammo = (gammw - p5)*dz(kk)*dzwr(kk)
endif
c-----------------------------------------------------------------------
c Compute the tri-linear interpolation coefficients
c-----------------------------------------------------------------------
wijk = (1. - alpho)*(1. - beto)*(1. - gammo)
wijpk = (1. - alpho)*beto*(1. - gammo)
wijkp = (1. - alpho)*(1. - beto)*gammo
wijpkp = (1. - alpho)*beto*gammo
wipjk = alpho*(1. - beto)*(1. - gammo)
wipjpk = alpho*beto*(1. - gammo)
wipjkp = alpho*(1. - beto)*gammo
wipjpkp = alpho*beto*gammo
c-----------------------------------------------------------------------
c Velocity components are already in local frame
c First do U
c-----------------------------------------------------------------------
U1 = UVEL(:,:,kk,curtime,1)
U2 = UVEL(:,:,kp,curtime,1)
c testing
c U1 = 0.
c U2 = 0.
c U1 = 100.*cos(ANGLE) ! zonal only
c U2 = 100.*cos(ANGLE) ! zonal only
c U1 = 100.*sin(ANGLE) ! meridional only
c U2 = 100.*sin(ANGLE) ! meridional only
c-----------------------------------------------------------------------
c Interpolate U velocity to the particle
c-----------------------------------------------------------------------
up = drifter_intp(U1, U2, i, ip, j, jp, wijk, wijpk,
& wijkp, wijpkp, wipjk, wipjpk, wipjkp, wipjpkp)
c-----------------------------------------------------------------------
c Now do V
c-----------------------------------------------------------------------
U1 = VVEL(:,:,kk,curtime,1)
U2 = VVEL(:,:,kp,curtime,1)
c testing
c U1 = 0.
c U2 = 0.
c U1 = 100.*sin(-ANGLE) ! zonal only
c U2 = 100.*sin(-ANGLE) ! zonal only
c U1 = 100.*cos(ANGLE) ! meridional only
c U2 = 100.*cos(ANGLE) ! meridional only
c-----------------------------------------------------------------------
c Interpolate V velocity to the particle
c-----------------------------------------------------------------------
vp = drifter_intp(U1, U2, i, ip, j, jp, wijk, wijpk,
& wijkp, wijpkp, wipjk, wipjpk, wipjkp, wipjpkp)
#if drifter_particles != 2
c-----------------------------------------------------------------------
c Now do W differently because it is located at cell bottom
c Note that the horizontal weights are the same as for U,V
c since WVEL is calculated in advu (still cell-bottom, though).
c Also, the vertical weights are different--they are simply
c gammw since the vertical index for the drifters is based on
c the W-grid.
c-----------------------------------------------------------------------
if(k.eq.0) then
U1 = DHU
else
U1 = WVEL(:,:,k)
endif
if(k.ge.km) then
U2 = 0. ! should not happen
else
U2 = WVEL(:,:,k + 1)
endif
c testing...
c U1 = 0.01
c U2 = 0.01
c-----------------------------------------------------------------------
c Compute the tri-linear interpolation coefficients
c-----------------------------------------------------------------------
wijk = (1. - alpho)*(1. - beto)*(1. - gammw)
wijpk = (1. - alpho)*beto*(1. - gammw)
wijkp = (1. - alpho)*(1. - beto)*gammw
wijpkp = (1. - alpho)*beto*gammw
wipjk = alpho*(1. - beto)*(1. - gammw)
wipjpk = alpho*beto*(1. - gammw)
wipjkp = alpho*(1. - beto)*gammw
wipjpkp = alpho*beto*gammw
c-----------------------------------------------------------------------
c Interpolate W velocity to the particle
c-----------------------------------------------------------------------
wp = drifter_intp(U1, U2, i, ip, j, jp, wijk, wijpk,
& wijkp,wijpkp, wipjk, wipjpk, wipjkp, wipjpkp)
#else
c-----------------------------------------------------------------------
c Drifters are drogued => set vertical velocity to zero
c-----------------------------------------------------------------------
wp = c0
#endif
c-----------------------------------------------------------------------
c Compute the derivative of the metric
c Remember that z is positive up and k increases down, so
c dgammo_dz is actually negative.
c-----------------------------------------------------------------------
dalpho_dx = (1. - alpho)/HUS(i,jp) + alpho/HUS(ip,jp)
dbeto_dy = (1. - beto)/HUW(ip,j) + beto/HUW(ip,jp)
c dgammo_dz = (1. - gammo)*dzr(kk) + gammo*dzr(kp) ! original
dgammo_dz = (1. - gammw)*dzwr(k) + gammw*dzwr(k+1)
dgammo_dz = - dgammo_dz
c-----------------------------------------------------------------------
c Move the drifters and correct for longitudinal periodicity and
c drifters moving through boundaries (ungracefully).
c-----------------------------------------------------------------------
xp = drifter_ijk(1,np) + dtpred*up*dalpho_dx
yp = drifter_ijk(2,np) + dtpred*vp*dbeto_dy
zp = drifter_ijk(3,np) + dtpred*wp*dgammo_dz
if(xp.gt.imt) xp = xp - imt
if(xp.lt.c0 ) xp = xp + imt
if(yp.gt.jmt) yp = jmt
if(yp.lt.c0 ) yp = 1.
if(zp.gt.(km + p5) ) zp = km + p5
if(zp.lt.p5 ) zp = p5
600 continue
c-----------------------------------------------------------------------
c Check to see if this is the predictor or corrector step
c-----------------------------------------------------------------------
if (predictor) then ! do corrector step
dtpred = 2.*dtpred
predictor = .false.
goto 10
endif
c-----------------------------------------------------------------------
c Done with this drifter=> save its position
c-----------------------------------------------------------------------
drifter_ijk(1,np) = xp
drifter_ijk(2,np) = yp
drifter_ijk(3,np) = zp
1000 continue ! go do next drifter
end subroutine drifter_move
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
c end file drifter_move.F
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
c begin file drifter_prop.F
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
subroutine drifter_prop 1,8
c-----------------------------------------------------------------------
c Calculate field quantities at location of all drifters using
c trilinear interpolation. Information is then write to the file
c specified by out_drifters in the in_files.dat file. this file is
c opened in drifter_read and remains open throughout the calculation.
c-----------------------------------------------------------------------
use grid
use time_management
use prognostic
implicit none
real drifter_intp
integer np,i,j,ip,jp,k,kk,kp
real dx,dy,zz,dlong,dlat,tan_alpha,dulong,zd
real drifter_long,drifter_lat,drifter_z,up,vp,wp,tempp,saltp
real wijk, wipjk, wijpk, wipjpk, wijkp, wipjkp, wijpkp,
& wipjpkp
real, dimension(imt,jmt) :: U1, U2, TEMP
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c Write the time to the out_drifters file
c-----------------------------------------------------------------------
write(37,*)tday
c-----------------------------------------------------------------------
c Cycle over all drifters. Include even those that are not
c deployed yet.
c-----------------------------------------------------------------------
do np = 1,ndrifters_total
c-----------------------------------------------------------------------
c Find each particle position in physical coordinates
c-----------------------------------------------------------------------
i = int(drifter_ijk(1,np))
j = int(drifter_ijk(2,np))
k = int(drifter_ijk(3,np))
dx = drifter_ijk(1,np) - i
dy = drifter_ijk(2,np) - j
zz = drifter_ijk(3,np) - k ! cannot use dz
ip = i + 1
if(ip.gt.imt) ip = 1
jp = j + 1
if(jp.gt.jmt) jp = j
kk = int(drifter_ijk(3,np) + p5)
kp = kk + 1
kk = max(kk,1)
kp = min(kp,km)
if(zz.lt.p5) then
zd = (p5*dz(kk) + zz*dz(kp))*dzwr(kk)
else
zd = (zz - p5)*dz(kk)*dzwr(kk)
endif
c-----------------------------------------------------------------------
c Check for crossing 0/360 longitude
c-----------------------------------------------------------------------
dulong = ULONG(ip,j) - ULONG(i,j)
if(dulong.lt.-5.0) then
dulong = ULONG(ip,j) - ULONG(i,j) + 2.*pi
elseif(dulong.gt.5.0) then
dulong = ULONG(ip,j) - ULONG(i,j) - 2.*pi
endif
c-----------------------------------------------------------------------
c First find latitude/longitude
c-----------------------------------------------------------------------
if(dx.ne.c0) then
tan_alpha = dy*HTE(i,jp)/(dx*HTN(ip,j))
dlong = dx*dulong*(c1 - tan(ANGLE(i,j))*tan_alpha)
if(dy.ne.c0) then
dlat = dy*(ULAT(i,jp) - ULAT(i,j))
& *(c1 + tan(ANGLE(i,j))/tan_alpha)
else
dlat = c0
endif
else
dlong = c0
dlat = dy*(ULAT(i,jp) - ULAT(i,j))
endif
drifter_long = ULONG(i,j) + dlong
drifter_lat = ULAT (i,j) + dlat
c-----------------------------------------------------------------------
c Now find vertical position
c-----------------------------------------------------------------------
if(k.eq.0) then
drifter_z = zz*dz(k+1)
else
drifter_z = zz*dz(k+1) + zw(k)
endif
c-----------------------------------------------------------------------
c Convert to degrees, meters
c-----------------------------------------------------------------------
drifter_long = 180.*drifter_long/pi
if(drifter_long.gt.360.0) drifter_long = drifter_long - 360.
if(drifter_long.lt.0.0) drifter_long = drifter_long + 360.
drifter_lat = 180.*drifter_lat/pi
drifter_z = .01*drifter_z
c-----------------------------------------------------------------------
c Calculate interpolation coefficients for calculating
c horizontal velocity at drifter location.
c
c Note that U and V are defined relative to the local coordinate
c frame so we need to rotate them to get true zonal
c and meridional velocities.
c-----------------------------------------------------------------------
wijk = (1. - dx)*(1. - dy)*(1. - zd)
wijpk = (1. - dx)*dy*(1. - zd)
wijkp = (1. - dx)*(1. - dy)*zd
wijpkp = (1. - dx)*dy*zd
wipjk = dx*(1. - dy)*(1. - zd)
wipjpk = dx*dy*(1. - zd)
wipjkp = dx*(1. - dy)*zd
wipjpkp = dx*dy*zd
c-----------------------------------------------------------------------
c Interpolate U and V velocity to the particle
c-----------------------------------------------------------------------
U1 = UVEL(:,:,kk,curtime,1)*cos(ANGLE) +
& VVEL(:,:,kk,curtime,1)*sin(-ANGLE)
U2 = UVEL(:,:,kp,curtime,1)*cos(ANGLE) +
& VVEL(:,:,kp,curtime,1)*sin(-ANGLE)
up = drifter_intp(U1, U2, i, ip, j, jp, wijk, wijpk,
& wijkp, wijpkp, wipjk, wipjpk, wipjkp, wipjpkp)
U1 = VVEL(:,:,kk,curtime,1)*cos(ANGLE) -
& UVEL(:,:,kk,curtime,1)*sin(-ANGLE)
U2 = VVEL(:,:,kp,curtime,1)*cos(ANGLE) -
& UVEL(:,:,kp,curtime,1)*sin(-ANGLE)
vp = drifter_intp(U1, U2, i, ip, j, jp, wijk, wijpk,
& wijkp, wijpkp, wipjk, wipjpk, wipjkp, wipjpkp)
#if drifter_particles != 2
c-----------------------------------------------------------------------
c Now do W differently because it is located at cell bottom
c Note that the horizontal weights are the same as for U,V
c since WVEL is calculated in advu (still cell-bottom, though).
c Also, the vertical weights are different--they are simply
c zz since the vertical index for the drifters is based on
c the W-grid.
c-----------------------------------------------------------------------
if(k.eq.0) then
U1 = DHU
else
U1 = WVEL(:,:,k)
endif
if(k.ge.km) then
U2 = 0. ! should not happen
else
U2 = WVEL(:,:,k + 1)
endif
c-----------------------------------------------------------------------
c Compute the tri-linear interpolation coefficients
c-----------------------------------------------------------------------
wijk = (1. - dx)*(1. - dy)*(1. - zz)
wijpk = (1. - dx)*dy*(1. - zz)
wijkp = (1. - dx)*(1. - dy)*zz
wijpkp = (1. - dx)*dy*zz
wipjk = dx*(1. - dy)*(1. - zz)
wipjpk = dx*dy*(1. - zz)
wipjkp = dx*(1. - dy)*zz
wipjpkp = dx*dy*zz
c-----------------------------------------------------------------------
c Interpolate W velocity to the particle
c-----------------------------------------------------------------------
wp = drifter_intp(U1, U2, i, ip, j, jp, wijk, wijpk,
& wijkp,wijpkp, wipjk, wipjpk, wipjkp, wipjpkp)
#else
c-----------------------------------------------------------------------
c Drifters are drogued so set W = 0
c-----------------------------------------------------------------------
wp = c0
#endif
c-----------------------------------------------------------------------
c Now do Temperature and Salinity. The horizontal weights need
c to be modified due to the staggered grid. Vertical remains
c the same as with U,V.
c-----------------------------------------------------------------------
i = int(drifter_ijk(1,np) + 0.5)
j = int(drifter_ijk(2,np) + 0.5)
if(j.gt.jmt) j = jmt
dx = drifter_ijk(1,np) - i + 0.5
dy = drifter_ijk(2,np) - j + 0.5
if(i.gt.imt) then
i = 1
dx = drifter_ijk(1,np) - imt + 0.5
endif
ip = i + 1
if(ip.gt.imt) ip = 1
jp = j + 1
if(jp.gt.jmt) jp = j
wijk = (1. - dx)*(1. - dy)*(1. - zd)
wijpk = (1. - dx)*dy*(1. - zd)
wijkp = (1. - dx)*(1. - dy)*zd
wijpkp = (1. - dx)*dy*zd
wipjk = dx*(1. - dy)*(1. - zd)
wipjpk = dx*dy*(1. - zd)
wipjkp = dx*(1. - dy)*zd
wipjpkp = dx*dy*zd
c-----------------------------------------------------------------------
c Interpolate Temperature and Salinity to the particle
c-----------------------------------------------------------------------
U1 = TRACER(:,:,kk,1,curtime,1)
U2 = TRACER(:,:,kp,1,curtime,1)
tempp = drifter_intp(U1, U2, i, ip, j, jp, wijk, wijpk,
& wijkp, wijpkp, wipjk, wipjpk, wipjkp, wipjpkp)
U1 = TRACER(:,:,kk,2,curtime,1)
U2 = TRACER(:,:,kp,2,curtime,1)
saltp = drifter_intp(U1, U2, i, ip, j, jp, wijk, wijpk,
& wijkp, wijpkp, wipjk, wipjpk, wipjkp, wipjpkp)
saltp = saltp*1000. ! convert to sigma units
c-----------------------------------------------------------------------
c Write data to out_drifters file
c-----------------------------------------------------------------------
write(37,'(i5,8e13.5)')np,drifter_long,drifter_lat,drifter_z,
& up,vp,wp,tempp,saltp
enddo
end subroutine drifter_prop
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
c end file drifter_prop.F
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
#endif