! $Id: init_floats.F 1458 2014-02-03 15:01:25Z gcambon $ ! !====================================================================== ! CROCO is a branch of ROMS developped at IRD and INRIA, in France ! The two other branches from UCLA (Shchepetkin et al) ! and Rutgers University (Arango et al) are under MIT/X style license. ! CROCO specific routines (nesting) are under CeCILL-C license. ! ! CROCO website : http://www.croco-ocean.org !====================================================================== ! #include "cppdefs.h" subroutine init_floats #ifdef FLOATS ! !================================================== John M. Klinck === ! Copyright (c) 2000 Rutgers/UCLA ! !================================================ Hernan G. Arango === ! ! ! This routine reads in and process initial float locations from ! ! input floats script. ! ! ! !===================================================================== ! implicit none # include "param.h" # include "grid.h" # include "ocean3d.h" # include "ncscrum_floats.h" # include "scalars.h" # include "floats.h" # include "init_floats.h" ! integer Ncount, i, icard, iunit, j, k, lstr, mc, nc integer index(Mfloats) integer lenstr real xfloat, yfloat, zfloat real Ip0(Mfloats), Jp0(Mfloats), lat(Mfloats), lon(Mfloats) character*35 frmt parameter (iunit=50) ! k=0 ! !--------------------------------------------------------------------- ! Read in initial float locations. !--------------------------------------------------------------------- ! lstr=lenstr(fposnam) open(iunit,file=fposnam(1:lstr),form='formatted', & status='old', err=195) c write(stdout,10) c 10 format(/,' FLOATS processing parameters:',/) ! ! Read input parameters according to their input card number. ! icard=0 do while (icard.lt.99) ! ! Read in floats identification title. ! if (icard.eq.1) then read(iunit,'(a)',err=70) Ftitle lstr=lenstr(Ftitle) write(stdout,20) Ftitle(1:lstr) 20 format(1x,'(',a,')',/) ! ! Read in initial floats location. ! elseif (icard.eq.2) then i=0 nfloats=0 c write(stdout,30) c30 format(1x,'Fcoor',2x,'Ftype',2x,'Fcount',2x,'Ft0', c & 2x,'Fx0',2x,'Fy0',2x,'Fz0',2x,'Fgrd',2x,'Fdt', c & 2x,'Fdx',2x,'Fdy',2x,'Fdz',/) do while (.true.) i=i+1 read(iunit,*,err=40) Ft0(i), Fx0(i), Fy0(i), Fz0(i), & Fgrd(i),Fcoor(i), Ftype(i), Fcount(i), & Fdt(i), Fdx(i), Fdy(i), Fdz(i) nfloats=nfloats+Fcount(i) c# if !defined SPHERICAL c if (Fcoor(i).eq.1) then c frmt='(i1,i2,i5,f10.4,2f8.0,f8.2,2x,i1,4f9.3)' c else c frmt='(i1,i2,i5,f10.4,3f8.2,2x,i1,4f9.3)' c endif c# else c frmt='(i1,i2,i5,f10.4,3f8.2,2x,i1,4f9.3)' c# endif c write(stdout,frmt) Fcoor(i), Ftype(i), Fcount(i), c & Ft0(i), Fx0(i), Fy0(i), Fz0(i), c & Fgrd(i),Fdt(i), Fdx(i), Fdy(i), Fdz(i) enddo 40 if (Ft0(i).ne.99.0) then write(stdout,50) icard, i, fposnam 50 format(/,' INIT_FLOATS - error while reading input card: ', & i2, ', floats location entry: ',i3,/,15x, & 'in input script: ',a) stop else Ncount=i-1 i_floats=i-1 goto 90 endif if (Mfloats.lt.nfloats) then write(stdout,60) Mfloats, nfloats 60 format(/,' INIT_FLOATS - too small dimension parameter,', & ' Mfloats',2i6,/,15x,'change file param.h and', & ' recompile.') stop endif endif ! ! Read last input card ID. ! read(iunit,*,err=70) icard enddo goto 90 ! ! Error while reading input parameters. ! 70 write(stdout,80) icard, fposnam 80 format(/,' INIT_FLOATS - error while reading input card: ', & i2,15x,'in input script: ',a) 90 close(iunit) write(stdout,100) nfloats 100 format(/,2x,i6,4x,'nfloats',t26, & 'Number of float trajectories to compute.',/) ! !--------------------------------------------------------------------- ! Set initial float location. !--------------------------------------------------------------------- ! ! Set time of float release (seconds after model initialization) and ! initial float horizontal positions (grid units). Fill the initial ! vertical level or depth position. ! mc=0 nc=0 do i=1,Ncount if (Fcount(i).eq.1) then nc=nc+1 Tinfo(itstr,nc)=(tdays+Ft0(i))*day2sec Tinfo(izgrd,nc)=Fz0(i) Tinfo(igrd,i)=FLOAT(Fgrd(i)) ! each float belongs to this grid initially if (Fcoor(i).eq.0) then Tinfo(ixgrd,nc)=MIN(MAX(0.5,Fx0(i)),FLOAT(LLm)+0.5) Tinfo(iygrd,nc)=MIN(MAX(0.5,Fy0(i)),FLOAT(MMm)+0.5) else mc=mc+1 lon(mc)=Fx0(i) lat(mc)=Fy0(i) index(mc)=nc endif elseif (Fcount(i).gt.1) then do j=1,Fcount(i) nc=nc+1 Tinfo(igrd,nc)=FLOAT(Fgrd(i)) ! each float belongs to this grid initially if (Fdt(i).gt.0.0) then Tinfo(itstr,nc)=(tdays+Ft0(i)+FLOAT(j-1)*Fdt(i))* & day2sec Tinfo(izgrd,nc)=Fz0(i) if (Fcoor(i).eq.0) then Tinfo(ixgrd,nc)=MIN(MAX(0.5,Fx0(i)),FLOAT(LLm)+0.5) Tinfo(iygrd,nc)=MIN(MAX(0.5,Fy0(i)),FLOAT(MMm)+0.5) else mc=mc+1 lon(mc)=Fx0(i) lat(mc)=Fy0(i) index(mc)=nc endif else Tinfo(itstr,nc)=(tdays+Ft0(i))*day2sec if (Fdz(i).eq.0.0) then Tinfo(izgrd,nc)=Fz0(i) else if (Fz0(i).gt.0.0) then zfloat=Fz0(i)+FLOAT(j-1)*Fdz(i) Tinfo(izgrd,nc)=MIN(MAX(0.0,zfloat),FLOAT(N)) else Tinfo(izgrd,nc)=Fz0(i)+FLOAT(j-1)*Fdz(i) endif endif if (Fcoor(i).eq.0) then xfloat=Fx0(i)+FLOAT(j-1)*Fdx(i) Tinfo(ixgrd,nc)=MIN(MAX(0.5,xfloat), & FLOAT(LLm)+0.5) yfloat=Fy0(i)+FLOAT(j-1)*Fdy(i) Tinfo(iygrd,nc)=MIN(MAX(0.5,yfloat), & FLOAT(MMm)+0.5) else mc=mc+1 index(mc)=nc lon(mc)=Fx0(i)+FLOAT(j-1)*Fdx(i) lat(mc)=Fy0(i)+FLOAT(j-1)*Fdy(i) endif endif enddo endif enddo ! ! Set number of floats trajectories to compute. ! nfloats=nc ! ! If applicable, convert floats initial (lon,lat) positions to grid ! units. ! if (mc.gt.0) then # ifdef SPHERICAL call hindices (Ip0,Jp0,lon,lat,mc,lonr,latr) # else call hindices (Ip0,Jp0,lon,lat,mc,xr,yr) # endif do i=1,mc nc=index(i) Tinfo(ixgrd,nc)=MIN(MAX(0.5,Ip0(i)),FLOAT(LLm)+0.5) Tinfo(iygrd,nc)=MIN(MAX(0.5,Jp0(i)),FLOAT(MMm)+0.5) enddo endif ! return 195 write(stdout,205) fposnam 205 format(/,'float file ',A,/,' not found => float initialization', & ' using restart file') call get_initial_floats ldefflt= .false. C! Test de l'initialisation des flotteurs xa C do nc=1,nfloats C write(stdout,53) nc, Tinfo(ixgrd,nc), Tinfo(iygrd,nc), C & Tinfo(izgrd,nc), NINT(Tinfo(igrd,nc)), Tinfo(itstr,nc) C enddo C! Fin test #endif /* FLOATS */ return end #ifdef FLOATS subroutine hindices (Ipos,Jpos,Xpos,Ypos,Npos,Xgrd,Ygrd) ! !================================================ Hernan G. Arango === ! Copyright (c) 2000 Rutgers/UCLA ! !======================================== Alexander F. Shchepetkin === ! ! ! Given position vectors Xpos and Ypos of size Npos, this routine ! ! finds the corresponding indices Ipos and Jpos of the model grid ! ! (Xgrd,Ygrd) cell containing each requested position. ! ! ! ! Calls: Try_Range ! ! ! !===================================================================== ! #ifdef AGRIF USE Agrif_Util #endif implicit none # include "param.h" ! logical found logical Try_Range integer Imax, Imin, Jmax, Jmin, Npos, i0, j0, k real Ipos(Npos), Jpos(Npos), Xpos(Npos), Ypos(Npos), & deltax, deltay, dx1, dy1, dx2, dy2, c1, c2 real Xgrd(GLOBAL_2D_ARRAY), & Ygrd(GLOBAL_2D_ARRAY) ! !--------------------------------------------------------------------- ! Determine grid cell indices containing requested position points. ! Then, interpolate to fractional cell position. !--------------------------------------------------------------------- ! ! Initialize all indices. ! do k=1,Npos Ipos(k)=0.0 Jpos(k)=0.0 enddo ! ! Check each position to find if it falls inside the whole domain. ! Once it is stablished that it inside, find the exact cell to which ! it belongs by successively dividing the domain by a half (binary ! search). ! do k=1,Npos found=Try_Range(0,LLm+1,0,MMm+1,Xpos(k),Ypos(k),Xgrd,Ygrd) if (found) then Imin=0 Imax=LLm+1 Jmin=0 Jmax=MMm+1 do while (((Imax-Imin).gt.1).or.((Jmax-Jmin).gt.1)) if ((Imax-Imin).gt.1) then i0=(Imin+Imax)/2 found=Try_Range(Imin,i0,Jmin,Jmax,Xpos(k),Ypos(k), & Xgrd,Ygrd) if (found) then Imax=i0 else Imin=i0 endif endif if ((Jmax-Jmin).gt.1) then j0=(Jmin+Jmax)/2 found=Try_Range(Imin,Imax,Jmin,j0,Xpos(k),Ypos(k), & Xgrd,Ygrd) if (found) then Jmax=j0 else Jmin=j0 endif endif enddo ! Improved interpolation block (02/06/03). Not totally ! accurate in case where the grid is strongly deformed. ! The nonalignment of xsi (eta) axis with longitude ! (latitude) direction is taken into account. dy1=Ygrd(Imin,Jmin+1)-Ygrd(Imin,Jmin) dx1=Xgrd(Imin,Jmin+1)-Xgrd(Imin,Jmin) dy2=Ygrd(Imin+1,Jmin)-Ygrd(Imin,Jmin) dx2=Xgrd(Imin+1,Jmin)-Xgrd(Imin,Jmin) c1=Xpos(k) *dy1-Ypos(k) *dx1 c2=Xgrd(Imin,Jmin)*dy2-Ygrd(Imin,Jmin)*dx2 deltax=(c1*dx2-c2*dx1)/(dx2*dy1-dy2*dx1) deltax=(deltax-Xgrd(Imin,Jmin))/dx2 Ipos(k)=FLOAT(Imin)+MIN(MAX(0.0,deltax),1.0) c1=Xgrd(Imin,Jmin)*dy1-Ygrd(Imin,Jmin)*dx1 c2=Xpos(k) *dy2-Ypos(k) *dx2 deltay=(c1*dy2-c2*dy1)/(dx2*dy1-dy2*dx1) deltay=(deltay-Ygrd(Imin,Jmin))/dy1 Jpos(k)=FLOAT(Jmin)+MIN(MAX(0.0,deltay),1.0) endif enddo return end function Try_Range (Imin,Imax,Jmin,Jmax,Xo,Yo,Xgrd,Ygrd) ! !================================================ Hernan G. Arango === ! Copyright (c) 2000 Rutgers/UCLA ! !======================================== Alexander F. Shchepetkin === ! ! ! Given a grided domain with matrix coordinates Xgrd and Ygrd, this ! ! function finds if the point (Xo,Yo) is inside the box defined by ! ! the requested corners (Imin,Jmin) and (Imax,Jmax). It will return ! ! logical switch Try_Range=.true. if (Xo,Yo) is inside, otherwise ! ! it will return false. ! ! ! ! Calls: inside ! ! ! !===================================================================== ! implicit none #include "param.h" ! logical Try_Range, inside integer Imax, Imin, Jmax, Jmin, Nb, NX, i, j, shft parameter(NX=2*(LLm0+2)+2*(MMm0+2)+1) real Xb(NX), Yb(NX), Xo, Yo real Xgrd(GLOBAL_2D_ARRAY), & Ygrd(GLOBAL_2D_ARRAY) ! !--------------------------------------------------------------------- ! Define closed polygon. !--------------------------------------------------------------------- ! ! Note that the last point (Xb(Nb),Yb(Nb)) does not repeat first ! point (Xb(1),Yb(1)). Instead, in function inside, it is implied ! that the closing segment is (Xb(Nb),Yb(Nb))-->(Xb(1),Yb(1)). In ! fact, function inside sets Xb(Nb+1)=Xb(1) and Yb(Nb+1)=Yb(1). ! Nb=2*(Jmax-Jmin+Imax-Imin) shft=1-Imin do i=Imin,Imax-1 Xb(i+shft)=Xgrd(i,Jmin) Yb(i+shft)=Ygrd(i,Jmin) enddo shft=1-Jmin+Imax-Imin do j=Jmin,Jmax-1 Xb(j+shft)=Xgrd(Imax,j) Yb(j+shft)=Ygrd(Imax,j) enddo shft=1+Jmax-Jmin+2*Imax-Imin do i=Imax,Imin+1,-1 Xb(shft-i)=Xgrd(i,Jmax) Yb(shft-i)=Ygrd(i,Jmax) enddo shft=1+2*Jmax-Jmin+2*(Imax-Imin) do j=Jmax,Jmin+1,-1 Xb(shft-j)=Xgrd(Imin,j) Yb(shft-j)=Ygrd(Imin,j) enddo ! !--------------------------------------------------------------------- ! Check if point (Xo,Yo) is inside of the defined polygon. !--------------------------------------------------------------------- ! Try_Range=inside(Xo,Yo,Xb,Yb,Nb) return end function inside (Xo,Yo,Xb,Yb,Nb) ! !================================================ Hernan G. Arango === ! Copyright (c) 2000 Rutgers/UCLA ! !======================================== Alexander F. Shchepetkin === ! ! ! Given the vectors Xb and Yb of size Nb, defining the coordinates ! ! of a closed polygon, this function find if the point (Xo,Yo) is ! ! inside the polygon. If the point (Xo,Yo) falls exactly on the ! ! boundary of the polygon, it still considered inside. ! ! ! ! This algorithm does not rely on the setting of Xb(Nb)=Xb(1) and ! ! Yb(Nb)=Yb(1). Instead, it assumes that the last closing segment ! ! is (Xb(Nb),Yb(Nb)) --> (Xb(1),Yb(1)). ! ! ! ! Reference: ! ! ! ! Reid, C., 1969: A long way from Euclid. Oceanography EMR, ! ! page 174. ! ! ! ! Algorithm: ! ! ! ! The decision whether the point is inside or outside the polygon ! ! is done by counting the number of crossings from the ray (Xo,Yo) ! ! to (Xo,-infinity), hereafter called meridian, by the boundary of ! ! the polygon. In this counting procedure, a crossing is counted ! ! as +2 if the crossing happens from "left to right" or -2 if from ! ! "right to left". If the counting adds up to zero, then the point ! ! is outside. Otherwise, it is either inside or on the boundary. ! ! ! ! This routine is a modified version of the Reid (1969) algorithm, ! ! where all crossings were counted as positive and the decision is ! ! made based on whether the number of crossings is even or odd. ! ! This new algorithm may produce different results in cases where ! ! Xo accidentally coinsides with one of the (Xb(k),k=1:Nb) points. ! ! In this case, the crossing is counted here as +1 or -1 depending ! ! of the sign of (Xb(k+1)-Xb(k)). Crossings are not counted if ! ! Xo=Xb(k)=Xb(k+1). Therefore, if Xo=Xb(k0) and Yo>Yb(k0), and if ! ! Xb(k0-1) < Xb(k0) < Xb(k0+1), the crossing is counted twice but ! ! with weight +1 (for segments with k=k0-1 and k=k0). Similarly if ! ! Xb(k0-1) > Xb(k0) > Xb(k0+1), the crossing is counted twice with ! ! weight -1 each time. If, on the other hand, the meridian only ! ! touches the boundary, that is, for example, Xb(k0-1) < Xb(k0)=Xo ! ! and Xb(k0+1) < Xb(k0)=Xo, then the crossing is counted as +1 for ! ! segment k=k0-1 and -1 for segment k=k0, resulting in no crossing. ! ! ! ! Note 1: (Explanation of the logical condition) ! ! ! ! Suppose that there exist two points (x1,y1)=(Xb(k),Yb(k)) and ! ! (x2,y2)=(Xb(k+1),Yb(k+1)), such that, either (x1 < Xo < x2) or ! ! (x1 > Xo > x2). Therefore, meridian x=Xo intersects the segment ! ! (x1,y1) -> (x2,x2) and the ordinate of the point of intersection ! ! is: ! ! ! ! y1*(x2-Xo) + y2*(Xo-x1) ! ! y = ----------------------- ! ! x2-x1 ! ! ! ! The mathematical statement that point (Xo,Yo) either coinsides ! ! with the point of intersection or lies to the north (Yo>=y) from ! ! it is, therefore, equivalent to the statement: ! ! ! ! Yo*(x2-x1) >= y1*(x2-Xo) + y2*(Xo-x1), if x2-x1 > 0 ! ! or ! ! Yo*(x2-x1) <= y1*(x2-Xo) + y2*(Xo-x1), if x2-x1 < 0 ! ! ! ! which, after noting that Yo*(x2-x1) = Yo*(x2-Xo + Xo-x1) may be ! ! rewritten as: ! ! ! ! (Yo-y1)*(x2-Xo) + (Yo-y2)*(Xo-x1) >= 0, if x2-x1 > 0 ! ! or ! ! (Yo-y1)*(x2-Xo) + (Yo-y2)*(Xo-x1) <= 0, if x2-x1 < 0 ! ! ! ! and both versions can be merged into essentially the condition ! ! that (Yo-y1)*(x2-Xo)+(Yo-y2)*(Xo-x1) has the same sign as x2-x1. ! ! That is, the product of these two must be positive or zero. ! ! ! !===================================================================== ! implicit none ! logical inside integer Nb, Nstep, crossings, i, inc, k, kk, nc parameter (Nstep=128) integer index(Nstep) real Xb(Nb+1), Yb(Nb+1), Xo, Yo, dx1, dx2, dxy ! !--------------------------------------------------------------------- ! Find intersections. !--------------------------------------------------------------------- ! ! Set crossings counter and close the contour of the polygon. ! crossings=0 Xb(Nb+1)=Xb(1) Yb(Nb+1)=Yb(1) ! ! The search is optimized. First select the indices of segments ! where Xb(k) is different from Xb(k+1) and Xo falls between them. ! Then, further investigate these segments in a separate loop. ! Doing it in two stages takes less time because the first loop is ! pipelined. ! do kk=0,Nb-1,Nstep nc=0 do k=kk+1,MIN(kk+Nstep,Nb) if (((Xb(k+1)-Xo)*(Xo-Xb(k)).ge.0.0).and. & (Xb(k).ne.Xb(k+1))) then nc=nc+1 index(nc)=k endif enddo do i=1,nc k=index(i) if (Xb(k).ne.Xb(k+1)) then dx1=Xo-Xb(k) dx2=Xb(k+1)-Xo dxy=dx2*(Yo-Yb(k))-dx1*(Yb(k+1)-Yo) inc=0 if ((Xb(k).eq.Xo).and.(Yb(k).eq.Yo)) then crossings=1 goto 10 elseif (((dx1.eq.0.0).and.(Yo.ge.Yb(k ))).or. & ((dx2.eq.0.0).and.(Yo.ge.Yb(k+1)))) then inc=1 elseif ((dx1*dx2.gt.0.0).and. ! See Note 1 & ((Xb(k+1)-Xb(k))*dxy.ge.0.0)) then inc=2 endif if (Xb(k+1).gt.Xb(k)) then crossings=crossings+inc else crossings=crossings-inc endif endif enddo enddo ! ! Determine if point (Xo,Yo) is inside of closed polygon. ! 10 if (crossings.eq.0) then inside=.false. else inside=.true. endif return end #endif /* FLOATS */