subroutine derivs0 (time,pos,vel,tflag) c c======================================================================= c === c This routine computes the velocity needed to compute a Lagrangian === c trajectory via a fourth order Runge-Kutta scheme. === c === c ------ === c Input: === c ------ === c === c POS the position vector at the last time (grid,grid,meters; === c real array) === c TIME new time to compute particle position (real) === c TFLAG indicates where in time to compute velocities (integer) === c TFLAG = 0 -> compute at last time === c TFLAG = 1 -> compute midway in time interval === c TFLAG = other -> compute at new time === c === c Common blocks: === c === c /FIELDS/ === c === c HV reciprocal of depth at velocity points (real array)=== c P transport streamfunction at advanced time (real array)=== c PB transport streamfunction at current time (real array)=== c === c /FULLWD/ === c === c NDISK index for unit number, current data. (integer) === c NDISKB index for unit number, previous data. (integer) === c LABS rolling unit numbers for virtual data. (integer array)=== c === c /ONEDIM/ === c === c CSR reciprocal of cosine metric factor. (real vector) === c DXU2R reciprocal of twice UV x-grid spacing. (real vector) === c DYU2R reciprocal of twice UV y-grid spacing. (real vector) === c === c /WORKSPB/ === c === c PTD change in transport streamfunction. (real array) === c === c ------- === c Output: === c ------- === c === c VEL the velocities needed (real array) === c === c Common blocks: === c === c /TRKSCL/ === c === c DTYP Drifter type. (character) === c FACTT time scale. (real; sec) === c FACTXY horizontal space scale. (real; cm) === c HD Diameter (height) of extended drifter. (real; cm) === c === #ifndef ext_tide c Calls: BES1D, OPICK, VAVGV, VNTRPV === # else c Calls: BES1D, DRFTDV, OPICK, VAVGV, VNTRPV === #endif c === c Note: The variables U, UB, V and VB in the common block /WORKSP/ === c are used as local workspace. === c === c======================================================================= c c----------------------------------------------------------------------- c Define global data. c----------------------------------------------------------------------- c #include #include #include #include #include #include #include #include #include c c----------------------------------------------------------------------- c Define local data. c----------------------------------------------------------------------- c logical first integer i,j,ix,jy,ibry,jbry,jwk,iwk,neu,nev,nsu,nsv, * tflag,iwkp1,jwkp1 #ifdef cyclic * ,ncyc #endif FLOAT * bes1d,diag1,diag2,diag1b,diag2b,dtc,dto,time,usclr, * uz,uzb,vsclr,vz,vzb,x,xmix,y,ymax,ymjy,z #ifndef cyclic * ,xmax #endif FLOAT * pos(3),vel(3),uwk1(4),uwk2(4),vwk1(4),vwk2(4) parameter (nsu=1+imtkm*nt,nsv=nsu+imtkm,neu=nsu-1+imtkm, * nev=nsv-1+imtkm) #ifndef cyclic save first,usclr,vsclr,xmax,ymax # else save first,usclr,vsclr,ymax #endif data first/.true./ c c----------------------------------------------------------------------- c Begin executable code. c----------------------------------------------------------------------- c c----------------------------------------------------------------------- c On first pass, set bounds for horizontal interpolations. c----------------------------------------------------------------------- c if(first) then #ifndef cyclic xmax=FLoaT(imtm2) #endif ymax=FLoaT(jmtm2) vsclr=factt/facty usclr=factt/factx first=.false. endif c c----------------------------------------------------------------------- c Set up parameters for linear interpolation in time. c----------------------------------------------------------------------- c if(tflag.eq.0) then dto=c1 dtc=c0 elseif(tflag.eq.1) then dto=p5 dtc=p5 else dto=c0 dtc=c1 endif c c----------------------------------------------------------------------- c Set up parameters for spatial interpolations. c----------------------------------------------------------------------- c x=pos(1)-p5 y=pos(2)-p5 z=max(pos(3)*m2cm,c0) c ix=int(x) xmix=x-FLoaT(ix) ibry=0 #ifndef cyclic if(x.lt.c2) then ibry=1 elseif(x.ge.xmax) then ibry=2 endif #endif c jy=int(y) ymjy=y-FLoaT(jy) jbry=0 if(y.lt.c2) then jbry=1 elseif(y.ge.ymax) then jbry=2 endif c c----------------------------------------------------------------------- c Linearly interpolate horizontal velocities in depth & time. Cubicly c interpolate the velocities in zonal direction. c----------------------------------------------------------------------- c do 20 j=1,4 jwk=max(min(jy+j-2,jmtm1),1) jwkp1=jwk+1 call opick(labs(ndisk),nslab,(jwk-1)*nslab+1,nsu,neu,u) call opick(labs(ndisk),nslab,(jwk-1)*nslab+1,nsv,nev,v) call opick(labs(ndiskb),nslab,(jwk-1)*nslab+1,nsu,neu,ub) call opick(labs(ndiskb),nslab,(jwk-1)*nslab+1,nsv,nev,vb) #ifdef ext_tide call drftdv (jwk,u,ub,v,vb) #endif c do 10 i=1,4 c #ifndef cyclic iwk=max(min(ix+i-2,imtm1),1) iwkp1=iwk+1 # else iwk=ix+i-2 ncyc=(iwk-1)/imtm2 if(iwk.lt.1) ncyc=ncyc-1 iwk=iwk-ncyc*imtm2 iwkp1=iwk+1 ncyc=iwk/imtm2 if(iwkp1.lt.1) ncyc=ncyc-1 iwkp1=iwkp1-ncyc*imtm2 #endif c if (dtyp.ne.'x') then call vntrpv (iwk,jwk,z,u,ub,v,vb,uz,uzb,vz,vzb) else call vavgv (iwk,jwk,z,hd,u,ub,v,vb,uz,uzb,vz,vzb) end if c diag1b=p(iwkp1,jwkp1)-p(iwk,jwk)+ptd(iwk,jwk)-ptd(iwkp1,jwkp1) diag2b=p(iwk,jwkp1)-p(iwkp1,jwk)+ptd(iwkp1,jwk)-ptd(iwk,jwkp1) diag1 =pb(iwkp1,jwkp1)-pb(iwk ,jwk) diag2 =pb(iwk ,jwkp1)-pb(iwkp1,jwk) c uwk1(i)=(dto*uzb+dtc*uz) * -((diag1+diag2)*dtc+(diag1b+diag2b)*dto)* * dyu2r(jwk)*hv(iwk,jwk) vwk1(i)=(dto*vzb+dtc*vz) * +((diag1-diag2)*dtc+(diag1b-diag2b)*dto)* * dxu2r(iwk)*csr(jwk)*hv(iwk,jwk) 10 continue c uwk2(j)=bes1d(xmix,uwk1(1),uwk1(2),uwk1(3),uwk1(4),ibry)* * csr(jwk) vwk2(j)=bes1d(xmix,vwk1(1),vwk1(2),vwk1(3),vwk1(4),ibry) 20 continue c c----------------------------------------------------------------------- c Interpolate horizontal velocities to desired meridional position c set vertical velocity to zero. c----------------------------------------------------------------------- c vel(1)=bes1d(ymjy,uwk2(1),uwk2(2),uwk2(3),uwk2(4),jbry)*usclr vel(2)=bes1d(ymjy,vwk2(1),vwk2(2),vwk2(3),vwk2(4),jbry)*vsclr vel(3)=c0 c return end