#if defined(ROW_LAND)
#define SEA_P .true.
#define SEA_U .true.
#define SEA_V .true.
#elif defined(ROW_ALLSEA)
#define SEA_P allip(j).or.ip(i,j).ne.0
#define SEA_U alliu(j).or.iu(i,j).ne.0
#define SEA_V alliv(j).or.iv(i,j).ne.0
#else
#define SEA_P ip(i,j).ne.0
#define SEA_U iu(i,j).ne.0
#define SEA_V iv(i,j).ne.0
#endif
module mod_floats
use mod_xc ! HYCOM communication interface
use mod_pipe ! HYCOM debugging interface
#if defined(STOKES)
use mod_stokes ! Stokes Drift Velocity Module
#endif
!
implicit none
!
! --- HYCOM synthetic floats, drifters and moorings
! --- See subroutine floats (below) for more information
!
integer, parameter, public :: nfldim=400 !maximum number of synthetic floats
real, allocatable, dimension(:,:,:), save, public :: &
wveli, & ! interface vertical velocity
uold2, &
vold2, &
wold2u, &
wold2d
real, allocatable, dimension(:,:), save, public :: &
dlondx, &
dlondy, &
dlatdx, &
dlatdy
real, allocatable, dimension(:,:,:), save, private :: &
wvelup, &
wveldn
real, allocatable, dimension(:,:), save, private :: &
dpdxup, &
dpdyup, &
dpdxdn, &
dpdydn
real, save :: flt(nfldim,13), &
deltfl,fldepm,tbvar,tdecri,dtturb,uturb0
integer, save :: kfloat(nfldim),iflnum(nfldim),ifltyp(nfldim), &
nflsam,nfldta,fltflg,nfladv,intpfl,iturbv,ismpfl, &
nflout,nfl,nflt,nstepfl
logical, save :: synflt,turbvel,samplfl,wvelfl,hadvfl,nonlatlon
character(48), save :: flnmflti,flnmflto,flnmfltio
contains
subroutine floats_init(m,n,time0)
use mod_cb_arrays ! HYCOM saved arrays
implicit none
integer m,n
real time0
!
! --- read initial float data
! --- initialize parameters and arrays required for floats
!
integer i,j,k,l,margin
integer nbcday,nsmday,ityp
!
# include "stmt_fns.h"
!
! --- initialize flags
!
! --- wvelfl: calculate vertical velocity?
wvelfl=.false.
!
! --- hadvfl: horizontally advect the float?
hadvfl=.false.
!
! --- nonlatlon: does the model grid cross latitude/longitude lines anywhere
! --- in the domain?
nonlatlon=.false.
!
! --- initialize file names
flnmflti = 'floats.input'
flnmflto = 'floats_out'
flnmfltio= 'floats.input_out'
!
! -------------------------
! --- set timing parameters
! -------------------------
!
! --- must have an integer number of baroclinic time steps per day
nbcday=nint(86400.0/baclin)
if (real(nbcday).ne.86400.0/baclin) then
if (mnproc.eq.1) then
write(lp,'(/ a /)') &
'error - need integer no. of bacl. time steps/day'
call flush(lp)
endif !1st tile
call xcstop('(floats_init)')
stop '(floats_init)'
endif
!
! --- parameter nfladv is entered in blkdat.input as the number of
! --- baroclinic time steps between updates of float position. at model
! --- time nstep when the float is to be advected, velocity fields saved
! --- at times nstep, nstep-nfladv/2, and nstep-nfladv are used to perform
! --- the runga-kutta time interpolation.
!
! --- nfladv must be set so the float is advected every 2-4 hours, but
! --- also must be no smaller than 4. This ensures that the runga-
! --- kutta time interpolation is performed with a delta-time of 1-2 hours.
!
! --- if all floats are stationary (synthetic moorings), then nfladv
! --- must still be no smaller than 4, but the above time restrection
! --- is not necessary. this allows the user to sample at very high
! --- frequency. for example, if baclin is 360 seconds, then setting
! --- nfladv to 10 will allow water properties to be sampled once per hour
!
! --- set variable nfldta to nfladv/2 so that it equals the number of time
! --- steps separating the velocity fields input into the runga-kutta
! --- interpolation
!
nfldta=nfladv/2
!
! --- make sure that the float will be advected at an integer number of
! --- temporal points per day
nsmday=nbcday/nfladv
if (int(real(nbcday)/real(2*nfldta)).ne.nsmday) then
if (mnproc.eq.1) then
write(lp,'(/ a /)') &
'error - need integer no. of adv. times per day'
call flush(lp)
endif !1st tile
call xcstop('(floats_init)')
stop '(floats_init)'
endif
!
! --- test to make sure that the advection time interval is between 2
! --- and 4 hours. only stop the program later if there are floats to
! --- be advected
if (nfladv.ne.4 .and. (nsmday.gt.12 .or. nsmday.lt.6)) then
nfladv=-nfladv
endif
!
! --- initialize number of time steps between float output
if (nflsam.gt.0) then
nflsam=nflsam*nbcday
elseif (nflsam.eq.0) then
nflsam=abs(nfladv)
endif
!
! --- calculate the time interval for the runga-kutta interpolation
! --- as 1/2 of the advection time interval
deltfl=nfldta*baclin
!
! --- set minimum float depth (m)
fldepm=1.0
!
! --- -------------------------------------------------------------
! --- initialize synthetic drifters, floats, and moored instruments
! --- -------------------------------------------------------------
!
open(unit=uoff+99,file=trim(flnminp)//flnmflti,status='old') !on all nodes
!
! --- first line in float input file is the number of floats
read(uoff+99,*) nflt
if (nflt.gt.nfldim) then
if (mnproc.eq.1) then
write(lp,'(/ a /)') &
'error - increase nfldim in mod_floats.F'
call flush(lp)
endif !1st tile
call xcstop('(floats_init)')
stop '(floats_init)'
endif
if (nflt.eq.0) then
synflt=.false.
endif
!
! --- second line in float input file is the initial float time step
! --- this number is set to zero for the first model run segment
read(uoff+99,*) nstepfl
!
do nfl=1,nflt
!
! --- read input file containing the following information for each float:
! ---
! --- column 1 - initial sequential float number
! --- column 2 - float type
! --- 1 = 3-d lagrangian (vertically advected by diagnosed w)
! --- 2 = isopycnic
! --- 3 = isobaric (surface drifter when released in sfc. layer)
! --- 4 = stationary (synthetic moored instrument)
! --- column 3 - deployment time (days from model start, 0.0 = immediate)
! --- column 4 - termination time (days from model start, 0.0 = forever)
! --- column 5 - initial longitude (must be between minimum and maximum
! --- longitudes defined in regional.grid.b)
! --- column 6 - initial latitude (must be between minimum and maximum
! --- latitudes defined in regional.grid.b)
! --- column 7 - initial depth (or reference sigma for isopycnic floats)
!
read(uoff+99,*) iflnum(nfl),ifltyp(nfl),flt(nfl,8),flt(nfl,9), &
flt(nfl,1),flt(nfl,2),flt(nfl,3)
if (ifltyp(nfl).eq.1 .or. ifltyp(nfl).eq.4) then
wvelfl=.true.
endif
if (ifltyp(nfl).ne.4) then
hadvfl=.true.
!
! --- since this float is to be advected, stop if the advection time interval
! --- is not between 2 and 4 hr
if (nfladv.lt.0) then
if (mnproc.eq.1) then
write(lp,'(/ a /)') &
'error - set nfladv to advect floats every 2-4 hr'
call flush(lp)
endif !1st tile
call xcstop('(floats_init)')
stop '(floats_init)'
endif
!
endif
!
flt(nfl,4)=0.0
flt(nfl,5)=0.0
flt(nfl,6)=0.0
flt(nfl,7)=0.0
if(ifltyp(nfl).eq.2) then
flt(nfl,7)=flt(nfl,3)
flt(nfl,3)=0.0
endif
flt(nfl,10)=0.0
flt(nfl,11)=0.0
flt(nfl,12)=0.0
flt(nfl,13)=0.0
!
kfloat(nfl)=-1
!
enddo
!
nfladv=abs(nfladv)
!
close(unit=uoff+99) !file='float.input'
!
if (mnproc.eq.1) then
write(lp,955) nflt
write(lp,956) nfladv
write(lp,957) nflsam
write(lp,958) min(nflt,10)
955 format(' read initial data for',i6,' floats, time step',i9)
956 format(' float advected every',i7,' baroclinic time steps')
957 format(' float sampled every',i8,' baroclinic time steps'/)
958 format(' input data for first',i3,' floats')
do nfl=1,min(nflt,10)
write(lp,959) nfl,(flt(nfl,i),i=1,9)
959 format(3x,i2,1p,5e13.5/5x,4e13.5/)
enddo
call flush(lp)
endif !1st tile
!
! ---------------------------------------------------------------------------
! --- allocate arrays for floats
! ---------------------------------------------------------------------------
!
allocate( wvelup(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,kdm), &
wveldn(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,kdm), &
dpdxup(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
dpdyup(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
dpdxdn(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
dpdydn(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) )
allocate( wveli( 1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy, kdm+1), &
uold2( 1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,2*kdm), &
vold2( 1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,2*kdm), &
wold2u(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,2*kdm), &
wold2d(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,2*kdm), &
dlondx(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
dlondy(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
dlatdx(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
dlatdy(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) )
call mem_stat_add( 11*(idm+2*nbdy)*(jdm+2*nbdy)*kdm )
call mem_stat_add( 9*(idm+2*nbdy)*(jdm+2*nbdy) )
!
if (synflt) then
!$OMP PARALLEL DO PRIVATE(j,i,k) &
!$OMP SCHEDULE(STATIC,jblk)
do j=1-nbdy,jj+nbdy
do i=1-nbdy,ii+nbdy
do k=1,kk
uold2( i,j,k )=hugel
uold2( i,j,k+kk)=hugel
vold2( i,j,k )=hugel
vold2( i,j,k+kk)=hugel
wold2u(i,j,k )=hugel
wold2u(i,j,k+kk)=hugel
wold2d(i,j,k )=hugel
wold2d(i,j,k+kk)=hugel
wveli( i,j,k) =hugel
enddo !k
wveli(i,j,kk+1)=hugel
enddo
enddo
!$OMP END PARALLEL DO
endif !synflt
!
! ---------------------------------------------------------------------------
! --- initialize old horizontal velocities for runga-kutta time interpolation
! --- calculate dlondx, dlondy, dlatdx, dlatdy required to convert u, v to
! --- longitude/time and latitude/time for float advection
! ---------------------------------------------------------------------------
!
margin = nbdy - 1
!
do j=1-margin,jj+margin
do i=1-margin,ii+margin
if (SEA_U) then
do k=1,kk
uold2(i,j,k )=u(i,j,k,n)+ubavg(i,j,n)
uold2(i,j,k+kk)=u(i,j,k,n)+ubavg(i,j,n)
#if defined(STOKES)
uold2(i,j,k )=uold2(i,j,k )+usd(i,j,k)
uold2(i,j,k+kk)=uold2(i,j,k+kk)+usd(i,j,k)
#endif
enddo !k
dlondx(i,j)=(plon(i ,j)*scpy(i ,j)- &
plon(i-1,j)*scpy(i-1,j))/scu2(i,j)
dlatdx(i,j)=(plat(i ,j)*scpy(i ,j)- &
plat(i-1,j)*scpy(i-1,j))/scu2(i,j)
if (plat(i,j)-plat(i-1,j).ne.0.0) then
nonlatlon=.true.
endif
endif !iu
if (SEA_V) then
do k=1,kk
vold2(i,j,k )=v(i,j,k,n)+vbavg(i,j,n)
vold2(i,j,k+kk)=v(i,j,k,n)+vbavg(i,j,n)
#if defined(STOKES)
vold2(i,j,k )=vold2(i,j,k )+vsd(i,j,k)
vold2(i,j,k+kk)=vold2(i,j,k+kk)+vsd(i,j,k)
#endif
enddo !k
dlondy(i,j)=(plon(i,j )*scpx(i,j )- &
plon(i,j-1)*scpx(i,j-1))/scv2(i,j)
dlatdy(i,j)=(plat(i,j )*scpx(i,j )- &
plat(i,j-1)*scpx(i,j-1))/scv2(i,j)
if (plon(i,j)-plon(i,j-1).ne.0.0) then
nonlatlon=.true.
endif
endif !iv
if (SEA_P) then
do k=1,kk
wold2u(i,j,k )=0.0
wold2d(i,j,k )=0.0
wold2u(i,j,k+kk)=0.0
wold2d(i,j,k+kk)=0.0
wveli (i,j,k )=0.0
enddo !k
wveli(i,j,kk+1)=0.0
endif !ip
enddo !i
enddo !j
!
! ----------------------------------------
! --- turbulent horizontal velocity option
! ----------------------------------------
!
! --- initialize turbulence constants
!
! --- tdecri is in inverse days and dtturb is the turbulent time step in days
!
! --- dtturb is set to the runga-kutta interpolation time step, which must
! --- be within 1 and 2 hours
!
if (turbvel) then
dtturb=deltfl/86400.0
uturb0=sqrt(2.*tbvar*tdecri*dtturb)
endif
!
return
end subroutine floats_init
subroutine floats_restart
implicit none
!
! -----------------------------
! --- output float restart file
! -----------------------------
!
integer i,k,l
!
if (mnproc.eq.1) then
!
! --- first remove floats that have run aground or left the domain
i=0
do l=1,nflt
if (flt(l,1).gt.-999.0) then
i=i+1
endif
enddo !1
!
write(lp,*) 'writing new float input file for restart'
call flush(lp)
!
open(unit=uoff+802,file=flnmfltio, &
form='formatted',status='new')
write(uoff+802,970) i
write(uoff+802,970) nstepfl
do l=1,nflt
if (flt(l,1).gt.-999.0) then
if (ifltyp(l).ne.2) then
write(uoff+802,971) iflnum(l),ifltyp(l),flt(l,8),flt(l,9), &
(flt(l,k),k=1,3)
else
write(uoff+802,971) iflnum(l),ifltyp(l),flt(l,8),flt(l,9), &
(flt(l,k),k=1,2),flt(l,7)
endif
endif
enddo !l
close(uoff+802)
endif !1st tile
970 format(i9)
971 format(i6,i2,2f9.2,3f15.9)
return
end subroutine floats_restart
subroutine floats(m,n,timefl,ioflag)
use mod_cb_arrays ! HYCOM saved arrays
implicit none
!
integer, parameter :: nfl_debug = -1 !no debugging
! integer, parameter :: nfl_debug = 6 !debug float nfl_debug
!
integer m,n,ioflag
real timefl
!
! --- local variables
integer i,j,k,l,margin
real tflt,sflt,thlo,thhi,thflt, &
uflt1,uflt2,uflt3,uflt4,uflt, &
vflt1,vflt2,vflt3,vflt4,vflt, &
wflt1,wflt2,wflt3,wflt4,wflt, &
vkflt,tkflt,skflt,xpos0,xpos1,xpos2,xpos3, &
ypos0,ypos1,ypos2,ypos3,depflt,plo,phi,q, &
depiso,wvhi,wvlo,uvfctr,qisop,ufltm,vfltm, &
alomin,alomax,alamin,alamax
real uturb,vturb,uturb1,vturb1,dlodx,dlody,dladx,dlady
real*4 rtab(200,2)
integer ifltll(3),jfltll(3),ngood(3)
integer kflt,k1,k2,k3,ier,ngrid,kold1w,kold2w,kold1,kold2, &
ntermn,i1,j1,ngoodi
integer ied,iede,ifladv
!
integer*4 seed,numran,inisee,iflag
integer*4 iseed1,iseed2
equivalence ( iseed1, iseed2, seed )
!
! --- local velocity component storage for synthetic moorings
real fltloc(nflt,3)
!
! --- 2-d local arrays used for the 16-point box horizontal interpolation
real varb2d(4,4),ptlon(4,4,3),ptlat(4,4,3)
logical maskpt(4,4,3),maskpi(4,4)
!
! --- arrays required for parallelization
real flt1(12*nflt),flt3(17*nflt)
integer iproc
!
real timer1, timer2, totime
real*8 wtime
!
# include "stmt_fns.h"
!
! -------------------------------------------------------------------------
! --- floats.f (hycom version 2.2)
! -------------------------------------------------------------------------
!
! --- synthetic floats, drifters, and mooring instruments
!
! --- optionally samples time series of dynamical/thermodynamical variables
! --- at the location of each float
!
! --- four float types are presently supported:
! --- 3-d lagrangian (vertically advected by diagnosed vertical velocity)
! --- isopycnic (remains on specified density surface)
! --- isobaric (remains at the pressure depth where it was released)
! --- stationary (synthetic instrument/mooring)
!
! --- horizontal advection, along with vertical advection of lagrangian floats,
! --- is performed using the MICOM runga-kutta-4 time interpolation algorithm
! --- developed by Zulema Garraffo
!
! --- works on any curvilinear grid
!
! --- float position is stored as longitude and latitude
!
! --- to horizontal advect the float, u and v are first converted to
! --- d(longitude)/dt and d(latitude)/dt as follows:
!
! --- u_lon = u*dlondx + v*dlondy
! --- v_lat = u*dlatdx + v*dlatdy
!
! --- to horizontally advect the float, terms u*dlondx and u*dlatdx are
! --- calculated on u grid points while terms v*dlondy and v*dlatdy are
! --- calculated on v grid ponts. these terms are then interpolated to
! --- the float locations from the surrounding u and v grid points,
! --- respectively
!
! --- since the terms v*dlondy and u*dlatdx are always zero when the model
! --- x,y axes are everywhere lines of constant latitude and longitude,
! --- respectively, logical variable 'nonlatlon' prevents horizontal
! --- interpolation of these terms in this case
!
! --- horizontal interpolation to the float location follows the MICOM
! --- procedure of Zulema Garraffo - second-order interpolation from the 16
! --- surrounding grid points is performed unless fewer than nptmin (presently
! --- set to 10 in subroutine intrph) water points are available, in which
! --- case nearest-neighbor interpolation is used.
!
! --- variables are horizontally interpolated from their native grid (p, u,
! --- or v)
!
! --- adapted for hycom by george halliwell
! --- code parallelized by remy baraille
!
! --- the second-order interpolation subroutine included here is the one
! --- included in the mariano and brown parameter matrix objective analysis
! --- algorithm to estimate the large scale trend surface - it is different
! --- from the algorithm used by Garraffo in MICOM
!
! --- float initialization is performed in floats_init - information
! --- about the input file containing initial float information is
! --- presented in that subroutine.
!
! --- variables in array flt(nfl,n) are:
!
! --- n = 1 longitude
! --- n = 2 latitude
! --- n = 3 float depth
! --- n = 4 water depth
! --- n = 5 temperature
! --- n = 6 salinity
! --- n = 7 water density
! --- n = 8 float start time (in days from start of model run)
! --- n = 9 float end time (in days from start of model run)
!
! --- time series of several variables are interpolated to the location of each
! --- float and saved every nflsam time steps when ioflag is set to 1
!
! --- variables output onto file float.out for each float are:
!
! --- 1. initial sequential float number
! --- 2. time (in days from start of model run)
! --- 3. model layer number that contains the float
! --- 4. longitude (u for synthetic moorings)
! --- 5. latitude (v for synthetic moorings)
! --- 6. float depth (w for synthetic moorings)
! --- 7. water depth
! --- 8. temperature
! --- 9. salinity
! --- 10. water density (remains constant for isopycnic floats)
!
call xctilr(dp( 1-nbdy,1-nbdy,1,n),1,kk, 6,6, halo_ps)
call xctilr(dpu( 1-nbdy,1-nbdy,1,n),1,kk, 6,6, halo_vs)
call xctilr(dpv( 1-nbdy,1-nbdy,1,n),1,kk, 6,6, halo_vs)
call xctilr(u( 1-nbdy,1-nbdy,1,n),1,kk, 6,6, halo_uv)
call xctilr(v( 1-nbdy,1-nbdy,1,n),1,kk, 6,6, halo_vv)
call xctilr(ubavg( 1-nbdy,1-nbdy, n),1, 1, 6,6, halo_uv)
call xctilr(vbavg( 1-nbdy,1-nbdy, n),1, 1, 6,6, halo_vv)
call xctilr(temp( 1-nbdy,1-nbdy,1,n),1,kk, 6,6, halo_ps)
call xctilr(saln( 1-nbdy,1-nbdy,1,n),1,kk, 6,6, halo_ps)
call xctilr(th3d( 1-nbdy,1-nbdy,1,n),1,kk, 6,6, halo_ps)
!
! --- calculate vertical velocity at interfaces as the sum of the
! --- vertical interface velocity estimated in cnuity.f and the
! --- advective component due to flow parallel to interfaces.
! --- wvelup and wveldn represent velocity at the top and bottom of
! --- layer k
!
margin = 6
!
! --- set pressure array at p points
do j=1-margin,jj+margin
do i=1-margin,ii+margin
if (SEA_P) then
do k=1,kk
k1=k+1
p(i,j,k+1)=p(i,j,k)+dp(i,j,k,n)
enddo !k
endif !ip
enddo !i
enddo !j
!
! --- set pressure array at u and v points
! --- calculate dpdx, dpdy at interfaces above and below layer k for
! --- estimating wvel within layer k
!
do k=1,kk
!
margin = 5
!
do j=1-margin,jj+margin
do i=1-margin,ii+margin
if (SEA_U) then
pu(i,j,k+1)=pu(i,j,k)+dpu(i,j,k,n)
if (wvelfl) then
dpdxup(i,j)= &
(p(i,j,k )*scpy(i,j)-p(i-1,j,k )*scpy(i-1,j)) &
/scu2(i,j)
dpdxdn(i,j)= &
(p(i,j,k+1)*scpy(i,j)-p(i-1,j,k+1)*scpy(i-1,j)) &
/scu2(i,j)
endif
endif !iu
if (SEA_V) then
pv(i,j,k+1)=pv(i,j,k)+dpv(i,j,k,n)
if (wvelfl) then
dpdyup(i,j)= &
(p(i,j,k )*scpx(i,j)-p(i,j-1,k )*scpx(i,j-1)) &
/scv2(i,j)
dpdydn(i,j)= &
(p(i,j,k+1)*scpx(i,j)-p(i,j-1,k+1)*scpx(i,j-1)) &
/scv2(i,j)
endif
endif !iv
enddo !i
enddo !j
!
! --- calculate the vertical velocity arrays
!
margin = 4
!
do j=1-margin,jj+margin
do i=1-margin,ii+margin
if (SEA_P) then
if (wvelfl) then
wveli(i,j,k+1)=deltfl*wveli(i,j,k+1)/delt1
if (dp(i,j,k,n).gt.tencm) then
wvelup(i,j,k)=-0.5*deltfl* &
#if defined(STOKES)
((u(i ,j ,k,n)+usd(i ,j ,k)+ &
ubavg(i ,j ,n))*dpdxup(i ,j )+ &
(u(i+1,j ,k,n)+usd(i+1,j ,k)+ &
ubavg(i+1,j ,n))*dpdxup(i+1,j )+ &
(v(i ,j ,k,n)+vsd(i ,j ,k)+ &
vbavg(i ,j ,n))*dpdyup(i ,j )+ &
(v(i ,j+1,k,n)+vsd(i ,j+1,k)+ &
vbavg(i ,j+1,n))*dpdyup(i ,j+1)) &
#else
((u(i ,j ,k,n)+ubavg(i ,j ,n))*dpdxup(i ,j )+ &
(u(i+1,j ,k,n)+ ubavg(i+1,j ,n))*dpdxup(i+1,j )+ &
(v(i ,j ,k,n)+vbavg(i ,j ,n) )*dpdyup(i ,j )+ &
(v(i ,j+1,k,n)+vbavg(i ,j+1,n))*dpdyup(i ,j+1)) &
#endif
/onem+wveli(i,j,k )
wveldn(i,j,k)=-0.5*deltfl* &
#if defined(STOKES)
((u(i ,j ,k,n)+usd(i ,j ,k)+ &
ubavg(i ,j ,n))*dpdxdn(i ,j )+ &
(u(i+1,j ,k,n)+usd(i+1,j ,k)+ &
ubavg(i+1,j,n))*dpdxdn(i+1,j )+ &
(v(i ,j ,k,n)+vsd(i ,j ,k)+ &
vbavg(i ,j ,n))*dpdydn(i ,j )+ &
(v(i ,j+1,k,n)+vsd(i ,j+1,k)+ &
vbavg(i ,j+1,n))*dpdydn(i ,j+1)) &
#else
((u(i ,j ,k,n)+ubavg(i ,j ,n))*dpdxdn(i ,j )+ &
(u(i+1,j ,k,n)+ubavg(i+1,j ,n))*dpdxdn(i+1,j )+ &
(v(i ,j ,k,n)+vbavg(i ,j ,n))*dpdydn(i ,j )+ &
(v(i ,j+1,k,n)+vbavg(i ,j+1,n))*dpdydn(i ,j+1)) &
#endif
/onem+wveli(i,j,k+1)
else
wvelup(i,j,k)=0.0
wveldn(i,j,k)=0.0
endif
endif
endif !ip
enddo !i
enddo !j
!
enddo !k
!
! call xctilr(wvelup(1-nbdy,1-nbdy,1),1, kk, 4,4, halo_ps)
! call xctilr(wveldn(1-nbdy,1-nbdy,1),1, kk, 4,4, halo_ps)
!
! --- set old velocity indices used for interpolation of float position
! --- perform full float update every two times this subroutine is called;
! --- at the intermediate times, just store old velocity fields for the next
! --- call
!
if (mod(nstepfl,nfladv).eq.nfldta) then
kold1 =kk
kold1w=kk+1
kold2 =0
kold2w=0
go to 10
else
kold1 =0
kold1w=0
kold2 =kk
kold2w=kk+1
endif
!
! --- turbulent horizontal velocity option
!
! --- to customize the parameter settings in blkdat.input that control
! --- the calculated turbulent velocity (tbvar, tdecri), refer to
! --- Griffa (1996), "Aplications of stochastic particles models to
! --- oceanographic problems" in "Stochatic Modelling in Physical
! --- Oceanography", pp. 114-140, Adler, Muller, and Rozovoskii, editors
!
if (turbvel) then
!
! --- initialize random seeds and create random number table
! --- if iflag=1 (iflag=0), then time() is (is not) used to generate seeds
!
iflag=1
call system_clock(inisee)
seed = 414957000-inisee
iede = iseed1
ied = iseed2
if (iede .lt. 0) iede = abs(iede)
if (ied .lt. 0) ied = abs(ied)
!
call rantab_ini(rtab,iseed1,iseed2,iflag)
endif
!
! --- get particle locations from processor 1
if (mnproc.eq.1) then
do nfl=1,nflt
do i=1,9
flt1( i+12*(nfl-1)) = flt(nfl,i)
enddo
flt1(10+12*(nfl-1)) = kfloat(nfl)
flt1(11+12*(nfl-1)) = iflnum(nfl)
flt1(12+12*(nfl-1)) = ifltyp(nfl)
enddo
endif !1st tile
call xcastr(flt1(1:12*nflt), 1)
if (mnproc.ne.1) then
do nfl=1,nflt
do i=1,9
flt(nfl,i) = flt1( i+12*(nfl-1))
enddo
kfloat(nfl) = flt1(10+12*(nfl-1))
iflnum(nfl) = flt1(11+12*(nfl-1))
ifltyp(nfl) = flt1(12+12*(nfl-1))
enddo
endif !1st tile
!
! ----------------------
! --- begin float loop
! ----------------------
!
do nfl=1,nflt
!
! --- skip if float has previously run aground or exceeded termination time
if(flt(nfl,1).le.-999.) then
go to 22
endif
!
! --- ier = error flag for horizontal interpolation
! --- ntermn = float termination flag
! --- -10 => reached termination time
! --- -5 => exited domain
! --- >0 => ran aground
! --- ifladv = float horizontal advection flag
ier=0
ntermn=0
ifladv=1
!
! --- suppress float advection if this is the first time step so that
! --- the initial position on the output file is exactly the position
! --- specified on the input file
if (nstepfl.eq.1) then
ifladv=0
endif
!
! --- check if model time has reached float deployment time
! --- when it does, again suppress float advection during first pass
if (timefl-flt(nfl,8).lt.-0.001) then
go to 22
endif
if (kfloat(nfl).lt.0) then
ifladv=0
endif
!
! --- check if model time has reached float termination time
if ( flt(nfl,9).gt.0.0 .and. &
timefl-flt(nfl,9).gt.0.001 ) then
ntermn=-10
endif
!
! ---------------------------------------------------------------------------
! --- search for the processor controlling the tile containing the grid point
! --- search for the surrounding 16 grid points on the p, u, and v grids
! ---------------------------------------------------------------------------
!
! --- processeur sur lequel va se derouler le calcul
!
iproc=0
!
margin = 0
!
do j=1-margin,jj+margin
do i=1-margin,ii+margin
if (SEA_P) then
!
! --- search for the surrounding p-grid points
if (i+i0.lt.itdm .and. j+j0.lt.jtdm) then
alomin=min(plon(i ,j),plon(i ,j+1), &
plon(i+1,j),plon(i+1,j+1))
alomax=max(plon(i ,j),plon(i ,j+1), &
plon(i+1,j),plon(i+1,j+1))
alamin=min(plat(i ,j),plat(i ,j+1), &
plat(i+1,j),plat(i+1,j+1))
alamax=max(plat(i ,j),plat(i ,j+1), &
plat(i+1,j),plat(i+1,j+1))
if (flt(nfl,1).ge.alomin.and.flt(nfl,1).lt.alomax.and. &
flt(nfl,2).ge.alamin.and.flt(nfl,2).lt.alamax) then
ifltll(1)=i
jfltll(1)=j
!
! --- determine if float has exited domain
if (i+i0.lt.2 .or. i+i0.gt.itdm-1 .or. &
j+j0.lt.2 .or. j+j0.gt.jtdm-1) then
ntermn=-5
endif
!
! --- search for the surrounding u-grid points
do i1=i-1,i+1
do j1=j-1,j+1
if (i1+i0.lt.itdm .and. j1+j0.lt.jtdm) then
alomin=min(ulon(i1 ,j1),ulon(i1 ,j1+1), &
ulon(i1+1,j1),ulon(i1+1,j1+1))
alomax=max(ulon(i1 ,j1),ulon(i1 ,j1+1), &
ulon(i1+1,j1),ulon(i1+1,j1+1))
alamin=min(ulat(i1 ,j1),ulat(i1 ,j1+1), &
ulat(i1+1,j1),ulat(i1+1,j1+1))
alamax=max(ulat(i1 ,j1),ulat(i1 ,j1+1), &
ulat(i1+1,j1),ulat(i1+1,j1+1))
if (flt(nfl,1).ge.alomin .and. &
flt(nfl,1).lt.alomax .and. &
flt(nfl,2).ge.alamin .and. &
flt(nfl,2).lt.alamax) then
ifltll(2)=i1
jfltll(2)=j1
!
! --- determine if float has exited domain
if (i1+i0.lt.2 .or. i1+i0.gt.itdm-1 .or. &
j1+j0.lt.2 .or. j1+j0.gt.jtdm-1) then
ntermn=-5
endif
go to 11
endif
endif
enddo !j1
enddo !i1
11 continue
!
! --- search for the surrounding v-grid points
do i1=i-1,i+1
do j1=j-1,j+1
if (i1+i0.lt.itdm .and. j1+j0.lt.jtdm) then
alomin=min(vlon(i1 ,j1),vlon(i1 ,j1+1), &
vlon(i1+1,j1),vlon(i1+1,j1+1))
alomax=max(vlon(i1 ,j1),vlon(i1 ,j1+1), &
vlon(i1+1,j1),vlon(i1+1,j1+1))
alamin=min(vlat(i1 ,j1),vlat(i1 ,j1+1), &
vlat(i1+1,j1),vlat(i1+1,j1+1))
alamax=max(vlat(i1 ,j1),vlat(i1 ,j1+1), &
vlat(i1+1,j1),vlat(i1+1,j1+1))
if (flt(nfl,1).ge.alomin .and. &
flt(nfl,1).lt.alomax .and. &
flt(nfl,2).ge.alamin .and. &
flt(nfl,2).lt.alamax) then
ifltll(3)=i1
jfltll(3)=j1
!
! --- determine if float has exited domain
if (i1+i0.lt.2 .or. i1+i0.gt.itdm-1 .or. &
j1+j0.lt.2 .or. j1+j0.gt.jtdm-1) then
ntermn=-5
endif
go to 12
endif
endif
enddo !j1
enddo !i1
12 continue
!
! --- set processor number for the tile containing the float and exit grid
! --- point search
iproc=mnproc
go to 13
endif
endif
endif !ip
enddo !i
enddo !j
!
13 continue
!
! --- float nfl is now updated by the processor running the tile containing
! --- the float
!
if (iproc.gt.0) then
!
if (nfl.eq.nfl_debug) then
write(lp,100) nfl,nflt,nstep,nstepfl
write(lp,101) nfl,ntermn,flt(nfl,1),flt(nfl,2), &
(ifltll(i)+i0,jfltll(i)+j0,i=1,3), &
mnproc, &
plon(ifltll(1),jfltll(1)), &
plat(ifltll(1),jfltll(1))
100 format(/'diagnostics for float',i6,' of',i6,', time step', &
i9/'float time step',i9)
101 format('float',i6,' ntermn:',i6,' position:',2(1pe12.4)/ &
'lower left points, p,u,v:',3(2i5,2x)/ &
'mnproc =',i4,' plon,plat =',1pe12.4,1pe12.4)
call flush(lp)
endif !nfl_debug
!
! --- if float has exited domain or run aground as determined previously,
! --- jump ahead to the float termination code bloci
if (ntermn.eq.-5 .or. ntermn.eq.1) then
go to 30
endif
!
! --- set ptlon, ptlat for all horizontal interpolations from each grid
ngrid=1
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
ptlon(i1,j1,ngrid)=plon(i,j)
ptlat(i1,j1,ngrid)=plat(i,j)
enddo
enddo
!
ngrid=2
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
ptlon(i1,j1,ngrid)=ulon(i,j)
ptlat(i1,j1,ngrid)=ulat(i,j)
enddo
enddo
!
ngrid=3
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
ptlon(i1,j1,ngrid)=vlon(i,j)
ptlat(i1,j1,ngrid)=vlat(i,j)
enddo
enddo
!
! --- the float is assumed to remain within the same layer that it was in
! --- during the previous update unless it is the first advection time step
! --- for the float, in which case it is initially assumed to be in layer 1
!
k=max(1,kfloat(nfl))
!
! --- mask p grid points that are on land or where the layer containing the
! --- float is a zero or nearly-zero thickness layer at the bottom
ngrid=1
ngood(ngrid)=0
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
if (depths(i,j).lt.0.01 .or. depths(i,j).eq.hugel .or. &
depths(i,j)*onem-p(i,j,k).lt.tencm) then
maskpt(i1,j1,ngrid)=.true.
else
maskpt(i1,j1,ngrid)=.false.
ngood(ngrid)=ngood(ngrid)+1
endif
enddo
enddo
!
if (nfl.eq.nfl_debug) then
write(lp,102) ngood(1)
102 format('initial masking: no. of good p points',i4)
endif !nfl_debug
!
! --------------------------------------------------
! --- determine the model layer containing the float
! --------------------------------------------------
!
! --- if this is not the first advection time step for the float, determine
! --- if the float has transited to another model layer
!
if (kfloat(nfl).ge.1) then
!
! --- check if the float is deeper than the interpolated bottom
! --- depth - if so, reset the depth to 0.1 m above the bottom - this
! --- prevents from running aground too frequently in coastal regions
if (ifltyp(nfl).ne.4) then
ngrid=1
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
varb2d(i1,j1)=p(i,j,kk+1)
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
flt(nfl,1),flt(nfl,2),maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,plo,ier)
flt(nfl,3)=max(fldepm,min(flt(nfl,3),(plo-tencm)/onem))
endif
!
! --- has the float remained in the same layer?
ngrid=1
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
varb2d(i1,j1)=p(i,j,k)
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
flt(nfl,1),flt(nfl,2),maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,phi,ier)
!
! --- terminate float because it has run aground
if (ier.eq.999) then
ntermn=2
go to 30
endif
!
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
varb2d(i1,j1)=p(i,j,k+1)
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
flt(nfl,1),flt(nfl,2),maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,plo,ier)
phi=phi/onem
plo=plo/onem
!
! --- if the float is in the same model layer, skip ahead
if (flt(nfl,3).gt.phi .and. flt(nfl,3).le.plo) then
kflt=kfloat(nfl)
go to 40
endif
!
! --- determine the new model layer containing the float
!
! --- did the float move up one layer?
if (flt(nfl,3).le.phi .and. kfloat(nfl).gt.1) then
k=kfloat(nfl)-1
plo=phi
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
varb2d(i1,j1)=p(i,j,k)
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
flt(nfl,1),flt(nfl,2),maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,phi,ier)
phi=phi/onem
if (flt(nfl,3).gt.phi .and. flt(nfl,3).le.plo) then
kflt=k
go to 40
endif
endif
!
! --- did the float move down one layer?
if (flt(nfl,3).gt.plo .and. kfloat(nfl).lt.kk) then
k=kfloat(nfl)+1
phi=plo
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
varb2d(i1,j1)=p(i,j,k)
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
flt(nfl,1),flt(nfl,2),maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,plo,ier)
plo=plo/onem
if (flt(nfl,3).gt.phi .and. flt(nfl,3).le.plo) then
kflt=k
go to 40
endif
endif
!
endif
!
! --- if the float did not move up or down one layer, or if this is the
! --- first time step after float initialization, search for the layer
! --- containg the float from top to bottom
ngrid=1
plo=0.
k2=0
k3=0
do k=1,kk
k1=min(k+3,kk+1)
!
! --- if the depth of interface k1 is not deeper than the float at any
! --- surrounding grid points, skip to end of k-loop
if (k2.eq.0) then
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
if (.not.maskpt(i1,j1,ngrid) .and. &
p(i,j,k1)/onem.gt.flt(nfl,3)) then
k2=1
endif
enddo
enddo
endif
if (k2.eq.0) then
go to 42
endif
!
! --- mask points where the layer being tested rests on the bottom
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
if (k2.eq.1) then
maskpi(i1,j1)=maskpt(i1,j1,ngrid)
endif
if (.not.maskpi(i1,j1) .and. &
depths(i,j)*onem-p(i,j,k).lt.tencm) then
maskpi(i1,j1)=.true.
ngood(ngrid)=ngood(ngrid)-1
endif
enddo
enddo
k2=2
!
! --- if too few good points remain, float has run aground
if (ngood(ngrid).le.2) then
ntermn=3
go to 30
endif
!
! --- interpolate to find plo
phi=plo
ngrid=1
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
varb2d(i1,j1)=p(i,j,k+1)
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
flt(nfl,1),flt(nfl,2),maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,plo,ier)
plo=plo/onem
if(plo-phi.gt.0.1) then
k3=k
else
kflt=k3
do i1=1,4
do j1=1,4
maskpt(i1,j1,ngrid)=maskpi(i1,j1)
enddo
enddo
go to 40
endif
if (flt(nfl,3).lt.plo) then
kflt=k
do i1=1,4
do j1=1,4
maskpt(i1,j1,ngrid)=maskpi(i1,j1)
enddo
enddo
go to 40
endif
42 continue
enddo
!
! --- if layer selection fails, assume float has run aground
print *,'warning - selection of float layer failed'
ntermn=4
go to 30
!
! --- we now know the model layer containing the float
40 kfloat(nfl)=kflt
!
! ------------------------------------------
! --- set grid masks for u and v grid points
! ------------------------------------------
!
k=kflt
ngrid=2
ngood(ngrid)=0
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
if (depthu(i,j).lt.onecm .or. depthu(i,j).eq.hugel .or. &
depthu(i,j)-pu(i,j,k).lt.tencm) then
maskpt(i1,j1,ngrid)=.true.
else
maskpt(i1,j1,ngrid)=.false.
ngood(ngrid)=ngood(ngrid)+1
endif
enddo
enddo
!
ngrid=3
ngood(ngrid)=0
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
if (depthv(i,j).lt.onecm .or. depthv(i,j).eq.hugel .or. &
depthv(i,j)-pv(i,j,k).lt.tencm) then
maskpt(i1,j1,ngrid)=.true.
else
maskpt(i1,j1,ngrid)=.false.
ngood(ngrid)=ngood(ngrid)+1
endif
enddo
enddo
!
! --- if too few good p, u, v grid points, assume float has run aground
if (ngood(1).le.2 .or. ngood(2).le.2 .or. ngood(3).le.2) then
ntermn=5
go to 30
endif
!
if (nfl.eq.nfl_debug) then
write(lp,103) nfl,k,phi,flt(nfl,3),plo,q
write(lp,104) ngood
103 format('nfl,k,phi,pflt,plo,q',i6,i3,1p,4e12.4)
104 format('number of good points, p,u,v',3i4)
endif !nfl_debug
!
! --- set vertical indices, and also set q for vertical interpolation
k=kfloat(nfl)
k1=max(1,k-1)
k2=min(kk,k+1)
q=(plo-flt(nfl,3))/max(0.001,plo-phi)
!
! ----------------------------------------------------------
! --- advect the float horizontally, then move it vertically
! --- do not move float if it is a mooring instrument
! ----------------------------------------------------------
!
! ----------------------------------
! --- horizontal advection of floats
! ----------------------------------
!
if (ifltyp(nfl).ne.4) then
!
! --- interpolate u,v to float location - velocities at the new time and
! --- at two earlier times are used for the runga-kutta time interpolation
!
xpos0=flt(nfl,1)
ypos0=flt(nfl,2)
!
ngrid=2
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
varb2d(i1,j1)=dlondx(i,j)*uold2(i,j,k+kold2)
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
xpos0,ypos0,maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,uflt1,ier)
!
! --- terminate float because it has run aground
if (ier.eq.999) then
ntermn=6
go to 30
endif
!
if (nonlatlon) then
ngrid=3
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
varb2d(i1,j1)=dlondy(i,j)*vold2(i,j,k+kold2)
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
xpos0,ypos0,maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,ufltm,ier)
uflt1=uflt1+ufltm
!
! --- terminate float because it has run aground
if (ier.eq.999) then
ntermn=7
go to 30
endif
endif
!
ngrid=3
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
varb2d(i1,j1)=dlatdy(i,j)*vold2(i,j,k+kold2)
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
xpos0,ypos0,maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,vflt1,ier)
!
! --- terminate float because it has run aground
if (ier.eq.999) then
ntermn=8
go to 30
endif
!
if (nonlatlon) then
ngrid=2
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
varb2d(i1,j1)=dlatdx(i,j)*uold2(i,j,k+kold2)
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
xpos0,ypos0,maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,vfltm,ier)
vflt1=vflt1+vfltm
!
! --- terminate float because it has run aground
if (ier.eq.999) then
ntermn=9
go to 30
endif
endif
!
xpos1=flt(nfl,1)+uflt1
ypos1=flt(nfl,2)+vflt1
!
ngrid=2
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
varb2d(i1,j1)=dlondx(i,j)*uold2(i,j,k+kold1)
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
xpos1,ypos1,maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,uflt2,ier)
!
if (nonlatlon) then
ngrid=3
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
varb2d(i1,j1)=dlondy(i,j)*vold2(i,j,k+kold1)
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
xpos0,ypos0,maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,ufltm,ier)
uflt2=uflt2+ufltm
endif
!
ngrid=3
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
varb2d(i1,j1)=dlatdy(i,j)*vold2(i,j,k+kold1)
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
xpos1,ypos1,maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,vflt2,ier)
!
if (nonlatlon) then
ngrid=2
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
varb2d(i1,j1)=dlatdx(i,j)*uold2(i,j,k+kold1)
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
xpos0,ypos0,maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,vfltm,ier)
vflt2=vflt2+vfltm
endif
!
xpos2=flt(nfl,1)+uflt2
ypos2=flt(nfl,2)+vflt2
ngrid=2
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
varb2d(i1,j1)=dlondx(i,j)*uold2(i,j,k+kold1)
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
xpos2,ypos2,maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,uflt3,ier)
!
if (nonlatlon) then
ngrid=3
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
varb2d(i1,j1)=dlondy(i,j)*vold2(i,j,k+kold1)
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
xpos0,ypos0,maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,ufltm,ier)
uflt3=uflt3+ufltm
endif
!
ngrid=3
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
varb2d(i1,j1)=dlatdy(i,j)*vold2(i,j,k+kold1)
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
xpos2,ypos2,maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,vflt3,ier)
!
if (nonlatlon) then
ngrid=2
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
varb2d(i1,j1)=dlatdx(i,j)*uold2(i,j,k+kold1)
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
xpos0,ypos0,maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,vfltm,ier)
vflt3=vflt3+vfltm
endif
!
endif
!
if (ifltyp(nfl).ne.4) then
xpos3=flt(nfl,1)+2.0*uflt3
ypos3=flt(nfl,2)+2.0*vflt3
else
xpos3=flt(nfl,1)
ypos3=flt(nfl,2)
endif
!
ngrid=2
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
if (ifltyp(nfl).ne.4) then
#if defined(STOKES)
varb2d(i1,j1)=dlondx(i,j)*(u(i,j,k,n)+usd(i,j,k)+ &
ubavg(i,j,n))
else
varb2d(i1,j1)=u(i,j,k,n)+usd(i,j,k)+ubavg(i,j,n)
#else
varb2d(i1,j1)=dlondx(i,j)*(u(i,j,k,n)+ubavg(i,j,n))
else
varb2d(i1,j1)=u(i,j,k,n)+ubavg(i,j,n)
#endif
endif
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
xpos3,ypos3,maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,uflt4,ier)
!
if (nonlatlon .and. ifltyp(nfl).ne.4) then
ngrid=3
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
#if defined(STOKES)
varb2d(i1,j1)=dlondy(i,j)*(v(i,j,k,n)+vsd(i,j,k)+ &
vbavg(i,j,n))
#else
varb2d(i1,j1)=dlondy(i,j)*(v(i,j,k,n)+vbavg(i,j,n))
#endif
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
xpos0,ypos0,maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,ufltm,ier)
uflt4=uflt4+ufltm
endif
!
ngrid=3
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
if (ifltyp(nfl).ne.4) then
#if defined(STOKES)
varb2d(i1,j1)=dlatdy(i,j)*(v(i,j,k,n)+vsd(i,j,k)+ &
vbavg(i,j,n))
else
varb2d(i1,j1)=v(i,j,k,n)+vsd(i,j,k)+vbavg(i,j,n)
#else
varb2d(i1,j1)=dlatdy(i,j)*(v(i,j,k,n)+vbavg(i,j,n))
else
varb2d(i1,j1)=v(i,j,k,n)+vbavg(i,j,n)
#endif
endif
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
xpos3,ypos3,maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,vflt4,ier)
!
if (nonlatlon .and. ifltyp(nfl).ne.4) then
ngrid=2
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
#if defined(STOKES)
varb2d(i1,j1)=dlatdx(i,j)*(u(i,j,k,n)+usd(i,j,k)+ &
ubavg(i,j,n))
#else
varb2d(i1,j1)=dlatdx(i,j)*(u(i,j,k,n)+ubavg(i,j,n))
#endif
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
xpos0,ypos0,maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,vfltm,ier)
vflt4=vflt4+vfltm
endif
!
if (ifltyp(nfl).ne.4) then
uflt=(uflt1+2.0*uflt2+2.0*uflt3+uflt4)/6.0
vflt=(vflt1+2.0*vflt2+2.0*vflt3+vflt4)/6.0
!
! --- add turbulent velocity to u and v if requested
!
if (turbvel) then
call rantab(rtab,iseed1,iseed2,numran)
uturb=uturb0*rtab(numran, 2)
call rantab(rtab,iseed1,iseed2,numran)
vturb=uturb0*rtab(numran, 2)
!
! --- convert uturb, vturb to d(longitude)/dt, d(latitude)/dt
! --- interpolate dlondx, dlondy, dlatdx, dlatdy to the float location
!
ngrid=2
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
varb2d(i1,j1)=dlondx(i,j)
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
xpos0,ypos0,maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,dlodx,ier)
if (nonlatlon) then
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
varb2d(i1,j1)=dlatdx(i,j)
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
xpos0,ypos0,maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,dladx,ier)
ngrid=3
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
varb2d(i1,j1)=dlondy(i,j)
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
xpos0,ypos0,maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,dlody,ier)
else
dladx=0.0
dlody=0.0
endif
ngrid=3
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
varb2d(i1,j1)=dlatdy(i,j)
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
xpos0,ypos0,maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,dlady,ier)
!
! --- add turbulent velocity to uflt, vflt
!
uturb1=uturb*dlodx+vturb*dlody
vturb1=uturb*dladx+vturb*dlady
!
if (nfl.eq.nfl_debug) then
write(lp,105) nfl,uflt,uturb1,vflt,vturb1
105 format('nfl,uflt,uturb1,vflt,vturb1',i5,1p,4e12.4)
endif !nfl_debug
!
uflt=(1.0-dtturb*tdecri)*uflt+uturb1
vflt=(1.0-dtturb*tdecri)*vflt+vturb1
!
endif
!
! --- update float horizontal position
!
if (ifladv.eq.1) then
flt(nfl,1)=flt(nfl,1)+2.0*deltfl*uflt
flt(nfl,2)=flt(nfl,2)+2.0*deltfl*vflt
endif
!
if (nfl.eq.nfl_debug) then
write(lp,106) nfl,flt(nfl,1),xpos1,xpos2,xpos3, &
uflt1,uflt2,uflt3,uflt4,2.0*deltfl*uflt
write(lp,107) nfl,flt(nfl,2),ypos1,ypos2,ypos3, &
vflt1,vflt2,vflt3,vflt4,2.0*deltfl*vflt
write(lp,108) nfl,uflt,vflt,flt(nfl,1),flt(nfl,2)
106 format('nfl,x-position',i6,1p,4e12.4/' u-velocity',6x,5e12.4)
107 format('nfl,y-position',i6,1p,4e12.4/' v-velocity',6x,5e12.4)
108 format('nfl,udeg,vdeg',i6,1p,2e13.5/'new position',2e13.5)
endif !nfl_debug
!
endif
!
! -----------------------------------------------
! --- vertical advection of 3-d lagrangian floats
! -----------------------------------------------
!
! --- horizontally interpolate diagnosed vertical velocity to the float
! --- location at the new time and at two earlier times to execute the
! --- runga-kutta time interpolation
!
if (wvelfl .and. ifltyp(nfl).eq.1) then
!
ngrid=1
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
wvhi=wold2u(i,j,k+kold2)
wvlo=wold2d(i,j,k+kold2)
varb2d(i1,j1)=q*wvhi+(1.0-q)*wvlo
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
xpos0,ypos0,maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,wflt1,ier)
!
ngrid=1
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
wvhi=wold2u(i,j,k+kold1)
wvlo=wold2d(i,j,k+kold1)
varb2d(i1,j1)=q*wvhi+(1.0-q)*wvlo
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
xpos1,ypos1,maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,wflt2,ier)
!
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
xpos2,ypos2,maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,wflt3,ier)
!
endif
!
if (wvelfl) then
ngrid=1
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
wvhi=wvelup(i,j,k)
wvlo=wveldn(i,j,k)
varb2d(i1,j1)=q*wvhi+(1.0-q)*wvlo
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
xpos3,ypos3,maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,wflt4,ier)
endif
!
if (ifltyp(nfl).eq.1) then
!
! --- dividing by 3 gives total vertical displacement over interval 2*deltfl
wflt=(wflt1+2.0*wflt2+2.0*wflt3+wflt4)/3.0
if (ifladv.eq.1) then
flt(nfl,3)=max(fldepm,flt(nfl,3)-wflt)
endif
if (nfl.eq.nfl_debug) then
write(lp,109) nfl,wflt,flt(nfl,3)
109 format('nfl,w,new depth',i6,1p,2e13.5)
endif !nfl_debug
!
elseif (ifltyp(nfl).eq.2) then
!
! ---------------------------------------
! --- vertical motion of isopycnic floats
! ---------------------------------------
!
! --- set float depth to the vertically interpolated depth of the specified
! --- density interface
!
! --- assume float runs aground if no water denser than the specified density
! --- exists in the water column
!
! --- limit float depth to fldepm if no water lighter than the specified
! --- density exists in the water column
!
! --- find smallest k where float density is exceeded by the density in
! --- layer k at at least one of the good 16 surrounding grid points
thflt=flt(nfl,7)-thbase
ngrid=1
ngoodi=ngood(ngrid)
do k=1,kk
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
maskpi(i1,j1)=maskpt(i1,j1,ngrid)
if (.not.maskpi(i1,j1)) then
if (th3d(i,j,k,n).gt.thflt) then
kflt=k
go to 33
endif
endif
enddo
enddo
enddo
33 continue
!
! --- move downward from layer kflt to find the two layers with densities
! --- bracketing the assigned float density
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
if (.not.maskpi(i1,j1)) then
if (kflt.eq.1) then
phi=0.
plo=.5*dp(i,j,kflt,n)/onem
thflt=flt(nfl,7)-thbase
if(thflt.le.th3d(i,j,kflt,n)) then
varb2d(i1,j1)=max(fldepm,plo)
go to 50
endif
endif
do k=max(2,kflt),kk
phi=p(i,j,k-1)/onem
plo=p(i,j,k )/onem
thhi=th3d(i,j,k-1,n)
thlo=th3d(i,j,k ,n)
if(thflt.gt.thhi .and. thflt.le.thlo) then
qisop=max(0.0,min(1.0,(thlo-thflt)/ &
max(epsil,thlo-thhi)))
varb2d(i1,j1)=max(fldepm,qisop*phi+(1.-qisop)*plo)
if (depths(i,j)-varb2d(i1,j1).lt.0.1) then
maskpi(i1,j1)=.true.
ngoodi=ngoodi-1
endif
go to 50
endif
enddo
50 continue
endif
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
flt(nfl,1),flt(nfl,2),maskpi, &
ngoodi,ngrid,intpfl,radian,depiso,ier)
flt(nfl,3)=max(fldepm,depiso)
if (nfl.eq.nfl_debug) then
write(lp,110) nfl,flt(nfl,3)
110 format('nfl,new isopycnic depth',i6,1p,e13.5)
endif !nfl_debug
!
elseif (ifltyp(nfl).eq.4) then
!
wflt=wflt4/deltfl
!
endif
!
! --------------------------------------------------------------
! --- interpolate water properties to the present float location
! --------------------------------------------------------------
!
if (ioflag.eq.1 .and. samplfl) then
!
! --- recalculate q. assume float does not leave the layer that it was
! --- originally diagnosed to be within.
!
q=max(0.0,min(1.0,(plo-flt(nfl,3))/max(0.001,plo-phi)))
!
! --- water depth interpolation
if (ifltyp(nfl).ne.4 .or. flt(nfl,4).eq.0.0) then
ngrid=1
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
varb2d(i1,j1)=depths(i,j)
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
flt(nfl,1),flt(nfl,2),maskpt(1,1,ngrid), &
16,ngrid,intpfl,radian,depflt,ier)
flt(nfl,4)=depflt
endif
!
! --- temperature interpolation
ngrid=1
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
if (ifltyp(nfl).eq.4) then
k1=max(1 ,k-1)
k2=min(kk,k+1)
varb2d(i1,j1)= &
0.5*( q *(temp(i,j,k1,n)+temp(i,j,k,n))+ &
(1.0-q)*(temp(i,j,k2,n)+temp(i,j,k,n)))
else
varb2d(i1,j1)=temp(i,j,k,n)
endif
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
flt(nfl,1),flt(nfl,2),maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,tflt,ier)
flt(nfl,5)=tflt
!
! --- salinity interpolation
ngrid=1
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
if (ifltyp(nfl).eq.4) then
k1=max(1 ,k-1)
k2=min(kk,k+1)
varb2d(i1,j1)= &
0.5*( q *(saln(i,j,k1,n)+saln(i,j,k,n))+ &
(1.0-q)*(saln(i,j,k2,n)+saln(i,j,k,n)))
else
varb2d(i1,j1)=saln(i,j,k,n)
endif
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
flt(nfl,1),flt(nfl,2),maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,sflt,ier)
flt(nfl,6)=sflt
!
! --- calculate density
if (ifltyp(nfl).ne.2) then
thflt=sig(tflt,sflt)
flt(nfl,7)=thflt
else
tflt=tofsig(flt(nfl,7),sflt)
flt(nfl,5)=tflt
endif
!
! --- interpolate model fields for float type 4 (stationary)
if (ifltyp(nfl).eq.4) then
!
! --- 3-d u and v interpolation
ngrid=2
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
k1=max(1 ,k-1)
k2=min(kk,k+1)
#if defined(STOKES)
varb2d(i1,j1)=ubavg(i,j,n)+ &
0.5*( q *(u(i,j,k1,n)+u(i,j,k,n)+ &
usd(i,j,k1)+usd(i,j,k))+ &
(1.0-q)*(u(i,j,k2,n)+u(i,j,k,n)+ &
usd(i,j,k2)+usd(i,j,k)))
#else
varb2d(i1,j1)=ubavg(i,j,n)+ &
0.5*( q *(u(i,j,k1,n)+u(i,j,k,n))+ &
(1.0-q)*(u(i,j,k2,n)+u(i,j,k,n)))
#endif
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
flt(nfl,1),flt(nfl,2),maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,uflt,ier)
!
ngrid=3
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
k1=max(1 ,k-1)
k2=min(kk,k+1)
#if defined(STOKES)
varb2d(i1,j1)=vbavg(i,j,n)+ &
0.5*( q *(v(i,j,k1,n)+v(i,j,k,n)+ &
vsd(i,j,k1)+vsd(i,j,k))+ &
(1.0-q)*(v(i,j,k2,n)+v(i,j,k,n)+ &
vsd(i,j,k2)+vsd(i,j,k)))
#else
varb2d(i1,j1)=vbavg(i,j,n)+ &
0.5*( q *(v(i,j,k1,n)+v(i,j,k,n))+ &
(1.0-q)*(v(i,j,k2,n)+v(i,j,k,n)))
#endif
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
flt(nfl,1),flt(nfl,2),maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,vflt,ier)
!
! --- wveli interpolation
ngrid=1
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
varb2d(i1,j1)=q*wveli(i,j,k)+(1.0-q)*wveli(i,j,k+1)
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
flt(nfl,1),flt(nfl,2),maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,wflt1,ier)
flt(nfl,10)=wflt1/deltfl
!
! --- vertical viscosity coefficient interpolation
ngrid=1
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
varb2d(i1,j1)=q*vcty(i,j,k)+(1.0-q)*vcty(i,j,k+1)
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
flt(nfl,1),flt(nfl,2),maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,vkflt,ier)
flt(nfl,11)=max(1.0e-4,vkflt)
!
! --- vertical temperature diffusivity coefficient interpolation
ngrid=1
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
varb2d(i1,j1)=q*dift(i,j,k)+(1.0-q)*dift(i,j,k+1)
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
flt(nfl,1),flt(nfl,2),maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,tkflt,ier)
flt(nfl,12)=max(1.0e-5,tkflt)
!
! --- vertical scalar diffusivity coefficient interpolation
ngrid=1
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
varb2d(i1,j1)=q*difs(i,j,k)+(1.0-q)*difs(i,j,k+1)
enddo
enddo
call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), &
flt(nfl,1),flt(nfl,2),maskpt(1,1,ngrid), &
ngood(ngrid),ngrid,intpfl,radian,skflt,ier)
flt(nfl,13)=max(1.0e-5,skflt)
!
endif ! float type 4
!
if (nfl.eq.nfl_debug) then
write(lp,111) nfl,tflt,sflt,thflt
111 format('new t,s,th, float',i5,2x,3f9.3)
endif !nfl_debug
!
endif ! ioflag.eq.1
!
! ----------------------------------------------------
! --- flag floats that have run aground or left domain
! ----------------------------------------------------
!
! --- jump ahead if the float has not reached its termination time
if (ntermn.gt.-10) go to 20
!
! --- the following code block is entered when the previous code has
! --- detected that the float needs to be terminated
30 flt(nfl,1)=-999.
flt(nfl,2)=-999.
if (ntermn.gt.0) then
write(lp,112) nfl,nstep,nstepfl,ier,ntermn,ngood
112 format('float',i6,' aground, time step',i9, &
' float time step',i9/' ier ',i3,' exit point', &
i3,' ngood 1-4',4i3)
elseif (ntermn.eq.-5) then
write(lp,113) nfl
113 format('float',i6,' exits domain')
elseif (ntermn.eq.-10) then
write(lp,114) nfl
114 format('float',i6,' reaches termination time')
endif
!
20 continue
!
! ----------------------------------------------------
! --- store velocity components for synthetic moorings
! ----------------------------------------------------
!
if (ifltyp(nfl).eq.4) then
fltloc(nfl,1)=uflt
fltloc(nfl,2)=vflt
fltloc(nfl,3)=wflt
else
fltloc(nfl,1)=0.0
fltloc(nfl,2)=0.0
fltloc(nfl,3)=0.0
endif
!
!--- store updated particle location
!
do i=1,13
flt3( i+17*(nfl-1)) = flt( nfl,i)
enddo
do i=1,3
flt3(13+i+17*(nfl-1)) = fltloc(nfl,i)
enddo
flt3(17+17*(nfl-1)) = kfloat(nfl)
!
! --- else du premier if de la boucle
else !float is not on this tile, so disable flt3
do i=1,17
flt3(i+17*(nfl-1)) = hugel
enddo
endif
!
go to 222
!
22 continue
!
! --- float is unchanged
do i=1,17
flt3(i+17*(nfl-1)) = hugel
enddo
!
222 continue
!
! ---------------------
! --- end of float loop
! ---------------------
!
enddo !nfl
!
! --- copy floats back to 1st processor
call xcminr(flt3(1:17*nflt), 1)
if (mnproc.eq.1) then
do nfl=1,nflt
if (flt3(1+17*(nfl-1)).ne.hugel) then
do i=1,13
flt( nfl,i) = flt3( i+17*(nfl-1))
enddo
do i=1,3
fltloc(nfl,i) = flt3(13+i+17*(nfl-1))
enddo
kfloat(nfl) = flt3(17+17*(nfl-1))
endif !updated float
enddo !nfl
endif !1st tile
!
! -----------------------------------
! --- output float array for analysis
! -----------------------------------
!
if (ioflag.eq.1 .and. mnproc.eq.1) then
!
write(lp,*) 'storing float parameters'
open(unit=uoff+801,file=flnmflto,status='unknown', &
form='formatted',position='append')
do nfl=1,nflt
if ( flt(nfl,1).gt.-999. .and. &
timefl-flt(nfl,8).ge.-0.001 ) then
!
if (samplfl) then
if (ifltyp(nfl).ne.4) then
write(uoff+801,901) iflnum(nfl),timefl,kfloat(nfl), &
(flt(nfl,j),j=1,7)
901 format(i6,f12.4,i3,2f10.4,f14.4,f8.2,3f7.3)
else
write(uoff+801,902) iflnum(nfl),timefl,kfloat(nfl), &
fltloc(nfl,1),fltloc(nfl,2), &
fltloc(nfl,3),flt(nfl,10), &
(flt(nfl,j),j=4,7), &
(flt(nfl,j),j=11,13)
902 format(i6,f12.4,i3,2f10.6,2f14.10,f8.2,3f7.3,1p,3e12.4)
endif
else
if (ifltyp(nfl).ne.4) then
write(uoff+801,901) iflnum(nfl),timefl,kfloat(nfl), &
(flt(nfl,j),j=1,3)
else
write(uoff+801,902) iflnum(nfl),timefl,kfloat(nfl), &
fltloc(nfl,1),fltloc(nfl,2), &
fltloc(nfl,3),flt(nfl,10)
endif
endif
endif
enddo
close(uoff+801)
endif
!
! ----------------------------------------------------
! --- store old velocities for next time interpolation
! ----------------------------------------------------
!
10 continue
!
if (hadvfl) then
!
margin = 6
!
do j=1-margin,jj+margin
do i=1-margin,ii+margin
if (SEA_U) then
do k=1,kk
uold2(i,j,k+kold2)=u(i,j,k,n)+ubavg(i,j,n)
enddo !k
endif !iu
if (SEA_V) then
do k=1,kk
vold2(i,j,k+kold2)=v(i,j,k,n)+vbavg(i,j,n)
enddo !k
endif !iv
if (SEA_P) then
do k=1,kk
wold2u(i,j,k+kold2)=wvelup(i,j,k)
wold2d(i,j,k+kold2)=wveldn(i,j,k)
enddo !k
endif !ip
enddo !i
enddo !j
!
endif
!
return
end subroutine floats
subroutine intrph(varb2d,ptlon,ptlat,fllon,fllat,maskpt, &
ngood,ngrid,intpfl,radian,vrbint,ier)
!
! ----------------------------------------------------
! --- 2-d polynomial or nearest neighbor interpolation
! ----------------------------------------------------
!
! --- polynomial interpolation code is adapted from the parameter matrix
! --- objective analysis algorithm of mariano and brown (1994), specifically
! --- the code used to calculate the large scale 2-d trend surface
!
! --- nptmin is presently set to 9
!
implicit none
real, dimension(4,4) :: varb2d,ptlon,ptlat
real fllon,fllat,fllon2,fllat2,vfbint
real wghtx,wghty,vrbint,suma,wghts,xcen,ycen,xo,yo,xsd,ysd,var
real dx,dy,dxflt,dyflt,det,denom,radian
real, dimension(16) :: at,xt,yt
real, dimension(12) :: parm
real, dimension(10,10) :: rtr,rtri
real, dimension(10) :: rta
integer iarea(200)
integer ngood,ngrid,intpfl,ier,iflag,ict,nptmin
integer i,j,nd,nds,np,na
logical maskpt(4,4)
!
! --- nptmin is the minimum number of points required for 2-d interpolation
! --- to be performed - otherwise, nearest neighbor interpolation is used
!
data nptmin /9/
!
! --- na is the row dimension of rtr
data na /10/
!
! --- set the 1-d arrays of field values and horizontal position
! --- find and remove central latitude and longitude of the good points
! --- correct horizontal position for earth's curvature to lowest order
!
ict=0
xcen=0.0
ycen=0.0
do i=1,4
do j=1,4
if (.not.maskpt(i,j) .or. ngood.eq.16) then
ict=ict+1
at(ict)=varb2d(i,j)
xt(ict)=ptlon(i,j)
xcen=xcen+xt(ict)
yt(ict)=ptlat(i,j)
ycen=ycen+yt(ict)
endif
enddo
enddo
!
np=ict
!
if(np.le.0) then
vrbint=0.0
ier=999
return
endif
!
! --- set interpolated value to zero if all input values are zero
iflag=0
do ict=1,np
if(at(ict).ne.0.) iflag=1
enddo
if(iflag.eq.0) then
vrbint=0.0
return
endif
!
xcen=xcen/np
ycen=ycen/np
denom=cos(ycen*radian)
!
do ict=1,np
xt(ict)=(xt(ict)-xcen)*cos(yt(ict)*radian)/denom
yt(ict)=(yt(ict)-ycen)*cos(0.5*(yt(ict)+ycen)*radian)/denom
!
!diag write(6,101) ngrid,ngood,ict,at(ict),xt(ict),yt(ict)
101 format('poly. - ngrid,ngood,ict,at,xt,yt:' &
/10x,3i4,1p,3e13.5)
!
enddo
!
fllon2=(fllon-xcen)*cos(fllat*radian)/denom
fllat2=(fllat-ycen)*cos(0.5*(fllat+ycen)*radian)/denom
!
! --- perform the interpolation
!
ier=0
!
if (intpfl.eq.0 .and. ngood.ge.nptmin) then
!
! --- 2-d, second order polynomial interpolation
!
! --- set parameters for quadratic polynomial and zero the parm array
nds=2
nd=6
do i=1,12
parm(i)=0.0
enddo
!
! --- calculate xo, the mean of xt, and yo, the mean of yt and their
! --- standard deviations from np observations
!
call f_stat(xt,np,xo,var,xsd)
call f_stat(yt,np,yo,var,ysd)
!
do i=1,nd
rta(i)=0.0
do j=1,nd
rtr(i,j)=0.0
enddo
enddo
!
! --- the next part of the code is a little messy, but we save alot of
! --- disk space by doing it this way
!
do 20 ict=1,np
dx=xt(ict)-xo
dy=yt(ict)-yo
rtr(1,2)=rtr(1,2) + dx
rtr(1,3)=rtr(1,3) + dy
rtr(2,2)=rtr(2,2) + dx**2
rtr(2,3)=rtr(2,3) + dx*dy
rtr(3,3)=rtr(3,3) + dy**2
rta(1)=rta(1) + at(ict)
rta(2)=rta(2) + dx*at(ict)
rta(3)=rta(3) + dy*at(ict)
if(nds.eq.1) go to 20
rtr(2,4)=rtr(2,4) + (dx**3)/2.
rtr(3,5)=rtr(3,5) + (dy**3)/2.
rtr(4,4)=rtr(4,4) + (dx**4)/4.
rtr(5,5)=rtr(5,5) + (dy**4)/4.
rtr(3,4)=rtr(3,4) + (dx**2)*dy/2.
rtr(2,5)=rtr(2,5) + (dy**2)*dx/2.
rtr(4,6)=rtr(4,6) + (dx**3)*dy/2.
rtr(5,6)=rtr(5,6) + (dy**3)*dx/2.
rtr(6,6)=rtr(6,6) + (dx**2)*(dy**2)
rta(4)=rta(4) + (dx**2)*at(ict)/2.
rta(5)=rta(5) + (dy**2)*at(ict)/2.
rta(6)=rta(6) + dx*dy*at(ict)
if(nds.ne.3) go to 20
rtr(4,7)=rtr(4,7) + (dx**5)/12.
rtr(4,8)=rtr(4,8) + (dx**2)*(dy**3)/12.
rtr(4,9)=rtr(4,9) + (dx**4)*dy/4.
rtr(4,10)=rtr(4,10) + (dx**3)*(dy**2)/4.
rtr(5,8)=rtr(5,8) + (dy**5)/12.
rtr(5,10)=rtr(5,10) + dx*(dy**4)/4.
rtr(6,7)=rtr(6,7) + dy*(dx**4)/6.
rtr(7,7)=rtr(7,7) + (dx**6)/36.
rtr(7,8)=rtr(7,8) + (dx**3)*(dy**3)/36.
rtr(7,9)=rtr(7,9) + (dx**5)*dy/12.
rtr(7,10)=rtr(7,10) + (dx**4)*(dy**2)/12.
rtr(8,8)=rtr(8,8) + (dy**6)/36.
rtr(8,9)=rtr(8,9) + (dx**2)*(dy**4)/12.
rtr(8,10)=rtr(8,10) + dx*(dy**5)/12.
rta(7)=rta(7) + (dx**3)*at(ict)/6.
rta(8)=rta(8) + (dy**3)*at(ict)/6.
rta(9)=rta(9) + (dx**2)*dy*at(ict)/2.
rta(10)=rta(10) + (dy**2)*dx*at(ict)/2.
20 continue
rtr(1,1)=real(np)
rtr(2,1)=rtr(1,2)
rtr(3,2)=rtr(2,3)
rtr(3,1)=rtr(1,3)
if(nds.eq.1) go to 21
rtr(1,4)=rtr(2,2)/2.
rtr(1,5)=rtr(3,3)/2.
rtr(1,6)=rtr(2,3)
rtr(2,6)=rtr(3,4)*2.
rtr(3,6)=rtr(2,5)*2.
rtr(4,5)=rtr(6,6)/4.
rtr(4,1)=rtr(1,4)
rtr(4,2)=rtr(2,4)
rtr(4,3)=rtr(3,4)
rtr(5,1)=rtr(1,5)
rtr(5,2)=rtr(2,5)
rtr(5,3)=rtr(3,5)
rtr(5,4)=rtr(4,5)
rtr(6,1)=rtr(1,6)
rtr(6,2)=rtr(2,6)
rtr(6,3)=rtr(3,6)
rtr(6,4)=rtr(4,6)
rtr(6,5)=rtr(5,6)
if(nds.ne.3) go to 21
rtr(1,7)=rtr(2,4)/3.
rtr(1,8)=rtr(3,5)/3.
rtr(1,9)=rtr(2,6)/2.
rtr(1,10)=rtr(2,5)
rtr(2,7)=2.*rtr(4,4)/3.
rtr(2,8)=rtr(5,6)/3.
rtr(2,9)=rtr(4,6)
rtr(2,10)=2.*rtr(4,5)
rtr(3,7)=rtr(4,6)/3.
rtr(3,8)=2.*rtr(5,5)/3.
rtr(3,9)=2.*rtr(4,5)
rtr(3,10)=3.*rtr(2,8)
rtr(5,7)=rtr(4,10)/3.
rtr(5,9)=rtr(4,8)/3.
rtr(6,8)=2.*rtr(5,10)/3.
rtr(6,9)=2.*rtr(4,10)
rtr(6,10)=2.*rtr(5,9)
rtr(7,1)=rtr(1,7)
rtr(7,2)=rtr(2,7)
rtr(7,3)=rtr(3,7)
rtr(7,4)=rtr(4,7)
rtr(7,5)=rtr(5,7)
rtr(7,6)=rtr(6,7)
rtr(8,1)=rtr(1,8)
rtr(8,2)=rtr(2,8)
rtr(8,3)=rtr(3,8)
rtr(8,4)=rtr(4,8)
rtr(8,5)=rtr(5,8)
rtr(8,6)=rtr(6,8)
rtr(8,7)=rtr(7,8)
rtr(9,1)=rtr(1,9)
rtr(9,2)=rtr(2,9)
rtr(9,3)=rtr(3,9)
rtr(9,4)=rtr(4,9)
rtr(9,5)=rtr(5,9)
rtr(9,6)=rtr(6,9)
rtr(9,7)=rtr(7,9)
rtr(9,8)=rtr(8,9)
rtr(9,9)=3.*rtr(7,10)
rtr(9,10)=9.*rtr(7,8)
rtr(10,1)=rtr(1,10)
rtr(10,2)=rtr(2,10)
rtr(10,3)=rtr(3,10)
rtr(10,4)=rtr(4,10)
rtr(10,5)=rtr(5,10)
rtr(10,6)=rtr(6,10)
rtr(10,7)=rtr(7,10)
rtr(10,8)=rtr(8,10)
rtr(10,9)=rtr(9,10)
rtr(10,10)=3.*rtr(8,9)
21 continue
!
! --- invert rtr and store in rtri, na is the row dimension of rtr
! --- and rtri, det is the determinant, nd is the number of coeff.
! --- to be estimated, i.e. the number of independent variables in
! --- rtr, iarea is a work area and ier is an error code. if ier
! --- is greater than 0, there is an inversion error
!
call f_invmtx(rtr,na,rtri,na,nd,det,iarea,ier)
if(ier.gt.0) print *, 'WARNING: INVERSION ERROR, ier= ',ier
!
! --- calculate the regression coeff. - parm
!
do i=1,nd
do j=1,nd
parm(i)=parm(i) + rtri(i,j)*rta(j)
enddo
enddo
!
parm(11)=xo
parm(12)=yo
85 continue
!
! --- interpolate the field to the float location
dxflt=fllon2-parm(11)
dyflt=fllat2-parm(12)
vrbint=parm(1)+parm(2)*dxflt+parm(3)*dyflt &
+parm(4)*(dxflt**2)/2.+parm(5)*(dyflt**2)/2.+ &
parm(6)*dxflt*dyflt
return
!
else
!
! --- 2-d nearest neighbor interpolation if the number of available water
! --- points is smaller than nptmin. variable xt is inverse distance
!
wghts=0.
do ict=1,np
xt(ict)=1.0/sqrt(max(1.0e-20,(xt(ict)-fllon2)**2+ &
(yt(ict)-fllat2)**2))
wghts=wghts+xt(ict)
!diag write(6,102) ngrid,ngood,i,j,ict,at(ict),xt(ict)
102 format('nearest neighbor - ngrid,ngood,i,j,ict,at,xt:' &
/10x,5i4,1p,2e13.5)
enddo
!
suma=0.
do ict=1,np
suma=suma+at(ict)*xt(ict)
enddo
vrbint=suma/wghts
return
endif
!
end subroutine intrph
!
!
SUBROUTINE F_INVMTX (A,NA,V,NV,N,D,IP,IER)
implicit none
!
! --- this matrix inversion code is used in the parameter matrix
! --- objective analysis algorithm of mariano and brown (1994), specifically
! --- in the code used to calculate the large scale 2-d trend surface
!
INTEGER NA,NV,N,IP(2*N),IER
REAL D
REAL, DIMENSION(NA,N) :: A
REAL, DIMENSION(NV,N) :: V
! DATA IEXMAX/75/
integer, parameter :: iexmax=75
integer j,m,i,iex, k, l, npm
real vmax,ch, pvt,pvtmx, hold, vh
115 FORMAT(28H0*MATRIX SINGULAR IN INVMTX*)
116 FORMAT(34H0*DETERMINANT TOO LARGE IN INVMTX*)
IER = IERINV_F(N,NA,NV)
IF (IER .NE. 0) RETURN
DO 102 J=1,N
IP(J) = 0
DO 101 I=1,N
V(I,J) = A(I,J)
101 CONTINUE
102 CONTINUE
D = 1.
IEX = 0
DO 110 M=1,N
VMAX = 0.
DO 104 J=1,N
IF (IP(J) .NE. 0) GO TO 104
DO 103 I=1,N
IF (IP(I) .NE. 0) GO TO 103
VH =ABS(V(I,J))
IF (VMAX .GT. VH) GO TO 103
VMAX = VH
K = I
L = J
103 CONTINUE
104 CONTINUE
IP(L) = K
NPM = N+M
IP(NPM) = L
D = D*V(K,L)
105 IF (ABS(D) .LE. 1.0) GO TO 106
D = D*0.1
IEX = IEX+1
GO TO 105
106 CONTINUE
PVT = V(K,L)
IF (M .EQ. 1) PVTMX = ABS(PVT)
IF (ABS(PVT/REAL(M))+PVTMX .EQ. PVTMX) GO TO 113
V(K,L) = 1.
DO 107 J=1,N
HOLD = V(K,J)
V(K,J) = V(L,J)
V(L,J) = HOLD/PVT
107 CONTINUE
DO 109 I=1,N
IF (I .EQ. L) GO TO 109
HOLD = V(I,L)
V(I,L) = 0.
DO 108 J=1,N
V(I,J) = V(I,J)-V(L,J)*HOLD
108 CONTINUE
109 CONTINUE
110 CONTINUE
M = N+N+1
DO 112 J=1,N
M = M-1
L = IP(M)
K = IP(L)
IF (K .EQ. L) GO TO 112
D = -D
DO 111 I=1,N
HOLD = V(I,L)
V(I,L) = V(I,K)
V(I,K) = HOLD
111 CONTINUE
112 CONTINUE
IF (IEX .GT. IEXMAX) GO TO 114
D = D*10.**IEX
RETURN
113 IER = 33
! PRINT 115
RETURN
114 IER = 1
D = REAL(IEX)
PRINT 116
RETURN
END SUBROUTINE F_INVMTX
!
INTEGER FUNCTION IERINV_F (N,NA,NV)
implicit none
INTEGER N, NA, NV
103 FORMAT(23H0* N .LT. 1 IN INVMTX *)
104 FORMAT(24H0* NA .LT. N IN INVMTX *)
105 FORMAT(24H0* NV .LT. N IN INVMTX *)
IERINV_F = 0
IF (N .GE. 1) GO TO 101
IERINV_F = 34
PRINT 103
RETURN
101 IF (NA .GE. N) GO TO 102
IERINV_F = 35
PRINT 104
RETURN
102 IF (NV .GE. N) RETURN
IERINV_F = 36
PRINT 105
RETURN
END FUNCTION IERINV_F
subroutine f_stat(ser,ls,amean,var,std)
implicit none
!
! --- computes mean, variance, standard deviation of data sequence
real, dimension(16) :: ser
integer ls
real amean,var,std
real sum,value
integer j
sum=0.0
do j=1,ls
sum=sum+ser(j)
enddo
amean=sum/ls
sum=0.0
do j=1,ls
value=ser(j)-amean
sum=sum+value*value
enddo
var=sum/ls
std=sqrt(var)
!
return
end subroutine f_stat
!
!
! --- following are subprograms provided by Geoff Samuels, and inserted
! --- by Nasseer Idrisi (10/5/01)
subroutine rantab_ini(rtab,iseed1,iseed2,iflag)
implicit none
!
real*4 rtab(200,2)
integer*4 iseed1,iseed2
integer*4 iflag
integer ic, inisee, itime
!
! --- get the two seeds
! --- iseed1 is for the 'usual' random number generator (udev_rt)
! --- iseed2 is for the 'gaussian' random number generator (grand_rt)
!
integer*4 seed
integer*2 isd2(2)
#if defined(NAGFOR)
!
! --- NAG Fortran - this routine won't work without the equivalence
!*****equivalence (isd2,seed)
!
call system_clock(inisee)
!
if (inisee.ne.-999) then !always .true.
call xcstop('(rantab_ini-NAG Fortran)')
stop '(rantab_ini-NAG Fortran)'
endif !.true.
#else
equivalence (isd2,seed)
!
call system_clock(inisee)
#endif
!
!diag write(*,*)'calc itime',time()
itime = inisee
seed = 123456789
iseed1 = isd2(1)
if (iseed1 .lt. 0) then
iseed1 = abs(iseed1)
endif
iseed1 = iseed1 + itime*iflag
!diag write(*,*)'iseed1',iseed1
seed = 987654321
iseed2 = isd2(2)
if (iseed2 .lt. 0) then
iseed2 = abs(iseed1)
endif
iseed2 = iseed2 + itime*iflag
!diag write(*,*)'iseed2',iseed2
!
! --- set seeds to negative values if iflag=1
!
if (iflag .eq. 1) then
iseed1 = -iseed1
iseed2 = -iseed2
endif
!
! --- load the table with a sequence of 200 random numbers.
do ic = 1,200
rtab(ic,1) = udev_rt(iseed1)
rtab(ic,2) = grand_rt(iseed2)
enddo
!diag write(*,*)'rtabs',rtab(100,1),rtab(100,2)
!
return
end subroutine rantab_ini
subroutine rantab (rtab,iseed1,iseed2,numran)
!
! --- each time this routine is called one of the 200 records in "rtab"
! --- is randomly selected and new random numbers are inserted
!
real*4 rtab(200,2)
integer*4 iseed1,iseed2,numran
!
! --- pick a number between 1 and 200
numran = nint(1+199*udev_rt(iseed1))
if (numran .gt. 200) then
numran = 200
endif
!
! --- get new random numbers
rtab(numran,1) = udev_rt(iseed1)
rtab(numran,2) = grand_rt(iseed2)
!
return
end subroutine rantab
!
real function grand_rt(idum)
implicit none
!
! provides random number from gaussian distribution, mean 0, varaince 1,
! for finding random component of turbulent u
!
! from arthur mariano using routines from numerical recipes
!
integer*4 idum
integer*4 iset
real*4 gset,v1,v2,r,fac
! data iset/0/
! save iset,gset
!
iset = 0
gset=0.0
if (iset.eq.0) then
1 v1=2.*udev_rt(idum)-1.0
v2=2.*udev_rt(idum)-1.0
r=v1**2+v2**2
if (r.ge.1.0)go to 1
fac=sqrt(-2.*log(r)/r)
gset=v1*fac
grand_rt=v2*fac
iset=1
else
grand_rt=gset
iset=0
endif
!
return
end function grand_rt
!
real function udev_rt(iseed)
integer*4 iseed
integer, parameter :: m=714025,ia=1366,ic=150889
real rm
integer, dimension(97) :: ir
integer iy, iff,j
!
iff = 0
rm = 1./m
if (iseed .lt. 0 .or. iff .eq. 0) then
iff = 1
iseed = mod(ic-iseed,m)
do j=1,97
iseed = mod(ia*iseed+ic,m)
ir(j) = iseed
enddo
iseed = mod(ia*iseed+ic,m)
else
iseed = mod(iseed,m)
endif
iy=iseed
j = min(97,1+(97*iy)/m)
iy = ir(j)
udev_rt = iy*rm
iseed = mod(ia*iseed+ic,m)
ir(j) = iseed
return
end function udev_rt
end module mod_floats
!
!
!> Revision history:
!>
!> Nov 2006 - 1st module version
!> Nov 2012 - implicit none added by Till Andreas Rasmussen
!> May 2014 - use land/sea masks (e.g. ip) to skip land
!> Aug 2017 - Added macro to allow NAG Fortran to compile