|
XBeach
|
00001 module interp 00002 implicit none 00003 private 00004 public linear_interp, interp_in_cyclic_function, interp_using_trapez_rule 00005 public grmap2, mkmap, grmap, linear_interp_2d, make_map, mkmap_step, trapezoidal 00006 save 00007 contains 00008 pure subroutine linear_interp_2d(X,nx,Y,ny,Z,xx,yy,zz,method,exception) 00009 00010 implicit none 00011 ! input/output 00012 integer,intent(in) :: nx,ny 00013 real*8,dimension(nx),intent(in) :: X 00014 real*8,dimension(ny),intent(in) :: Y 00015 real*8,dimension(nx,ny),intent(in) :: Z 00016 real*8,intent(in) :: xx,yy 00017 real*8,intent(out) :: zz 00018 character(len=*),intent(in) :: method 00019 real*8,intent(in) :: exception 00020 ! internal 00021 integer,dimension(4) :: ind 00022 real*8,dimension(2) :: yint 00023 real*8 :: modx,mody,disx,disy 00024 logical :: interpX,interpY 00025 00026 ! does the interpolation point fall within the data? 00027 if (xx>=minval(X) .and. xx<=maxval(X)) then 00028 interpX = .true. 00029 else 00030 interpX = .false. 00031 endif 00032 if (yy>=minval(Y) .and. yy<=maxval(Y)) then 00033 interpY = .true. 00034 else 00035 interpY = .false. 00036 endif 00037 00038 if (interpX .and. interpY) then 00039 ! find rank position xx in X direction 00040 ind(1) = minloc(X,1,X>=xx) 00041 ind(2) = maxloc(X,1,X<=xx) 00042 ! find rank position yy in Y direction 00043 ind(3) = minloc(Y,1,Y>=yy) 00044 ind(4) = maxloc(Y,1,Y<=yy) 00045 ! distance between X points and Y points 00046 disx = X(ind(1))-X(ind(2)) 00047 disy = Y(ind(3))-Y(ind(4)) 00048 ! relative position of (xx,yy) on disx,disy 00049 if (disx>0.d0) then 00050 modx = (xx-X(ind(2)))/disx 00051 else 00052 modx = 0.d0 ! xx corresponds exactly to X point 00053 endif 00054 if (disy>0.d0) then 00055 mody = (yy-Y(ind(4)))/disy 00056 else 00057 mody = 0.d0 ! yy corresponds exactly to Y point 00058 endif 00059 ! interpolate the correct Y value, based on two nearest X intersects 00060 ! if disy==0 then only single interpolation needed. This could also be done for X, 00061 ! but since this is used by waveparams and the interp angles (Y) are more often equal 00062 ! this is probably faster. 00063 if (disy>0.d0) then 00064 yint(1) = (1.d0-modx)*Z(ind(2),ind(3))+modx*Z(ind(1),ind(3)) 00065 yint(2) = (1.d0-modx)*Z(ind(2),ind(4))+modx*Z(ind(1),ind(4)) 00066 zz = (1.d0-mody)*yint(2)+mody*yint(1) 00067 else 00068 zz = (1.d0-modx)*Z(ind(2),ind(3))+modx*Z(ind(1),ind(3)) 00069 endif 00070 else 00071 select case (method) 00072 case ('interp') 00073 zz = exception 00074 case ('extendclosest') 00075 ! find closest X point 00076 ind(1) = minloc(abs(X-xx),1) 00077 ! find closest Y point 00078 ind(3) = minloc(abs(Y-yy),1) 00079 ! external value given that of closest point 00080 zz = Z(ind(1),ind(3)) 00081 case default 00082 zz = exception 00083 endselect 00084 endif 00085 00086 end subroutine linear_interp_2d 00087 00088 ! 00089 ! NAME 00090 ! linear_interp 00091 ! SYNOPSIS 00092 ! Interpolate linearly into an array Y given the value of X. 00093 ! Return both the interpolated value and the position in the 00094 ! array. 00095 ! 00096 ! ARGUMENTS 00097 ! - X: independent array (sorted in ascending order) 00098 ! - Y: array whose values are a function of X 00099 ! - XX: specified value, interpolating in first array 00100 ! - YY: interpolation result 00101 ! - INDINT: (optional) index in array (x(indint) <= xx <= x(indint+1)) 00102 ! 00103 ! SOURCE 00104 ! 00105 pure subroutine linear_interp(x, y, n, xx, yy, indint) 00106 integer, intent(in) :: n 00107 real*8, dimension(n), intent(in) :: x 00108 real*8, dimension(n), intent(in) :: y 00109 real*8, intent(in) :: xx 00110 real*8, intent(out) :: yy 00111 integer, intent(out), optional :: indint 00112 !**** 00113 ! 00114 ! CODE: linear interpolation 00115 ! 00116 ! 00117 ! 00118 real*8 :: a, b, dyy 00119 integer :: j 00120 00121 yy = 0.0d0 00122 if ( present(indint) ) then 00123 indint = 0 00124 endif 00125 00126 if (n.le.0) return 00127 ! 00128 ! *** N GREATER THAN 0 00129 ! 00130 if (n.eq.1) then 00131 yy = y(1) 00132 return 00133 endif 00134 00135 call binary_search( x, n, xx, j ) 00136 00137 if ( j .le. 0 ) then 00138 yy = y(1) 00139 elseif ( j .ge. n ) then 00140 yy = y(n) 00141 else 00142 a = x (j+1) 00143 b = x (j) 00144 if ( a .eq. b ) then 00145 dyy = 0.0d0 00146 else 00147 dyy = (y(j+1) - y(j)) / (a - b) 00148 endif 00149 yy = y(j) + (xx - x(j)) * dyy 00150 endif 00151 00152 if ( present(indint) ) then 00153 indint = j 00154 endif 00155 00156 return 00157 00158 end subroutine linear_interp 00159 00160 !****f* Interpolation/binary_search 00161 ! 00162 ! NAME 00163 ! binary_search 00164 ! SYNOPSIS 00165 ! Perform a binary search in an ordered real array 00166 ! to find the largest entry equal or lower than a given value: 00167 ! Given an array XX of length N, given value X, return a value J 00168 ! such that X is between XX(J) en XX (J+1) 00169 ! 00170 ! XX must be monotonic, either decreasing or increasing 00171 ! J=0 or J=N indicates X is out of range. 00172 ! 00173 ! ARGUMENTS 00174 ! - XX: ordered array of values 00175 ! - X: value to be found 00176 ! - J: index such that X is between XX(J) and XX(J+1) 00177 ! 00178 ! SOURCE 00179 ! 00180 pure subroutine binary_search(xx, n, x, j) 00181 integer, intent(in) :: N 00182 real*8, dimension(N), intent(in) :: xx 00183 real*8, intent(in) :: x 00184 integer, intent(out) :: j 00185 !**** 00186 ! 00187 ! CODE: binary search in (real) arrays 00188 ! 00189 ! Requirement: 00190 ! Parameter wp set to the proper kind 00191 ! 00192 ! Subroutine from 'Numerical recipes' Fortran edition. 00193 ! Given an array XX of length N, given value X, return a value J 00194 ! such that X is between XX(J) en XX (J+1) 00195 ! XX must be monotonic, either decreasing or increasin 00196 ! J=0 or J=N indicates X is out of range. 00197 00198 00199 ! 00200 ! Local variables 00201 ! 00202 integer jl, ju, jm 00203 logical l1, l2 00204 00205 jl = 0 00206 ju = n+1 00207 10 if (ju-jl .gt. 1) then 00208 jm = (ju+jl)/2 00209 l1 = xx(n) .gt. xx(1) 00210 l2 = x .gt. xx(jm) 00211 if ( (l1.and.l2) .or. (.not. (l1 .or. l2)) ) then 00212 jl = jm 00213 else 00214 ju = jm 00215 endif 00216 goto 10 00217 endif 00218 00219 j = jl 00220 00221 return 00222 00223 end subroutine binary_search 00224 00225 subroutine mkmap(code ,x1 ,y1 ,m1 ,n1 , & 00226 & x2 ,y2 ,n2 ,xs ,ys , & 00227 & nrx ,nry ,iflag ,nrin ,w , & 00228 & iref ,iprint ,covered ,xymiss) 00229 !!--description----------------------------------------------------------------- 00230 ! 00231 ! 00232 ! SUBROUTINE MKMAP 00233 ! Interpolation of curvilinear, numerically ordered grid (grid1) 00234 ! to random points, with weighting points (grid2). 00235 ! 00236 ! J.A. Roelvink 00237 ! Deltares 00238 ! 24-2-1992 (MKMAP 00239 ! 00240 ! Given: numerically ordered grid M1*N1 00241 ! with coordinates X1 (1:M1,1:N1) 00242 ! and Y1 (1:M1,1:N1) 00243 ! 00244 ! Also given: random points X2(1:N2) 00245 ! and Y2(1:N2) 00246 ! 00247 ! To be determined: weighting factors and pointers for bilinear interpolation 00248 ! Weighting factors and locations of the requested (random) points in the 00249 ! ordered grid are saved in resp. 00250 ! W(1:4,1:N2) and Iref(1:4,1:N2) 00251 ! 00252 !!--pseudo code and references-------------------------------------------------- 00253 ! NONE 00254 !!--declarations---------------------------------------------------------------- 00255 ! 00256 implicit none 00257 ! 00258 ! Global variables 00259 ! 00260 integer , intent(in) :: iprint 00261 integer , intent(in) :: m1 00262 integer , intent(in) :: n1 00263 integer :: n2 00264 integer , dimension(4 , n2), intent(out) :: iref 00265 integer , dimension(m1, n1), intent(in) :: code 00266 integer , dimension(n2) :: covered ! 0: target point is not covered by source grid (default) 00267 ! 1: target point is covered by valid points of source grid 00268 ! -1: target point is covered by invalid points of source grid 00269 integer , dimension(n2) :: iflag 00270 integer , dimension(n2) :: nrin 00271 integer , dimension(n2) :: nrx 00272 integer , dimension(n2) :: nry 00273 real , intent(in) :: xymiss ! missing value in grid 1 00274 real*8 , dimension(4 , n2) :: w ! Remains single precision! 00275 real*8 , dimension(m1, n1), intent(in) :: x1 00276 real*8 , dimension(m1, n1), intent(in) :: y1 00277 real*8 , dimension(n2) :: x2 00278 real*8 , dimension(n2) :: xs 00279 real*8 , dimension(n2) :: y2 00280 real*8 , dimension(n2) :: ys 00281 ! 00282 ! Local variables 00283 ! 00284 integer :: i 00285 integer :: i1 00286 integer :: i2 00287 integer :: ier 00288 integer :: iin 00289 integer :: inout 00290 integer :: ip 00291 integer :: j1 00292 integer :: lomaxx 00293 integer :: lominx 00294 integer :: lomnx 00295 integer :: lomny 00296 integer :: m1max 00297 integer :: m1min 00298 integer :: n1max 00299 integer :: n1min 00300 integer :: nin 00301 real*8 :: eps 00302 real*8 :: xpmax 00303 real*8 :: xpmean 00304 real*8 :: xpmin 00305 real*8 :: ypmax 00306 real*8 :: ypmean 00307 real*8 :: ypmin 00308 real*8 , dimension(5) :: xp 00309 real*8 , dimension(5) :: yp 00310 real*8 , dimension(4) :: hbuff 00311 ! 00312 !! executable statements ------------------------------------------------------- 00313 ! 00314 eps = 0.00001 00315 ! 00316 ! Initialise tables 00317 ! 00318 if (iprint==1) write (*, *) 'in mkmap', m1, n1, n2 00319 lomnx = 1 00320 lomny = 1 00321 m1min = 1 00322 m1max = m1 00323 n1min = 1 00324 n1max = n1 00325 nrx=0 00326 nry=0 00327 iflag=0 00328 nrin=0 00329 xs=0.d0 00330 ys=0.d0 00331 iref=0 00332 w=0.d0 00333 ! 00334 ! Sort X2 en Y2 00335 ! 00336 call sort(n2 ,x2 ,xs ,nrx ) 00337 call sort(n2 ,y2 ,ys ,nry ) 00338 ! 00339 ! Loop over all cels of grid1 00340 ! 00341 ! 00342 do j1 = n1min, n1max - 1 00343 do i1 = m1min, m1max - 1 00344 ! 00345 ! Cell definition 00346 ! 00347 xp(1) = x1(i1, j1) 00348 xp(2) = x1(i1 + 1, j1) 00349 xp(3) = x1(i1 + 1, j1 + 1) 00350 xp(4) = x1(i1, j1 + 1) 00351 yp(1) = y1(i1, j1) 00352 yp(2) = y1(i1 + 1, j1) 00353 yp(3) = y1(i1 + 1, j1 + 1) 00354 yp(4) = y1(i1, j1 + 1) 00355 if ( (xp(1) > xymiss-eps .and. xp(1) < xymiss+eps) & 00356 & .or. (xp(2) > xymiss-eps .and. xp(2) < xymiss+eps) & 00357 & .or. (xp(3) > xymiss-eps .and. xp(3) < xymiss+eps) & 00358 & .or. (xp(4) > xymiss-eps .and. xp(4) < xymiss+eps) & 00359 & .or. (yp(1) > xymiss-eps .and. yp(1) < xymiss+eps) & 00360 & .or. (yp(2) > xymiss-eps .and. yp(2) < xymiss+eps) & 00361 & .or. (yp(3) > xymiss-eps .and. yp(3) < xymiss+eps) & 00362 & .or. (yp(4) > xymiss-eps .and. yp(4) < xymiss+eps) ) cycle 00363 ! 00364 ! Determine minimum and maximum X and Y of the cell 00365 ! 00366 xpmin = 1.e10 00367 xpmax = -1.e10 00368 ypmin = 1.e10 00369 ypmax = -1.e10 00370 do ip = 1, 4 00371 xpmin = min(xp(ip), xpmin) 00372 xpmax = max(xp(ip), xpmax) 00373 ypmin = min(yp(ip), ypmin) 00374 ypmax = max(yp(ip), ypmax) 00375 enddo 00376 xpmean = .5*(xpmin + xpmax) 00377 ypmean = .5*(ypmin + ypmax) 00378 ! 00379 ! First selection of points of grid2 that can be located in the cell 00380 ! 00381 ! Find centre of the cell in tables Xs and Ys 00382 ! 00383 call hunt(xs ,n2 ,xpmean ,lomnx ) 00384 call hunt(ys ,n2 ,ypmean ,lomny ) 00385 ! 00386 ! For points with X-values between Xpmin and Xpmax set: iFlag(i)=1 00387 ! 00388 lominx = lomnx 00389 lomaxx = lomnx 00390 do i = lomnx, 1, -1 00391 if (xs(i)>=xpmin) then 00392 lominx = i 00393 iflag(nrx(i)) = 1 00394 else 00395 exit 00396 endif 00397 enddo 00398 do i = lomnx + 1, n2 00399 if (xs(i)<=xpmax) then 00400 lomaxx = i 00401 iflag(nrx(i)) = 1 00402 else 00403 exit 00404 endif 00405 enddo 00406 ! 00407 ! For the points with Y-values between Ypmin and Ypmax, 00408 ! that also lie between Xpmin and Xpmax: Save them in NrIn 00409 ! 00410 iin = 1 00411 do i = lomny, 1, -1 00412 if (ys(i)>=ypmin) then 00413 nrin(iin) = nry(i)*iflag(nry(i)) 00414 iin = iin + iflag(nry(i)) 00415 else 00416 exit 00417 endif 00418 enddo 00419 do i = lomny + 1, n2 00420 if (ys(i)<=ypmax) then 00421 nrin(iin) = nry(i)*iflag(nry(i)) 00422 iin = iin + iflag(nry(i)) 00423 else 00424 exit 00425 endif 00426 enddo 00427 nin = iin - 1 00428 ! 00429 ! Put iFlag back to 0 00430 ! 00431 do i = lominx, lomaxx 00432 if (i/=0) iflag(nrx(i)) = 0 00433 enddo 00434 ! 00435 ! Check whether selected points of grid2 lie within the cell 00436 ! using subroutine IPON; if so, determine weights W of the surrounding 00437 ! values in grid1 using subroutine INTRP4. Save the weights in Wtab 00438 ! The reference to grid1 is saved in arrays Iref and Jref. 00439 ! 00440 do iin = 1, nin 00441 i2 = nrin(iin) 00442 inout = -1 00443 call ipon(xp ,yp ,4 ,x2(i2) , & 00444 & y2(i2) ,inout ) 00445 if (inout>=0) then 00446 ! 00447 ! Check on point with nonsense information on source grid 1; 00448 ! This depends on agreements 00449 ! 00450 if ( (code(i1 , j1 )==1 .or. code(i1 , j1 )==3) & 00451 & .and. (code(i1 + 1, j1 )==1 .or. code(i1 + 1, j1 )==3) & 00452 & .and. (code(i1 + 1, j1 + 1)==1 .or. code(i1 + 1, j1 + 1)==3) & 00453 & .and. (code(i1 , j1 + 1)==1 .or. code(i1 , j1 + 1)==3) ) then 00454 ! 00455 ! grid2 point is covered by 4 valid grid1 points 00456 ! 00457 covered(i2) = 1 00458 call bilin5(xp, yp, x2(i2), y2(i2), hbuff, ier) 00459 w(:, i2) = real(hbuff(:),4) 00460 ! 00461 iref(1, i2) = i1 + (j1 - 1)*m1 00462 iref(2, i2) = i1 + 1 + (j1 - 1)*m1 00463 iref(3, i2) = i1 + 1 + j1*m1 00464 iref(4, i2) = i1 + j1*m1 00465 elseif (code(i1, j1)>0 .and. code(i1 + 1, j1)>0 .and. & 00466 & code(i1 + 1, j1 + 1)>0 .and. code(i1, j1 + 1)>0) then 00467 ! 00468 ! Grid2 point is covered by 4 valid grid1 points, containing 00469 ! one or more boundary points: 00470 ! - do not use grid1 values 00471 ! - use extrapolation values 00472 ! 00473 covered(i2) = 0 00474 else 00475 ! 00476 ! Grid2 point is not covered by 4 valid grid1 points: 00477 ! - do not use grid1 values 00478 ! - do not use extrapolation values 00479 ! 00480 covered(i2) = -1 00481 endif 00482 endif 00483 enddo 00484 enddo 00485 enddo 00486 end subroutine mkmap 00487 00488 subroutine make_map(code, x1, y1, m1, n1, x2, y2, n2, w, iref, & 00489 & iprint ,covered ,xymiss) 00490 !!--description----------------------------------------------------------------- 00491 ! 00492 ! 00493 ! SUBROUTINE MKMAP 00494 ! Interpolation of curvilinear, numerically ordered grid (grid1) 00495 ! to random points, with weighting points (grid2). 00496 ! 00497 ! J.A. Roelvink 00498 ! Deltares 00499 ! 24-2-1992 (MKMAP 00500 ! 00501 ! Given: numerically ordered grid M1*N1 00502 ! with coordinates X1 (1:M1,1:N1) 00503 ! and Y1 (1:M1,1:N1) 00504 ! 00505 ! Also given: random points X2(1:N2) 00506 ! and Y2(1:N2) 00507 ! 00508 ! To be determined: weighting factors and pointers for bilinear interpolation 00509 ! Weighting factors and locations of the requested (random) points in the 00510 ! ordered grid are saved in resp. 00511 ! W(1:4,1:N2) and Iref(1:4,1:N2) 00512 ! 00513 !!--pseudo code and references-------------------------------------------------- 00514 ! NONE 00515 !!--declarations---------------------------------------------------------------- 00516 ! 00517 implicit none 00518 ! 00519 ! Global variables 00520 ! 00521 integer , intent(in) :: iprint 00522 integer , intent(in) :: m1 00523 integer , intent(in) :: n1 00524 integer , intent(in) :: n2 00525 integer , dimension(4 , n2), intent(out) :: iref 00526 real*8 , dimension(4 , n2), intent(out) :: w ! Remains single precision! 00527 integer , dimension(m1, n1), intent(in) :: code 00528 integer , dimension(n2) , intent(out) :: covered ! 0: target point is not covered by source grid (default) 00529 ! 1: target point is covered by valid points of source grid 00530 ! -1: target point is covered by invalid points of source grid 00531 real , intent(in) :: xymiss ! missing value in grid 1 00532 real*8 , dimension(m1, n1), intent(in) :: x1 00533 real*8 , dimension(m1, n1), intent(in) :: y1 00534 real*8 , dimension(n2), intent(in) :: x2 00535 real*8 , dimension(n2), intent(in) :: y2 00536 ! 00537 ! Local variables 00538 ! 00539 integer, dimension(:),allocatable :: iflag 00540 integer, dimension(:),allocatable :: nrin 00541 integer, dimension(:),allocatable :: nrx 00542 integer, dimension(:),allocatable :: nry 00543 real*8, dimension(:),allocatable :: xs 00544 real*8, dimension(:),allocatable :: ys 00545 00546 integer :: i 00547 integer :: i1 00548 integer :: i2 00549 integer :: ier 00550 integer :: iin 00551 integer :: inout 00552 integer :: ip 00553 integer :: j1 00554 integer :: lomaxx 00555 integer :: lominx 00556 integer :: lomnx 00557 integer :: lomny 00558 integer :: m1max 00559 integer :: m1min 00560 integer :: n1max 00561 integer :: n1min 00562 integer :: nin 00563 real*8 :: eps 00564 real*8 :: xpmax 00565 real*8 :: xpmean 00566 real*8 :: xpmin 00567 real*8 :: ypmax 00568 real*8 :: ypmean 00569 real*8 :: ypmin 00570 real*8 , dimension(5) :: xp 00571 real*8 , dimension(5) :: yp 00572 real*8 , dimension(4) :: hbuff 00573 ! 00574 allocate(iflag(n2)) 00575 allocate(nrin(n2)) 00576 allocate(nrx(n2)) 00577 allocate(nry(n2)) 00578 allocate(xs(n2)) 00579 allocate(ys(n2)) 00580 !! executable statements ------------------------------------------------------- 00581 ! 00582 eps = 0.00001 00583 ! 00584 ! Initialise tables 00585 ! 00586 if (iprint==1) write (*, *) 'in mkmap', m1, n1, n2 00587 lomnx = 1 00588 lomny = 1 00589 m1min = 1 00590 m1max = m1 00591 n1min = 1 00592 n1max = n1 00593 nrx=0 00594 nry=0 00595 iflag=0 00596 nrin=0 00597 xs=0.d0 00598 ys=0.d0 00599 iref=0 00600 w=0.d0 00601 ! 00602 ! Sort X2 en Y2 00603 ! 00604 call sort(n2 ,x2 ,xs ,nrx ) 00605 call sort(n2 ,y2 ,ys ,nry ) 00606 ! 00607 ! Loop over all cels of grid1 00608 ! 00609 ! 00610 do j1 = n1min, n1max - 1 00611 do i1 = m1min, m1max - 1 00612 ! 00613 ! Cell definition 00614 ! 00615 xp(1) = x1(i1, j1) 00616 xp(2) = x1(i1 + 1, j1) 00617 xp(3) = x1(i1 + 1, j1 + 1) 00618 xp(4) = x1(i1, j1 + 1) 00619 yp(1) = y1(i1, j1) 00620 yp(2) = y1(i1 + 1, j1) 00621 yp(3) = y1(i1 + 1, j1 + 1) 00622 yp(4) = y1(i1, j1 + 1) 00623 if ( (xp(1) > xymiss-eps .and. xp(1) < xymiss+eps) & 00624 & .or. (xp(2) > xymiss-eps .and. xp(2) < xymiss+eps) & 00625 & .or. (xp(3) > xymiss-eps .and. xp(3) < xymiss+eps) & 00626 & .or. (xp(4) > xymiss-eps .and. xp(4) < xymiss+eps) & 00627 & .or. (yp(1) > xymiss-eps .and. yp(1) < xymiss+eps) & 00628 & .or. (yp(2) > xymiss-eps .and. yp(2) < xymiss+eps) & 00629 & .or. (yp(3) > xymiss-eps .and. yp(3) < xymiss+eps) & 00630 & .or. (yp(4) > xymiss-eps .and. yp(4) < xymiss+eps) ) cycle 00631 ! 00632 ! Determine minimum and maximum X and Y of the cell 00633 ! 00634 xpmin = 1.e10 00635 xpmax = -1.e10 00636 ypmin = 1.e10 00637 ypmax = -1.e10 00638 do ip = 1, 4 00639 xpmin = min(xp(ip), xpmin) 00640 xpmax = max(xp(ip), xpmax) 00641 ypmin = min(yp(ip), ypmin) 00642 ypmax = max(yp(ip), ypmax) 00643 enddo 00644 xpmean = .5*(xpmin + xpmax) 00645 ypmean = .5*(ypmin + ypmax) 00646 ! 00647 ! First selection of points of grid2 that can be located in the cell 00648 ! 00649 ! Find centre of the cell in tables Xs and Ys 00650 ! 00651 call hunt(xs ,n2 ,xpmean ,lomnx ) 00652 call hunt(ys ,n2 ,ypmean ,lomny ) 00653 ! 00654 ! For points with X-values between Xpmin and Xpmax set: iFlag(i)=1 00655 ! 00656 lominx = lomnx 00657 lomaxx = lomnx 00658 do i = lomnx, 1, -1 00659 if (xs(i)>=xpmin) then 00660 lominx = i 00661 iflag(nrx(i)) = 1 00662 else 00663 exit 00664 endif 00665 enddo 00666 do i = lomnx + 1, n2 00667 if (xs(i)<=xpmax) then 00668 lomaxx = i 00669 iflag(nrx(i)) = 1 00670 else 00671 exit 00672 endif 00673 enddo 00674 ! 00675 ! For the points with Y-values between Ypmin and Ypmax, 00676 ! that also lie between Xpmin and Xpmax: Save them in NrIn 00677 ! 00678 iin = 1 00679 do i = lomny, 1, -1 00680 if (ys(i)>=ypmin) then 00681 nrin(iin) = nry(i)*iflag(nry(i)) 00682 iin = iin + iflag(nry(i)) 00683 else 00684 exit 00685 endif 00686 enddo 00687 do i = lomny + 1, n2 00688 if (ys(i)<=ypmax) then 00689 nrin(iin) = nry(i)*iflag(nry(i)) 00690 iin = iin + iflag(nry(i)) 00691 else 00692 exit 00693 endif 00694 enddo 00695 nin = iin - 1 00696 ! 00697 ! Put iFlag back to 0 00698 ! 00699 do i = lominx, lomaxx 00700 if (i/=0) iflag(nrx(i)) = 0 00701 enddo 00702 ! 00703 ! Check whether selected points of grid2 lie within the cell 00704 ! using subroutine IPON; if so, determine weights W of the surrounding 00705 ! values in grid1 using subroutine INTRP4. Save the weights in Wtab 00706 ! The reference to grid1 is saved in arrays Iref and Jref. 00707 ! 00708 do iin = 1, nin 00709 i2 = nrin(iin) 00710 inout = -1 00711 call ipon(xp ,yp ,4 ,x2(i2) , & 00712 & y2(i2) ,inout ) 00713 if (inout>=0) then 00714 ! 00715 ! Check on point with nonsense information on source grid 1; 00716 ! This depends on agreements 00717 ! 00718 if ( (code(i1 , j1 )==1 .or. code(i1 , j1 )==3) & 00719 & .and. (code(i1 + 1, j1 )==1 .or. code(i1 + 1, j1 )==3) & 00720 & .and. (code(i1 + 1, j1 + 1)==1 .or. code(i1 + 1, j1 + 1)==3) & 00721 & .and. (code(i1 , j1 + 1)==1 .or. code(i1 , j1 + 1)==3) ) then 00722 ! 00723 ! grid2 point is covered by 4 valid grid1 points 00724 ! 00725 covered(i2) = 1 00726 call bilin5(xp, yp, x2(i2), y2(i2), hbuff, ier) 00727 w(:, i2) = real(hbuff(:),4) 00728 ! 00729 iref(1, i2) = i1 + (j1 - 1)*m1 00730 iref(2, i2) = i1 + 1 + (j1 - 1)*m1 00731 iref(3, i2) = i1 + 1 + j1*m1 00732 iref(4, i2) = i1 + j1*m1 00733 elseif (code(i1, j1)>0 .and. code(i1 + 1, j1)>0 .and. & 00734 & code(i1 + 1, j1 + 1)>0 .and. code(i1, j1 + 1)>0) then 00735 ! 00736 ! Grid2 point is covered by 4 valid grid1 points, containing 00737 ! one or more boundary points: 00738 ! - do not use grid1 values 00739 ! - use extrapolation values 00740 ! 00741 covered(i2) = 0 00742 else 00743 ! 00744 ! Grid2 point is not covered by 4 valid grid1 points: 00745 ! - do not use grid1 values 00746 ! - do not use extrapolation values 00747 ! 00748 covered(i2) = -1 00749 endif 00750 endif 00751 enddo 00752 enddo 00753 enddo 00754 00755 end subroutine make_map 00756 00757 subroutine mkmap_step (x1 ,y1 ,m1 ,n1 , & 00758 & alfaz ,dsu ,dnv , & 00759 & x2 ,y2 ,n2 ,w , iref ) 00760 !!--description----------------------------------------------------------------- 00761 ! 00762 ! 00763 ! SUBROUTINE MKMAP_STEP 00764 ! Interpolation of curvilinear, numerically ordered grid (grid1) 00765 ! to random points, with weighting points (grid2). 00766 ! 00767 ! J.A. Roelvink 00768 ! UNESCO-IHE / Deltares 00769 ! 8 May 2013 00770 ! 00771 ! Given: numerically ordered grid M1*N1 00772 ! with (new)coordinates X1 (1:M1,1:N1) 00773 ! and Y1 (1:M1,1:N1) 00774 ! 00775 ! Also given: random points X2(1:N2) 00776 ! and Y2(1:N2) 00777 ! previous iref 00778 ! 00779 ! To be determined: weighting factors and updated pointers for bilinear interpolation 00780 ! Weighting factors and locations of the requested (random) points in the 00781 ! ordered grid are saved in resp. 00782 ! W(1:4,1:N2) and Iref(1:4,1:N2) 00783 ! 00784 !!--pseudo code and references-------------------------------------------------- 00785 ! NONE 00786 !!--declarations---------------------------------------------------------------- 00787 ! 00788 implicit none 00789 ! 00790 ! Global variables 00791 ! 00792 integer , intent(in) :: m1 00793 integer , intent(in) :: n1 00794 integer :: n2 00795 integer , dimension(4 , n2) :: iref 00796 real*8 , dimension(4 , n2) :: w 00797 real*8 , dimension(m1, n1), intent(in) :: x1 00798 real*8 , dimension(m1, n1), intent(in) :: y1 00799 real*8 , dimension(m1, n1), intent(in) :: alfaz 00800 real*8 , dimension(m1, n1), intent(in) :: dsu 00801 real*8 , dimension(m1, n1), intent(in) :: dnv 00802 real*8 , dimension(n2) :: x2 00803 real*8 , dimension(n2) :: y2 00804 real*8 :: dx 00805 real*8 :: dy 00806 real*8 :: ds 00807 real*8 :: dn 00808 real*8 :: dsrel 00809 real*8 :: dnrel 00810 integer :: i1,j1,i2 00811 00812 do i2=1,n2 00813 if (iref(1,i2)>0) then 00814 i1=mod(iref(1,i2)-1,m1)+1 00815 j1=(iref(1,i2)-i1)/m1+1 00816 dx=x2(i2)-x1(i1,j1) 00817 dy=y2(i2)-y1(i1,j1) 00818 ds= dx*cos(alfaz(i1,j1))+dy*sin(alfaz(i1,j1)) 00819 dn=-dx*sin(alfaz(i1,j1))+dy*cos(alfaz(i1,j1)) 00820 dsrel=ds/dsu(i1,j1) 00821 dnrel=dn/dnv(i1,j1) 00822 i1=i1+floor(dsrel) 00823 j1=j1+floor(dnrel) 00824 dsrel=dsrel-floor(dsrel) 00825 dnrel=dnrel-floor(dnrel) 00826 iref(1, i2) = i1 + (j1 - 1)*m1 00827 iref(2, i2) = i1 + 1 + (j1 - 1)*m1 00828 iref(3, i2) = i1 + 1 + j1*m1 00829 iref(4, i2) = i1 + j1*m1 00830 w(1, i2) = (1.d0-dsrel)*(1.d0-dnrel) 00831 w(2, i2) = dsrel *(1.d0-dnrel) 00832 w(3, i2) = dsrel * dnrel 00833 w(4, i2) = (1.d0-dsrel)* dnrel 00834 endif 00835 enddo 00836 00837 end subroutine mkmap_step 00838 00839 subroutine grmap(f1 ,n1 ,f2 ,n2 ,iref , & 00840 & w ,np ,iprint ) 00841 !!--description----------------------------------------------------------------- 00842 ! 00843 ! compute interpolated values for all points on grid 2 00844 ! 00845 ! special treatment of points on grid 2 that are outside 00846 ! grid 1; in that case iref(1,i2)=0 AND w(ip,i2)=0 for all ip 00847 ! 00848 ! Iref(1,i2) i1 ifac F2(i2)*ifac Result 00849 ! 00850 ! 0 1 1 F2(i2) Old value is kept 00851 ! j,j>0 j 0 0. F2 is initialized 00852 ! 00853 !!--pseudo code and references-------------------------------------------------- 00854 ! NONE 00855 !!--declarations---------------------------------------------------------------- 00856 implicit none 00857 ! 00858 ! Global variables 00859 ! 00860 integer , intent(in) :: iprint 00861 integer , intent(in) :: n1 00862 integer , intent(in) :: n2 00863 integer , intent(in) :: np 00864 integer, dimension(np, n2), intent(in) :: iref 00865 real*8 , dimension(n1) , intent(in) :: f1 00866 real*8 , dimension(n2) :: f2 00867 real*8 , dimension(np, n2), intent(in) :: w 00868 ! 00869 ! Local variables 00870 ! 00871 integer :: i 00872 integer :: i1 00873 integer :: i2 00874 integer :: ip 00875 ! 00876 !! executable statements ------------------------------------------------------- 00877 ! 00878 if (iprint==1) write (*, *) 'in grmap n1 n2', n1, n2 00879 do i2 = 1, n2 00880 i = iref(1, i2) 00881 if (i>0) then 00882 f2(i2)=0.d0 00883 ! i1 = max(i, 1) 00884 ! ifac = 1 - i/i1 00885 ! f2(i2) = f2(i2)*ifac 00886 ! 00887 ! Function values at grid 2 are expressed as weighted average 00888 ! of function values in Np surrounding points of grid 1 00889 ! 00890 if (iprint==1 .and. i2<=n2) & 00891 & write (*, '(1X,A,I6,4(1X,E10.4))') ' i2 w ', i2, (w(ip, i2) , ip = 1, & 00892 & np) 00893 do ip = 1, np 00894 i = iref(ip, i2) 00895 i1 = max(i, 1) 00896 if (iprint==1 .and. i2<=n2) write (*, *) ' i1,f1(i1) ', i1, f1(i1) 00897 f2(i2) = f2(i2) + w(ip, i2)*f1(i1) 00898 enddo 00899 endif 00900 enddo 00901 end subroutine grmap 00902 00903 subroutine grmap2(f1, cellsz1i , n1 ,f2 ,cellsz2, n2 ,iref , & 00904 & w ,np ) 00905 !!--description----------------------------------------------------------------- 00906 ! 00907 ! compute interpolated values for all points on grid 1 given reference table 00908 ! for grid 2; this works the other way round from GRMAP. Assumption is that 00909 ! grid 2 is much finer than grid 1. For each point in grid 2 we know the 00910 ! surrounding points in grid 1 and the related weights. Instead of using this to 00911 ! interpolate from 1 to 2 we now integrate from 2 to 1 using the same weights. 00912 ! 00913 ! 00914 !!--pseudo code and references-------------------------------------------------- 00915 ! NONE 00916 !!--declarations---------------------------------------------------------------- 00917 implicit none 00918 ! 00919 ! Global variables 00920 ! 00921 integer , intent(in) :: n1 00922 integer , intent(in) :: n2 00923 integer , intent(in) :: np 00924 integer, dimension(np, n2), intent(in) :: iref 00925 real*8 , dimension(n1) :: f1 00926 real*8 , dimension(n1) , intent(in) :: cellsz1i !array with 1/cell size 00927 real*8 , intent(in) :: cellsz2 00928 real*8 , dimension(n2) , intent(in) :: f2 00929 real*8 , dimension(np, n2), intent(in) :: w 00930 ! 00931 ! Local variables 00932 ! 00933 integer :: i1 00934 integer :: i2 00935 integer :: ip 00936 ! 00937 !! executable statements ------------------------------------------------------- 00938 ! 00939 do ip = 1, np 00940 do i2=1,n2 00941 i1 = iref(ip, i2) 00942 if (i1>0) then 00943 !f2(i2) = f2(i2) + w(ip, i2)*f1(i1) 00944 f1(i1) = f1(i1) + w(ip, i2)*f2(i2)*cellsz2*cellsz1i(i1) 00945 endif 00946 enddo 00947 enddo 00948 end subroutine grmap2 00949 00950 subroutine ipon(xq ,yq ,n ,xp ,yp ,inout ) 00951 !--description---------------------------------------------------------------- 00952 ! 00953 ! Deltares * 00954 ! AUTHOR : J.A.ROELVINK * 00955 ! DATE : 22-12-1988 * 00956 ! DETERMINE WHETHER POINT (xp,yp) LIES IN POLYGON (x,y) OF n POINTS * 00957 ! POINT n+1 IS SET EQUAL TO POINT 1 * 00958 ! (ARRAY MUST HAVE DIMENSION n+1 IN MAIN PROGRAMME * 00959 ! inpout = -1 : OUTSIDE POLYGON * 00960 ! inpout = 0 : ON EDGE OF POLYGON * 00961 ! inpout = 1 : INSIDE POLYGON * 00962 ! USED METHOD : - DRAW A VERTICAL LINE THROUGH (xp,yp) * 00963 ! - DETERMINE NUMBER OF INTERSECTIONS WITH POLYGON * 00964 ! UNDER yp : nunder * 00965 ! - IF nunder IS EVEN, THEN THE POINT LIES OUTSIDE * 00966 ! THE POLYGON, OTHERWISE IT LIES INSIDE * 00967 ! - THE EDGE IS TREATED SEPARATELY * 00968 ! 00969 !--pseudo code and references------------------------------------------------- 00970 ! NONE 00971 !--declarations--------------------------------------------------------------- 00972 ! 00973 implicit none 00974 ! 00975 ! Global variables 00976 ! 00977 integer , intent(out) :: inout 00978 integer , intent(in) :: n 00979 real*8 , intent(in) :: xp 00980 real*8 , intent(in) :: yp 00981 real*8 , dimension(*) :: xq 00982 real*8 , dimension(*) :: yq 00983 ! 00984 ! Local variables 00985 ! 00986 integer :: i 00987 integer :: ierr 00988 integer :: nunder 00989 real*4 :: ysn 00990 real*4 , dimension(:), allocatable :: x 00991 real*4 , dimension(:), allocatable :: y 00992 ! 00993 ! executable statements ------------------------------------------------------ 00994 ! 00995 allocate(x(n+1)) 00996 allocate(y(n+1)) 00997 do i = 1, n 00998 x(i) = real( xq(i)-xp , 4) 00999 y(i) = real( yq(i)-yp , 4) 01000 enddo 01001 x(n + 1) = x(1) 01002 y(n + 1) = y(1) 01003 nunder = 0 01004 do i = 1, n 01005 if ((x(i)<0. .and. x(i + 1)>=0.).or.(x(i + 1)<0. .and. x(i)>=0.)) then 01006 if (y(i)<0. .and. y(i + 1)<0.) then 01007 nunder = nunder + 1 01008 elseif ((y(i)<=0. .and. y(i + 1)>=0.) .or. & 01009 & (y(i + 1)<=0. .and. y(i)>=0.)) then 01010 ysn = (y(i)*x(i + 1) - x(i)*y(i + 1))/(x(i + 1) - x(i)) 01011 if (ysn<0.) then 01012 nunder = nunder + 1 01013 elseif (ysn<=0.) then 01014 ! 01015 ! Edge 01016 ! 01017 inout = 0 01018 goto 100 01019 else 01020 endif 01021 else 01022 endif 01023 elseif (abs(x(i))<1.0E-8 .and. abs(x(i + 1))<1.0E-8) then 01024 if ((y(i)<=0. .and. y(i + 1)>=0.).or.(y(i + 1)<=0..and.y(i)>=0.)) & 01025 & then 01026 ! 01027 ! Edge 01028 ! 01029 inout = 0 01030 goto 100 01031 endif 01032 else 01033 endif 01034 enddo 01035 if (mod(nunder, 2)==0) then 01036 ! 01037 ! Outside 01038 ! 01039 inout = -1 01040 else 01041 ! 01042 ! Inside 01043 ! 01044 inout = 1 01045 endif 01046 100 continue 01047 deallocate(x, stat=ierr) 01048 deallocate(y, stat=ierr) 01049 end subroutine ipon 01050 01051 subroutine hunt(xx ,n ,x ,jlo ) 01052 !!--description----------------------------------------------------------------- 01053 ! NONE 01054 !!--pseudo code and references-------------------------------------------------- 01055 ! NONE 01056 !!--declarations---------------------------------------------------------------- 01057 ! 01058 implicit none 01059 ! 01060 ! Global variables 01061 ! 01062 integer :: jlo 01063 integer , intent(in) :: n 01064 real*8 , intent(in) :: x 01065 real*8, dimension(n), intent(in) :: xx 01066 ! 01067 ! Local variables 01068 ! 01069 integer :: inc 01070 integer :: jhi 01071 integer :: jm 01072 logical :: ascnd 01073 ! 01074 !! executable statements ------------------------------------------------------- 01075 ! 01076 ascnd = xx(n)>=xx(1) 01077 if (jlo<=0 .or. jlo>n) then 01078 jlo = 0 01079 jhi = n + 1 01080 goto 3 01081 endif 01082 inc = 1 01083 if (x>=xx(jlo) .eqv. ascnd) then 01084 1 continue 01085 jhi = jlo + inc 01086 if (jhi>n) then 01087 jhi = n + 1 01088 elseif (x>=xx(jhi) .eqv. ascnd) then 01089 jlo = jhi 01090 inc = inc + inc 01091 goto 1 01092 else 01093 endif 01094 else 01095 jhi = jlo 01096 2 continue 01097 jlo = jhi - inc 01098 if (jlo<1) then 01099 jlo = 0 01100 elseif (x<xx(jlo) .eqv. ascnd) then 01101 jhi = jlo 01102 inc = inc + inc 01103 goto 2 01104 else 01105 endif 01106 endif 01107 3 continue 01108 if (jhi - jlo==1) then 01109 return 01110 endif 01111 jm = (jhi + jlo)/2 01112 if (x>xx(jm) .eqv. ascnd) then 01113 jlo = jm 01114 else 01115 jhi = jm 01116 endif 01117 goto 3 01118 end subroutine hunt 01119 01120 subroutine indexx(n ,arrin ,indx ) 01121 !----- GPL --------------------------------------------------------------------- 01122 ! 01123 ! Copyright (C) Stichting Deltares, 2011. 01124 ! 01125 ! This program is free software: you can redistribute it and/or modify 01126 ! it under the terms of the GNU General Public License as published by 01127 ! the Free Software Foundation version 3. 01128 ! 01129 ! This program is distributed in the hope that it will be useful, 01130 ! but WITHOUT ANY WARRANTY; without even the implied warranty of 01131 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 01132 ! GNU General Public License for more details. 01133 ! 01134 ! You should have received a copy of the GNU General Public License 01135 ! along with this program. If not, see <http://www.gnu.org/licenses/>. 01136 ! 01137 ! contact: delft3d.support@deltares.nl 01138 ! Stichting Deltares 01139 ! P.O. Box 177 01140 ! 2600 MH Delft, The Netherlands 01141 ! 01142 ! All indications and logos of, and references to, "Delft3D" and "Deltares" 01143 ! are registered trademarks of Stichting Deltares, and remain the property of 01144 ! Stichting Deltares. All rights reserved. 01145 ! 01146 !------------------------------------------------------------------------------- 01147 ! $Id: interp.F90 5598 2019-12-09 10:21:29Z ridde_mo $ 01148 ! $HeadURL: https://svn.oss.deltares.nl/repos/xbeach/trunk/src/xbeachlibrary/interp.F90 $ 01149 !!--description----------------------------------------------------------------- 01150 ! NONE 01151 !!--pseudo code and references-------------------------------------------------- 01152 ! NONE 01153 !!--declarations---------------------------------------------------------------- 01154 ! 01155 implicit none 01156 ! 01157 ! Global variables 01158 ! 01159 integer , intent(in) :: n 01160 integer, dimension(n) :: indx 01161 real*8 , dimension(n), intent(in) :: arrin 01162 ! 01163 ! Local variables 01164 ! 01165 integer :: i 01166 integer :: indxt 01167 integer :: ir 01168 integer :: j 01169 integer :: l 01170 real*8 :: q 01171 ! 01172 !! executable statements ------------------------------------------------------- 01173 ! 01174 do j = 1, n 01175 indx(j) = j 01176 enddo 01177 l = n/2 + 1 01178 ir = n 01179 10 continue 01180 if (l>1) then 01181 l = l - 1 01182 indxt = indx(l) 01183 q = arrin(indxt) 01184 else 01185 indxt = indx(ir) 01186 q = arrin(indxt) 01187 indx(ir) = indx(1) 01188 ir = ir - 1 01189 if (ir==1) then 01190 indx(1) = indxt 01191 return 01192 endif 01193 endif 01194 i = l 01195 j = l + l 01196 20 continue 01197 if (j<=ir) then 01198 if (j<ir) then 01199 if (arrin(indx(j))<arrin(indx(j + 1))) j = j + 1 01200 endif 01201 if (q<arrin(indx(j))) then 01202 indx(i) = indx(j) 01203 i = j 01204 j = j + j 01205 else 01206 j = ir + 1 01207 endif 01208 goto 20 01209 endif 01210 indx(i) = indxt 01211 goto 10 01212 end subroutine indexx 01213 01214 subroutine bilin5(xa, ya, x0, y0, w, ier) 01215 !!--description----------------------------------------------------------------- 01216 ! NONE 01217 !!--pseudo code and references-------------------------------------------------- 01218 ! 01219 ! Author: H. Petit 01220 ! 01221 !!--declarations---------------------------------------------------------------- 01222 ! 01223 implicit none 01224 ! 01225 ! Global variables 01226 ! 01227 integer , intent(out) :: ier 01228 real*8 , intent(in) :: x0 01229 real*8 , intent(in) :: y0 01230 real*8 , dimension(4), intent(out) :: w 01231 real*8 , dimension(4), intent(in) :: xa 01232 real*8 , dimension(4), intent(in) :: ya 01233 ! 01234 ! Local variables 01235 ! 01236 real*8 :: a 01237 real*8 :: a21 01238 real*8 :: a22 01239 real*8 :: a31 01240 real*8 :: a32 01241 real*8 :: a41 01242 real*8 :: a42 01243 real*8 :: b 01244 real*8 :: c 01245 real*8 :: det 01246 real*8 :: discr 01247 real*8 :: eta 01248 real*8 :: x 01249 real*8 :: x1 01250 real*8 :: x2 01251 real*8 :: x3 01252 real*8 :: x3t 01253 real*8 :: x4 01254 real*8 :: xi 01255 real*8 :: xt 01256 real*8 :: y 01257 real*8 :: y1 01258 real*8 :: y2 01259 real*8 :: y3 01260 real*8 :: y3t 01261 real*8 :: y4 01262 real*8 :: yt 01263 ! 01264 !! executable statements ------------------------------------------------------- 01265 ! 01266 ! read(12,*)x1,y1,f1 01267 x1 = xa(1) 01268 y1 = ya(1) 01269 ! read(12,*)x2,y2,f2 01270 x2 = xa(2) 01271 y2 = ya(2) 01272 ! read(12,*)x3,y3,f3 01273 x3 = xa(3) 01274 y3 = ya(3) 01275 ! read(12,*)x4,y4,f4 01276 x4 = xa(4) 01277 y4 = ya(4) 01278 x = x0 01279 y = y0 01280 ! 01281 ! The bilinear interpolation problem is first transformed 01282 ! to the quadrangle with nodes (0,0),(1,0),(x3t,y3t),(0,1) 01283 ! and required location (xt,yt) 01284 ! 01285 a21 = x2 - x1 01286 a22 = y2 - y1 01287 a31 = x3 - x1 01288 a32 = y3 - y1 01289 a41 = x4 - x1 01290 a42 = y4 - y1 01291 det = a21*a42 - a22*a41 01292 if (abs(det) < 1.0e-20) then 01293 ier = 1 01294 goto 99999 01295 endif 01296 x3t = ( a42*a31 - a41*a32 ) / det 01297 y3t = ( -a22*a31 + a21*a32 ) / det 01298 xt = ( a42*(x - x1) - a41*(y - y1)) / det 01299 yt = ( -a22*(x - x1) + a21*(y - y1)) / det 01300 if ((x3t < 0.0) .or. (y3t < 0.0)) then 01301 ! write (*, *) 'distorted quadrangle' 01302 ier = 1 01303 goto 99999 01304 endif 01305 if (abs(x3t - 1.0d0) < 1.0e-7) then 01306 xi = xt 01307 if (abs(y3t - 1.0d0) < 1.0e-7) then 01308 eta = yt 01309 elseif (abs(1.0d0 + (y3t - 1.0d0)*xt) < 1.0e-6) then 01310 ! write (*, *) 'extrapolation over too large a distance' 01311 ier = 1 01312 goto 99999 01313 else 01314 eta = yt / (1.0d0 + (y3t - 1.0d0)*xt) 01315 endif 01316 elseif (abs(y3t - 1.0d0) < 1.0e-6) then 01317 eta = yt 01318 if (abs(1.0d0 + (x3t - 1.0d0)*yt) < 1.0e-6) then 01319 ! write (*, *) 'extrapolation over too large a distance' 01320 ier = 1 01321 goto 99999 01322 else 01323 xi = xt / (1.0d0 + (x3t - 1.0d0)*yt) 01324 endif 01325 else 01326 a = y3t - 1.0d0 01327 b = 1.0d0 + (x3t - 1.0d0)*yt - (y3t - 1.0d0)*xt 01328 c = -xt 01329 discr = b*b - 4.0d0*a*c 01330 if (discr < 1.0e-6) then 01331 ! write (*, *) 'extrapolation over too large a distance' 01332 ier = 1 01333 goto 99999 01334 endif 01335 xi = (-b + sqrt(discr)) / (2.0d0*a) 01336 eta = ((y3t - 1.0d0)*(xi - xt) + (x3t - 1.0d0)*yt) / (x3t - 1.0d0) 01337 endif 01338 w(1) = (1.0d0-xi) * (1.0d0-eta) 01339 w(2) = xi * (1.0d0-eta) 01340 w(3) = xi * eta 01341 w(4) = eta * (1.0d0-xi ) 01342 return 01343 99999 continue 01344 end subroutine bilin5 01345 01346 subroutine sort(n ,ra ,wksp ,iwksp ) 01347 !!--description----------------------------------------------------------------- 01348 ! Sorts an array, routine from Numerical Recipes 01349 !!--pseudo code and references-------------------------------------------------- 01350 ! NONE 01351 !!--declarations---------------------------------------------------------------- 01352 ! 01353 implicit none 01354 ! 01355 ! Global variables 01356 ! 01357 integer :: n 01358 integer, dimension(n) :: iwksp 01359 real*8 , dimension(n) :: ra 01360 real*8 , dimension(n), intent(out) :: wksp 01361 ! 01362 ! Local variables 01363 ! 01364 integer :: j 01365 ! 01366 !! executable statements ------------------------------------------------------- 01367 ! 01368 call indexx(n ,ra ,iwksp ) 01369 do j = 1, n 01370 wksp(j) = ra(iwksp(j)) 01371 enddo 01372 end subroutine sort 01373 01374 subroutine trapezoidal(x,y,n,x1_in,x2_in,integ) 01375 ! Compute integral over function y(x) from x1 to x2 01376 ! x is equidistant 01377 !(c)2014 Dano Roelvink 01378 implicit none 01379 integer, intent(in) :: n 01380 real*8, dimension(n), intent(in) :: x,y ! arrays x and y(x) 01381 real*8, intent(in) :: x1_in,x2_in ! integration limits 01382 real*8, intent(out) :: integ ! integral 01383 real*8 :: x1,y1,x2,y2,Ifirst,Imid,Ilast,dx 01384 integer :: i1,i2,i,i1p1,i2p1 01385 if (x1_in>x(n).or.x2_in<x(1)) then 01386 integ=0 01387 else 01388 01389 x1=max(x1_in,x(1)+1d-60) 01390 x2=min(x2_in,x(n)-1d-60) 01391 dx=(x(n)-x(1))/(n-1) 01392 i1=floor((x1-x(1))/dx)+1 01393 i2=floor((x2-x(1))/dx)+1 01394 i1p1=min(i1+1,n) 01395 i2p1=min(i2+1,n) 01396 ! first partial trapezoid 01397 y1=y(i1)+(x1-x(i1))/dx*(y(i1p1)-y(i1)) 01398 Ifirst=.5*(x(i1p1)-x1)*(y(i1p1)+y1) 01399 ! middle part 01400 Imid=.5*y(i1p1) 01401 do i=i1+2,i2-1 01402 Imid=Imid+y(i) 01403 end do 01404 Imid=Imid+.5*y(i2) 01405 Imid=Imid*dx 01406 ! last partial trapezoid 01407 y2=y(i2)+(x2-x(i2))/dx*(y(i2p1)-y(i2)) 01408 Ilast=.5*(x2-x(i2))*(y2+y(i2)) 01409 integ=Ifirst+Imid+Ilast 01410 endif 01411 end subroutine trapezoidal 01412 01413 subroutine trapezoidal_cyclic(x,y,n,xcycle,x1,x2,integ) 01414 ! Compute integral over function y(x) from x1 to x2 01415 ! x is equidistant, function is cyclic so y(x+k*xcycle)=y(x) 01416 !(c)2014 Dano Roelvink 01417 implicit none 01418 integer, intent(in) :: n 01419 real*8, dimension(n), intent(in) :: x,y ! arrays x and y(x) 01420 real*8, intent(in) :: x1,x2,xcycle ! integration limits,cycle length 01421 real*8, intent(out) :: integ ! integral 01422 real*8, dimension(:),allocatable :: xp,yp 01423 real*8 :: dx 01424 integer :: ip,indt,np 01425 01426 dx=x(2)-x(1) 01427 np=floor(x2/dx)-(floor(x1/dx)+1)+2 01428 allocate(xp(np)) 01429 allocate(yp(np)) 01430 xp(1)=x1 01431 do ip=2,np-1 01432 xp(ip)=(floor(x1/dx)+ip-1)*dx 01433 end do 01434 xp(np)=x2 01435 if (xcycle>0) then 01436 call interp_in_cyclic_function(x,y,n,xcycle,xp,np,yp) 01437 else 01438 do ip=1,np 01439 call linear_interp (x,y,n,xp(ip),yp(ip),indt) 01440 end do 01441 endif 01442 integ=0.d0 01443 do ip=1,np-1 01444 integ=integ+.5*(xp(ip+1)-xp(ip))*(yp(ip+1)+yp(ip)) 01445 end do 01446 01447 deallocate(xp) 01448 deallocate(yp) 01449 01450 end subroutine trapezoidal_cyclic 01451 01452 subroutine interp_using_trapez_rule(x,y,n,xcycle,xtarg,ytarg,ntarg) 01453 implicit none 01454 ! Compute integral over function y(x) from x1 to x2 01455 ! x is equidistant, function is cyclic so y(x+k*xcycle)=y(x) 01456 !(c)2014 Dano Roelvink 01457 integer, intent(in) :: n,ntarg 01458 real*8, dimension(n), intent(in) :: x,y ! arrays x and y(x) 01459 real*8, intent(in) :: xcycle ! cycle length 01460 real*8, dimension(ntarg), intent(in) :: xtarg ! x values to interpolate to 01461 real*8, dimension(ntarg), intent(out) :: ytarg ! y values to interpolate to 01462 real*8 :: dx,x1,x2,integ 01463 integer :: itarg 01464 01465 dx=xtarg(2)-xtarg(1) 01466 do itarg=1,ntarg 01467 x1=xtarg(itarg)-.5*dx 01468 x2=xtarg(itarg)+.5*dx 01469 call trapezoidal_cyclic(x,y,n,xcycle,x1,x2,integ) 01470 ytarg(itarg)=integ/dx 01471 end do 01472 01473 end subroutine interp_using_trapez_rule 01474 01475 subroutine interp_in_cyclic_function(x,y,n,xcycle,xp,np,yp) 01476 implicit none 01477 integer, intent(in) :: n 01478 real*8, dimension(n), intent(in) :: x,y ! arrays x and y(x) 01479 real*8, dimension(:), allocatable :: xc,yc ! complemented cyclic arrays 01480 real*8, intent(in) :: xcycle ! cycle length 01481 integer, intent(in) :: np 01482 real*8, dimension(np), intent(in) :: xp ! points to interpolate to 01483 real*8, dimension(np), intent(out) :: yp ! interpolated yp values 01484 integer :: icycle,ip,ileft,iright,nc,i 01485 real*8 :: dx,yleft,yright,facleft,facright 01486 01487 dx=x(2)-x(1) 01488 icycle=nint(xcycle/dx) 01489 allocate(xc(icycle+1)) 01490 allocate(yc(icycle+1)) 01491 if (n>icycle+1) then 01492 ! nonsense; error 01493 elseif (n==icycle+1) then 01494 ! e.g. 0, 30, 60,...360 01495 xc(1:n)=x 01496 yc(1:n)=y 01497 elseif (n==icycle) then 01498 ! e.g. 30, 60 ...360 or 15, 45 ...345 with n=12 and icycle=12 01499 ! interpolate between n-th and 1st value 01500 xc(1:n)=x 01501 yc(1:n)=y 01502 xc(n+1)=xc(n)+dx 01503 yc(n+1)=yc(1) 01504 elseif (n<icycle) then 01505 ! do not interpolate between n-th and 1st value; set values in gap to 0 01506 xc(1:n)=x 01507 yc(1:n)=y 01508 do i=n+1,icycle 01509 xc(i)=xc(i-1)+dx 01510 yc(i)=0.d0 01511 end do 01512 xc(icycle+1)=xc(icycle)+dx 01513 yc(icycle+1)=yc(1) 01514 endif 01515 01516 nc=icycle+1 01517 01518 do ip=1,np 01519 ileft=floor((xp(ip)-xc(1))/dx)+1 01520 do while (ileft<1) 01521 ileft=ileft+icycle 01522 end do 01523 do while (ileft>nc) 01524 ileft=ileft-icycle 01525 end do 01526 if (ileft>nc .or. ileft<1) then 01527 yleft=0 01528 else 01529 yleft=yc(ileft) 01530 endif 01531 iright=ileft+1 01532 do while (iright<1) 01533 iright=iright+icycle 01534 end do 01535 do while (iright>nc) 01536 iright=iright-icycle 01537 end do 01538 if (iright>nc .or. iright<1) then 01539 yright=0 01540 else 01541 yright=yc(iright) 01542 endif 01543 facright=mod(xp(ip),dx)/dx 01544 facleft=1.d0-facright 01545 yp(ip)=facleft*yleft+facright*yright 01546 01547 end do 01548 deallocate(xc) 01549 deallocate(yc) 01550 01551 end subroutine interp_in_cyclic_function 01552 01553 end module interp
1.7.6.1
