subroutine xtr_vel(lon,lat,dpth,map_v,imax,kmax,xcoord,xrlngd, * xrlatd,xthetad) c c======================================================================= c === c This routine extracts internal velocities from the model grid === c and interpolates to the sub-domain grid (lon,lat,dpth). After === c interpolation the velocity vectors are rotated to conform to === c the sub-domain grid. === c === c On Input: === c === c LON sub-domain V-point longitude (degrees west, real array)=== c LAT sub-domain V-point latitude (degrees north, real array)=== c DPTH sub-domain depths (centimeter, real array) === c KMAX number of points in the z-direction to interpolate === c (integer) === c IMAX number of points in the x-direction to interpolate === c (integer) === c XCOORD coordinate type of the sub-domain to extract (integer) === c XRLNGD centroid longitude of sub-domain to extract (degrees, === c real) === c XRLATD centroid latitude of sub-domain to extract (degrees, === c real) === c XTHETAD rotation angle of the sub-domain to extract (degrees, === c real) === c === c On Output: === c === c MAP_V sub-domain internal mode velocity (real array) === c === c Calls: DEPTHSLAB, LL2XY, OPICK, ROTANGLE === c === c======================================================================= c c----------------------------------------------------------------------- c Define global data. c----------------------------------------------------------------------- c #include #include #include #include #include #include c c----------------------------------------------------------------------- c Define local data. c----------------------------------------------------------------------- c integer i,ic,imax,ip,j,jold,k,kkc,kmax,l,nsu,nsv,neu,nev,xcoord #ifndef barotropic * ,kc #endif integer icell(4) FLOAT * angle,cos_rot,sin_rot,u_rot,v_rot,x,xdis,xrlatd,xrlngd, * xthetad,y,ydis,z FLOAT * rotangle FLOAT * depth(imt,km,2),dpth(ximtkm),lat(ximt),lon(ximt), * map_v(ximtkm,2),ui(imt,km,2),ucell(4),uzcell(km,4), * vi(imt,km,2),vcell(4),vzcell(km,4),zcell(km,4) c c======================================================================= c Begin executable code. c======================================================================= c c----------------------------------------------------------------------- c Extract and linearly interpolate internal mode velocity to c sub-domain grid. c----------------------------------------------------------------------- c jold=0 do 60 l=1,imax c c Set rotation coefficients. (First call MUST be with non-extracted c system to correctly position stored variables in rotangle). c angle=rotangle(lon(l),lat(l),coord,rlngd,rlatd,thetad) angle=rotangle(lon(l),lat(l),xcoord,xrlngd,xrlatd,xthetad) * -angle sin_rot=sin(angle) cos_rot=cos(angle) c c Find indices for horizontal grid cell. c call ll2xy (lon(l),lat(l),coord,imt,jmt,gridx,gridy,rlngd,rlatd, * delx,dely,thetad,x,y) i=int(x) j=int(y) xdis=x-FLoaT(i) ydis=y-FLoaT(j) c icell(1)=i icell(2)=i+1 icell(3)=i+1 icell(4)=i c c Get slab of data containing grid cell. c if(j.ne.jold) then nsu=1+imtkm*nt nsv=nsu+imtkm neu=nsv-1 nev=neu+imtkm call opick(labs(ndisk),nslab,(j-1)*nslab+1,nsu,neu,ui(1,1,1)) call opick(labs(ndisk),nslab,(j )*nslab+1,nsu,neu,ui(1,1,2)) call opick(labs(ndisk),nslab,(j-1)*nslab+1,nsv,nev,vi(1,1,1)) call opick(labs(ndisk),nslab,(j )*nslab+1,nsv,nev,vi(1,1,2)) c c Get depths at slabs u-points. c call depthslab(j ,vgrid,depth(1,1,1)) call depthslab(j+1,vgrid,depth(1,1,2)) j=jold endif c c Pass information from slabs to a vertical column. c do 10 k=1,km zcell(k,1)=depth(icell(1),k,1) zcell(k,2)=depth(icell(2),k,1) zcell(k,3)=depth(icell(3),k,2) zcell(k,4)=depth(icell(4),k,2) uzcell(k,1)=ui(icell(1),k,1) uzcell(k,2)=ui(icell(2),k,1) uzcell(k,3)=ui(icell(3),k,2) uzcell(k,4)=ui(icell(4),k,2) vzcell(k,1)=vi(icell(1),k,1) vzcell(k,2)=vi(icell(2),k,1) vzcell(k,3)=vi(icell(3),k,2) vzcell(k,4)=vi(icell(4),k,2) 10 continue c c Interpolation. c do 50 k=1,kmax c c Vertical interpolation. c ip=k+(l-1)*kmax z=abs(dpth(ip))*m2cm do 40 ic=1,4 if(z.le.zcell(1,ic)) then ucell(ic)=uzcell(1,ic) vcell(ic)=vzcell(1,ic) elseif(z.ge.zcell(km,ic)) then ucell(ic)=uzcell(km,ic) vcell(ic)=vzcell(km,ic) else #ifndef barotropic do 20 kc=1,km-1 if((zcell(kc,ic).lt.z).and. * (z.le.zcell(kc+1,ic))) then kkc=kc goto 30 endif 20 continue 30 continue #endif ucell(ic)=((z-zcell(kkc,ic))*uzcell(kkc+1,ic)+ * (zcell(kkc+1,ic)-z)*uzcell(kkc,ic))/ * (zcell(kkc+1,ic)-zcell(kkc,ic)) vcell(ic)=((z-zcell(kkc,ic))*vzcell(kkc+1,ic)+ * (zcell(kkc+1,ic)-z)*vzcell(kkc,ic))/ * (zcell(kkc+1,ic)-zcell(kkc,ic)) endif 40 continue c c Horizontal interpolation. c map_v(ip,1)=(c1-ydis)*((c1-xdis)*ucell(1)+xdis*ucell(2))+ * ydis*((c1-xdis)*ucell(4)+xdis*ucell(3)) map_v(ip,2)=(c1-ydis)*((c1-xdis)*vcell(1)+xdis*vcell(2))+ * ydis*((c1-xdis)*vcell(4)+xdis*vcell(3)) c c Rotate velocity vector to conform to sub-domain orientation. c u_rot=map_v(ip,1)*cos_rot-map_v(ip,2)*sin_rot v_rot=map_v(ip,1)*sin_rot+map_v(ip,2)*cos_rot map_v(ip,1)=u_rot map_v(ip,2)=v_rot 50 continue 60 continue return end