ROMS
Loading...
Searching...
No Matches
Data Types | Functions/Subroutines
interpolate.F File Reference
#include "cppdefs.h"
Include dependency graph for interpolate.F:

Go to the source code of this file.

Data Types

type  roms_interp_type
 

Functions/Subroutines

program __interpolate_f__
 
subroutine linterp2d (ng, lbx, ubx, lby, uby, xinp, yinp, finp, lbi, ubi, lbj, ubj, istr, iend, jstr, jend, iout, jout, xout, yout, fout, minval, maxval)
 
subroutine cinterp2d (ng, lbx, ubx, lby, uby, xinp, yinp, finp, lbi, ubi, lbj, ubj, istr, iend, jstr, jend, iout, jout, xout, yout, fout, minval, maxval)
 
subroutine hindices (ng, lbi, ubi, lbj, ubj, is, ie, js, je, angler, xgrd, ygrd, lbm, ubm, lbn, ubn, ms, me, ns, ne, xpos, ypos, ipos, jpos, ijspv, rectangular)
 
logical function try_range (ng, lbi, ubi, lbj, ubj, xgrd, ygrd, imin, imax, jmin, jmax, xo, yo)
 
logical function inside (nb, xb, yb, xo, yo)
 
subroutine roms_datum_interp_2d (s, fld_src, datum, method, obsvetting)
 
subroutine roms_datum_interp_3d (s, fld_src, z_src, z_locs, datum, method, obsvetting)
 
subroutine roms_datum_column_interp (s, fld_src, datum, method)
 
subroutine roms_horiz_interp_2d (s, fld_src, fld_dst, method)
 
subroutine roms_horiz_interp_3d (s, fld_src, fld_dst, method)
 
subroutine roms_interp_create (ng, tile, gtype, s)
 
subroutine roms_interp_delete (s)
 
subroutine roms_interp_fractional (s)
 

Function/Subroutine Documentation

◆ __interpolate_f__()

program __interpolate_f__

Definition at line 3 of file interpolate.F.

References strings_mod::founderror(), and mod_iounits::stdout.

Here is the call graph for this function:

◆ cinterp2d()

subroutine __interpolate_f__::cinterp2d ( integer, intent(in) ng,
integer, intent(in) lbx,
integer, intent(in) ubx,
integer, intent(in) lby,
integer, intent(in) uby,
real(r8), dimension(lbx:ubx,lby:uby), intent(in) xinp,
real(r8), dimension(lbx:ubx,lby:uby), intent(in) yinp,
real(r8), dimension(lbx:ubx,lby:uby), intent(in) finp,
integer, intent(in) lbi,
integer, intent(in) ubi,
integer, intent(in) lbj,
integer, intent(in) ubj,
integer, intent(in) istr,
integer, intent(in) iend,
integer, intent(in) jstr,
integer, intent(in) jend,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) iout,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) jout,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) xout,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) yout,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) fout,
real(r8), intent(out), optional minval,
real(r8), intent(out), optional maxval )
private

Definition at line 281 of file interpolate.F.

288!
289!=======================================================================
290! !
291! Given any gridded 2D field, Finp, at locations (Xinp,Yinp) this !
292! routine performs bicubic interpolation at locations (Xout,Yout). !
293! To facilitate the interpolation within any irregularly gridded !
294! field, the fractional grid cell indices (Iout,Jout) with respect !
295! Finp are needed at input. Notice that the routine "hindices" can !
296! be used to compute these indices. !
297! !
298! On Input: !
299! !
300! ng Nested grid number. !
301! LBx I-dimension lower bound of gridded field, Finp. !
302! UBx I-dimension upper bound of gridded field, Finp. !
303! LBy J-dimension lower bound of gridded field, Finp. !
304! UBy J-dimension upper bound of gridded field, Finp. !
305! Xinp X-locations of gridded field, Finp. !
306! Yinp Y-locations of gridded field, Finp. !
307! Finp 2D field to interpolate from. !
308! LBi I-dimension Lower bound of data to interpolate, Fout. !
309! UBi I-dimension Upper bound of data to interpolate, Fout. !
310! LBj J-dimension Lower bound of data to interpolate, Fout. !
311! UBj J-dimension Upper bound of data to interpolate, Fout. !
312! Istr Starting data I-index to interpolate, Fout. !
313! Iend Ending data I-index to interpolate, Fout. !
314! Jstr Starting data J-index to interpolate, Fout. !
315! Jend Ending data J-index to interpolate, Fout. !
316! Iout I-fractional Xinp grid cell containing Xout. !
317! Jout J-fractional Yinp grid cell containing Yout. !
318! Xout X-locations to interpolate, Fout. !
319! Yout Y-locations to interpolate, Fout. !
320! !
321! On Output: !
322! !
323! Fout Interpolated 2D field. !
324! MinVal Interpolated field minimum value (OPTIONAL). !
325! MaxVal Interpolated field maximum value (OPTIONAL). !
326! !
327!=======================================================================
328!
329! Imported variable declarations.
330!
331 integer, intent(in) :: ng, LBx, UBx, LBy, UBy
332 integer, intent(in) :: LBi, UBi, LBj, UBj
333 integer, intent(in) :: Istr, Iend, Jstr, Jend
334!
335 real(r8), intent(in) :: Xinp(LBx:UBx,LBy:UBy)
336 real(r8), intent(in) :: Yinp(LBx:UBx,LBy:UBy)
337 real(r8), intent(in) :: Finp(LBx:UBx,LBy:UBy)
338
339 real(r8), intent(in) :: Iout(LBi:UBi,LBj:UBj)
340 real(r8), intent(in) :: Jout(LBi:UBi,LBj:UBj)
341 real(r8), intent(in) :: Xout(LBi:UBi,LBj:UBj)
342 real(r8), intent(in) :: Yout(LBi:UBi,LBj:UBj)
343
344 real(r8), intent(out) :: Fout(LBi:UBi,LBj:UBj)
345!
346 real(r8), intent(out), optional :: MinVal
347 real(r8), intent(out), optional :: MaxVal
348!
349! Local variable declarations.
350!
351 integer i, ic, iter, i1, i2, j, jc, j1, j2
352!
353 real(r8) :: a11, a12, a21, a22
354 real(r8) :: e11, e12, e21, e22
355 real(r8) :: cff, d1, d2, dfc, dx, dy, eta, xi, xy, yx
356 real(r8) :: f0, fx, fxx, fxxx, fxxy, fxy, fxyy, fy, fyy, fyyy
357!
358 real(r8), parameter :: C01 = 1.0_r8/48.0_r8
359 real(r8), parameter :: C02 = 1.0_r8/32.0_r8
360 real(r8), parameter :: C03 = 0.0625_r8 ! 1/16
361 real(r8), parameter :: C04 = 1.0_r8/6.0_r8
362 real(r8), parameter :: C05 = 0.25_r8
363 real(r8), parameter :: C06 = 0.5_r8
364 real(r8), parameter :: C07 = 0.3125_r8 ! 5/16
365 real(r8), parameter :: C08 = 0.625_r8 ! 5/8
366 real(r8), parameter :: C09 = 1.5_r8
367 real(r8), parameter :: C10 = 13.0_r8/24.0_r8
368
369 real(r8), parameter :: LIMTR = 3.0_r8
370 real(r8), parameter :: spv = 0.0_r8 ! HGA need work
371
372 real(r8) :: Fmin, Fmax
373!
374 real(r8), dimension(-1:2,-1:2) :: dfx, dfy, f
375!
376!-----------------------------------------------------------------------
377! Interpolates requested field locations (Xout,Yout).
378!-----------------------------------------------------------------------
379!
380 fmin=1.0e+35_r8
381 fmax=-1.0e+35_r8
382 DO j=jstr,jend
383 DO i=istr,iend
384 i1=int(iout(i,j))
385 i2=i1+1
386 j1=int(jout(i,j))
387 j2=j1+1
388 IF (((lbx.le.i1).and.(i1.le.ubx)).and. &
389 & ((lby.le.j1).and.(j1.le.uby))) THEN
390!
391! Determine local fractional coordinates (xi,eta) corresponding to
392! the target point (Xout,Yout) on the grid (Xinp,Yinp). Here, "xi"
393! and "eta" are defined, in such a way, that xi=eta=0 corresponds
394! to the middle of the cell (i1:i1+1,j1:j1+1), while xi=+/-1/2 and
395! eta=+/-1/2 (any combination +/- signs) corresponds to the four
396! corner points of the cell. Inside the cell it is assumed that
397! (Xout,Yout) are expressed via bi-linear functions of (xi,eta),
398! where term proportional to xi*eta does not vanish because
399! coordinate transformation may be at least weakly non-orthogonal
400! due to discretization errors. The associated non-linear system
401! is solved by iterative method of Newton.
402!
403 xy=xinp(i2,j2)-xinp(i1,j2)-xinp(i2,j1)+xinp(i1,j1)
404 yx=yinp(i2,j2)-yinp(i1,j2)-yinp(i2,j1)+yinp(i1,j1)
405 dx=xout(i,j)-0.25_r8*(xinp(i2,j2)+xinp(i1,j2)+ &
406 & xinp(i2,j1)+xinp(i1,j1))
407 dy=yout(i,j)-0.25_r8*(yinp(i2,j2)+yinp(i1,j2)+ &
408 & yinp(i2,j1)+yinp(i1,j1))
409!
410! The coordinate transformation matrix:
411!
412! e11 e12
413! e21 e22
414!
415! contains derivatives of (Xinp,Yinp) with respect to (xi,eta). Because
416! the coordinates may be non-orthogonal (at least due to discretization
417! errors), the nonlinear system
418!
419! e11*xi+e12*eta+xy*xi*eta=dx
420! e21*xi+e22*eta+yx*xi*eta=dy
421!
422! needs to be solved in order to retain symmetry.
423!
424 e11=0.5_r8*(xinp(i2,j2)-xinp(i1,j2)+xinp(i2,j1)-xinp(i1,j1))
425 e12=0.5_r8*(xinp(i2,j2)+xinp(i1,j2)-xinp(i2,j1)-xinp(i1,j1))
426 e21=0.5_r8*(yinp(i2,j2)-yinp(i1,j2)+yinp(i2,j1)-yinp(i1,j1))
427 e22=0.5_r8*(yinp(i2,j2)+yinp(i1,j2)-yinp(i2,j1)-yinp(i1,j1))
428!
429 cff=1.0_r8/(e11*e22-e12*e21)
430 xi=cff*(e22*dx-e12*dy)
431 eta=cff*(e11*dy-e21*dx)
432!
433 DO iter=1,4
434 d1=dx-e11*xi-e12*eta-xy*xi*eta
435 d2=dy-e21*xi-e22*eta-yx*xi*eta
436 a11=e11+xy*eta
437 a12=e12+xy*xi
438 a21=e21+yx*eta
439 a22=e22+yx*xi
440 cff=1.0_r8/(a11*a22-a12*a21)
441 xi =xi +cff*(a22*d1-a12*d2)
442 eta=eta+cff*(a11*d2-a21*d1)
443 END DO
444
445#ifndef CUBIC_MASKED
446!
447! Genuinely two-dimensional, isotropic cubic interpolation scheme
448! using 12-point stencil. In the code below the interpolated field,
449! Fout, is expanded into two-dimensional Taylor series of local
450! fractional coordinates "xi" and "eta", retaining all terms of
451! combined power up to third order (that is, xi, eta, xi^2, eta^2,
452! xi*eta, xi^3, eta^3, xi^2*eta, and xi*eta^2), with all
453! coefficients (i.e, derivatives) computed via x x
454! two-dimensional finite difference expressions | |
455! of "natural" order of accuracy: 4th-order for x--x--x--x
456! the field itself and its first derivatives in | |
457! both directions; and 2nd-order for all higher- x--x--x--x
458! order derivatives. The permissible range of | |
459! of coordinates is -1/2 < xi,eta < +1/2, which x--x
460! covers the central cell on the stencil, while
461! xi=eta=0 corresponds to its center. This interpolation scheme has
462! the property that if xi,eta=+/-1/2 (any combination of +/- signs)
463! it reproduces exactly value of the function at the corresponding
464! corner of the central "working" cell. However, it does not pass
465! exactly through the extreme points of the stencil, where either
466! xi=+/-3/2 or eta+/-3/2. And, unlike a split-directional scheme,
467! when interpolating along the line eta=+/-1/2 (similarly xi=+/-1/2),
468! it has non-zero contribution from points on the side from the line,
469! except if xi=-1/2; 0; +1/2 (similarly eta=-1/2; 0; +1/2).
470!
471 DO jc=-1,2
472 DO ic=-1,2
473 f(ic,jc)=finp(max(1,min(ubx,i1+ic)), &
474 & max(1,min(uby,j1+jc)))
475 END DO
476 END DO
477
478 f0=c07*(f(1,1)+f(1,0)+f(0,1)+f(0,0))- &
479 & c02*(f(2,0)+f(2,1)+f(1,2)+f(0,2)+ &
480 & f(-1,1)+f(-1,0)+f(0,-1)+f(1,-1))
481
482 fx=c08*(f(1,1)+f(1,0)-f(0,1)-f(0,0))- &
483 & c01*(f(2,1)+f(2,0)-f(-1,1)-f(-1,0))- &
484 & c03*(f(1,2)-f(0,2)+f(1,-1)-f(0,-1))
485
486 fy=c08*(f(1,1)-f(1,0)+f(0,1)-f(0,0))- &
487 & c01*(f(1,2)+f(0,2)-f(1,-1)-f(0,-1))- &
488 & c03*(f(2,1)-f(2,0)+f(-1,1)-f(-1,0))
489
490 fxy=f(1,1)-f(1,0)-f(0,1)+f(0,0)
491
492 fxx=c05*(f(2,1)-f(1,1)-f(0,1)+f(-1,1)+ &
493 & f(2,0)-f(1,0)-f(0,0)+f(-1,0))
494
495 fyy=c05*(f(1,2)-f(1,1)-f(1,0)+f(1,-1)+ &
496 & f(0,2)-f(0,1)-f(0,0)+f(0,-1))
497
498 fxxx=c06*(f(2,1)+f(2,0)-f(-1,1)-f(-1,0))- &
499 & c09*(f(1,1)+f(1,0)-f(0,1)-f(0,0))
500
501 fyyy=c06*(f(1,2)+f(0,2)-f(1,-1)-f(0,-1))- &
502 & c09*(f(1,1)-f(1,0)+f(0,1)-f(0,0))
503
504 fxxy=c06*(f(2,1)-f(1,1)-f(0,1)+f(-1,1)- &
505 & f(2,0)+f(1,0)+f(0,0)-f(-1,0))
506
507 fxyy=c06*(f(1,2)-f(1,1)-f(1,0)+f(1,-1)- &
508 & f(0,2)+f(0,1)+f(0,0)-f(0,-1))
509#else
510!
511! Algorithm below is equivalent to the one above, except that special
512! care is taken to avoid interpolation accross land. This is achieved
513! by shortening the stencil and reducing order of polynomial, if
514! extreme points of the stencil touch land. This is achieved by
515! expressing all f0,fx,fy,...,fxyy in terms of values of interpolated
516! field at the four corners of central cell (which already checked to
517! stay away from land), and eight one-sided differences dfx,dfy (see
518! below) in such a way that field values at the extreme points of the
519! 12-point stencil do not participate directly into f0,fx,...,fxyy.
520! Should an extreme point of the stencil touch land, thus making it
521! impossible to compute the corresponding one-sided difference, this
522! difference is retracted toward the center of the stencil.
523!
524! Optionally, a slope-limiting algorithm may be employed to prevent
525! spurious oscillations of the interpolant. This is a valuable property,
526! if dealing with rough data, however, as a side effect, it turns off
527! high-order interpolation in the vicinity of extrema.
528!
529! The slope-limiting algorithm employed here checks that two consecutive
530! elementary differences, "dfx" and "dfc" have the same sign and differ
531! in magnitude by no more than factor of 3.
532!
533 f(0,0)=finp(i1,j1)
534 f(1,0)=finp(i2,j1)
535 f(0,1)=finp(i1,j2)
536 f(1,1)=finp(i2,j2)
537!
538 dfc=f(1,1)-f(0,1)
539 IF (i1+2.gt.ubx) THEN
540 dfx(1,1)=dfc
541 ELSE IF (finp(i1+2,j2).eq.spv) THEN
542 dfx(1,1)=dfc
543 ELSE
544 dfx(1,1)=finp(i1+2,j2)-f(1,1)
545# ifdef LIMTR
546 IF ((dfx(1,1)*dfc).lt.0.0_r8) THEN
547 dfx(1,1)=0.0_r8
548 ELSE IF (abs(dfx(1,1)).gt.(limtr*abs(dfc))) THEN
549 dfx(1,1)=limtr*dfc
550 END IF
551# endif
552 END IF
553!
554 dfc=f(1,0)-f(0,0)
555 IF ((i1+2).gt.ubx) THEN
556 dfx(1,0)=dfc
557 ELSE IF (finp(i1+2,j1).eq.spv) THEN
558 dfx(1,0)=dfc
559 ELSE
560 dfx(1,0)=finp(i1+2,j1)-f(1,0)
561# ifdef LIMTR
562 IF ((dfx(1,0)*dfc).lt.0.0_r8) THEN
563 dfx(1,0)=0.0_r8
564 ELSE IF (abs(dfx(1,0)).gt.(limtr*abs(dfc))) THEN
565 dfx(1,0)=limtr*dfc
566 END IF
567# endif
568 END IF
569!
570 dfc=f(1,1)-f(0,1)
571 IF (i1-1.lt.1) THEN
572 dfx(0,1)=dfc
573 ELSE IF (finp(i1-1,j2).eq.spv) THEN
574 dfx(0,1)=dfc
575 ELSE
576 dfx(0,1)=f(0,1)-finp(i1-1,j2)
577# ifdef LIMTR
578 IF ((dfx(0,1)*dfc).lt.0.0_r8) THEN
579 dfx(0,1)=0.0_r8
580 ELSE IF (abs(dfx(0,1)).gt.(limtr*abs(dfc))) THEN
581 dfx(0,1)=limtr*dfc
582 END IF
583# endif
584 END IF
585!
586 dfc=f(1,0)-f(0,0)
587 IF (i1-1.lt.1) THEN
588 dfx(0,0)=dfc
589 ELSE IF (finp(i1-1,j1).eq.spv) THEN
590 dfx(0,0)=dfc
591 ELSE
592 dfx(0,0)=f(0,0)-finp(i1-1,j1)
593# ifdef LIMTR
594 IF ((dfx(0,0)*dfc).lt.0.0_r8) THEN
595 dfx(0,0)=0.0_r8
596 ELSE IF (abs(dfx(0,0)).gt.(limtr*abs(dfc))) THEN
597 dfx(0,0)=limtr*dfc
598 END IF
599# endif
600 END IF
601!
602 dfc=f(1,1)-f(1,0)
603 IF (j1+2.gt.uby) THEN
604 dfy(1,1)=dfc
605 ELSE IF (finp(i2,j1+2).eq.spv) THEN
606 dfy(1,1)=dfc
607 ELSE
608 dfy(1,1)=finp(i2,j1+2)-f(1,1)
609# ifdef LIMTR
610 IF ((dfy(1,1)*dfc).lt.0.0_r8) THEN
611 dfy(1,1)=0.0_r8
612 ELSEIF (abs(dfy(1,1)).gt.(limtr*abs(dfc))) THEN
613 dfy(1,1)=limtr*dfc
614 END IF
615# endif
616 END IF
617!
618 dfc=f(0,1)-f(0,0)
619 IF (j1+2.gt.uby) THEN
620 dfy(0,1)=dfc
621 ELSE IF (finp(i1,j1+2).eq.spv) THEN
622 dfy(0,1)=dfc
623 ELSE
624 dfy(0,1)=finp(i1,j1+2)-f(0,1)
625# ifdef LIMTR
626 IF ((dfy(0,1)*dfc).lt.0.0_r8) THEN
627 dfy(0,1)=0.0_r8
628 ELSE IF (abs(dfy(0,1)).gt.(limtr*abs(dfc))) THEN
629 dfy(0,1)=limtr*dfc
630 END IF
631# endif
632 END IF
633!
634 dfc=f(1,1)-f(1,0)
635 IF (j1-1.lt.1) THEN
636 dfy(1,0)=dfc
637 ELSE IF (finp(i2,j1-1).eq.spv) THEN
638 dfy(1,0)=dfc
639 ELSE
640 dfy(1,0)=f(1,0)-finp(i2,j1-1)
641# ifdef LIMTR
642 IF ((dfy(1,0)*dfc).lt.0.0_r8) THEN
643 dfy(1,0)=0.0_r8
644 ELSE IF (abs(dfy(1,0)).gt.(limtr*abs(dfc))) THEN
645 dfy(1,0)=limtr*dfc
646 END IF
647# endif
648 END IF
649!
650 dfc=f(0,1)-f(0,0)
651 IF (j1-1.lt.1) THEN
652 dfy(0,0)=dfc
653 ELSE IF (finp(i1,j1-1).eq.spv) THEN
654 dfy(0,0)=dfc
655 ELSE
656 dfy(0,0)=f(0,0)-finp(i1,j1-1)
657# ifdef LIMTR
658 IF ((dfy(0,0)*dfc).lt.0.0_r8) THEN
659 dfy(0,0)=0.0_r8
660 ELSEIF (abs(dfy(0,0)).gt.(limtr*abs(dfc))) THEN
661 dfy(0,0)=limtr*dfc
662 END IF
663# endif
664 END IF
665!
666 f0=c05*(f(1,1)+f(1,0)+f(0,1)+f(0,0))- &
667 & c02*(dfx(1,1)+dfx(1,0)-dfx(0,1)-dfx(0,0)+ &
668 & dfy(1,1)-dfy(1,0)+dfy(0,1)-dfy(0,0))
669
670 fx=c10*(f(1,1)-f(0,1)+f(1,0)-f(0,0))- &
671 & c01*(dfx(1,1)+dfx(1,0)+dfx(0,1)+dfx(0,0))- &
672 & c03*(dfy(1,1)-dfy(0,1)-dfy(1,0)+dfy(0,0))
673
674 fy=c10*(f(1,1)-f(1,0)+f(0,1)-f(0,0))- &
675 & c01*(dfy(1,1)+dfy(0,1)+dfy(1,0)+dfy(0,0))- &
676 & c03*(dfx(1,1)-dfx(1,0)-dfx(0,1)+dfx(0,0))
677
678 fxy=f(1,1)-f(1,0)-f(0,1)+f(0,0)
679
680 fxx=c05*(dfx(1,1)-dfx(0,1)+dfx(1,0)-dfx(0,0))
681
682 fyy=c05*(dfy(1,1)-dfy(1,0)+dfy(0,1)-dfy(0,0))
683
684 fxxx=c06*(dfx(1,1)+dfx(1,0)+dfx(0,1)+dfx(0,0))- &
685 & f(1,1)+f(0,1)-f(1,0)+f(0,0)
686
687 fyyy=c06*(dfy(1,1)+dfy(0,1)+dfy(1,0)+dfy(0,0))- &
688 & f(1,1)+f(1,0)-f(0,1)+f(0,0)
689
690 fxxy=c06*(dfx(1,1)-dfx(0,1)-dfx(1,0)+dfx(0,0))
691
692 fxyy=c06*(dfy(1,1)-dfy(1,0)-dfy(0,1)+dfy(0,0))
693#endif
694 fout(i,j)=f0+ &
695 & fx*xi+ &
696 & fy*eta+ &
697 & c06*fxx*xi*xi+ &
698 & fxy*xi*eta+ &
699 & c06*fyy*eta*eta+ &
700 & c04*fxxx*xi*xi*xi+ &
701 & c06*fxxy*xi*xi*eta+ &
702 & c04*fyyy*eta*eta*eta+ &
703 & c06*fxyy*xi*eta*eta
704 fmin=min(fmin,fout(i,j))
705 fmax=max(fmax,fout(i,j))
706 END IF
707 END DO
708 END DO
709!
710! Return minimum and maximum values.
711!
712 IF (PRESENT(minval)) THEN
713 minval=fmin
714 END IF
715 IF (PRESENT(maxval)) THEN
716 maxval=fmax
717 END IF
718!
719 RETURN

Referenced by regrid_mod::regrid_nf90(), regrid_mod::regrid_pio(), roms_horiz_interp_2d(), and roms_horiz_interp_3d().

Here is the caller graph for this function:

◆ hindices()

subroutine __interpolate_f__::hindices ( integer, intent(in) ng,
integer, intent(in) lbi,
integer, intent(in) ubi,
integer, intent(in) lbj,
integer, intent(in) ubj,
integer, intent(in) is,
integer, intent(in) ie,
integer, intent(in) js,
integer, intent(in) je,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) angler,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) xgrd,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) ygrd,
integer, intent(in) lbm,
integer, intent(in) ubm,
integer, intent(in) lbn,
integer, intent(in) ubn,
integer, intent(in) ms,
integer, intent(in) me,
integer, intent(in) ns,
integer, intent(in) ne,
real(r8), dimension(lbm:ubm,lbn:ubn), intent(in) xpos,
real(r8), dimension(lbm:ubm,lbn:ubn), intent(in) ypos,
real(r8), dimension(lbm:ubm,lbn:ubn), intent(out) ipos,
real(r8), dimension(lbm:ubm,lbn:ubn), intent(out) jpos,
real(r8), intent(in) ijspv,
logical, intent(in) rectangular )
private

Definition at line 722 of file interpolate.F.

729!
730!=======================================================================
731! !
732! Given any geographical locations Xpos and Ypos, this routine finds !
733! the corresponding array cell indices (Ipos, Jpos) of gridded data !
734! Xgrd and Ygrd containing each requested location. This indices are !
735! usually used for interpolation. !
736! !
737! On Input: !
738! !
739! ng Nested grid number. !
740! LBi I-dimension Lower bound of gridded data. !
741! UBi I-dimension Upper bound of gridded data. !
742! LBj J-dimension Lower bound of gridded data. !
743! UBj J-dimension Upper bound of gridded data. !
744! Is Starting gridded data I-index to search. !
745! Ie Ending gridded data I-index to search. !
746! Js Starting gridded data J-index to search. !
747! Je Ending gridded data J-index to search. !
748! angler Gridded data angle between X-axis and true EAST !
749! (radians). !
750! Xgrd Gridded data X-locations (usually, longitude). !
751! Ygrd Gridded data Y-locations (usually, latitude). !
752! LBm I-dimension Lower bound of requested locations. !
753! UBm I-dimension Upper bound of requested locations. !
754! LBn J-dimension Lower bound of requested locations. !
755! UBn J-dimension Upper bound of requested locations. !
756! Ms Starting requested locations I-index to search. !
757! Me Ending requested locations I-index to search. !
758! Ns Starting requested locations J-index to search. !
759! Ne Ending requested locations J-index to search. !
760! Xpos Requested X-locations to process (usually longitude).!
761! Ypos Requested Y-locations to process (usually latitude). !
762! IJspv Unbounded special value to assign. !
763! rectangular Logical switch indicating that gridded data has a !
764! plaid distribution. !
765! !
766! On Output: !
767! !
768! Ipos Fractional I-cell index containing locations in data. !
769! Jpos Fractional J-cell index containing locations in data. !
770! !
771! Calls: Try_Range !
772! !
773!=======================================================================
774!
775! Imported variable declarations.
776!
777 logical, intent(in) :: rectangular
778!
779 integer, intent(in) :: ng
780 integer, intent(in) :: LBi, UBi, LBj, UBj, Is, Ie, Js, Je
781 integer, intent(in) :: LBm, UBm, LBn, UBn, Ms, Me, Ns, Ne
782
783 real(r8), intent(in) :: IJspv
784!
785 real(r8), intent(in) :: angler(LBi:UBi,LBj:UBj)
786 real(r8), intent(in) :: Xgrd(LBi:UBi,LBj:UBj)
787 real(r8), intent(in) :: Ygrd(LBi:UBi,LBj:UBj)
788
789 real(r8), intent(in) :: Xpos(LBm:UBm,LBn:UBn)
790 real(r8), intent(in) :: Ypos(LBm:UBm,LBn:UBn)
791
792 real(r8), intent(out) :: Ipos(LBm:UBm,LBn:UBn)
793 real(r8), intent(out) :: Jpos(LBm:UBm,LBn:UBn)
794!
795! Local variable declarations.
796!
797 logical :: found, foundi, foundj
798
799 integer :: Imax, Imin, Jmax, Jmin, i, i0, j, j0, mp, np
800
801 real(r8) :: aa2, ang, bb2, diag2, dx, dy, phi
802 real(r8) :: xfac, xpp, yfac, ypp
803!
804!-----------------------------------------------------------------------
805! Determine grid cell indices containing requested position points.
806! Then, interpolate to fractional cell position.
807!-----------------------------------------------------------------------
808!
809 DO np=ns,ne
810 DO mp=ms,me
811 ipos(mp,np)=ijspv
812 jpos(mp,np)=ijspv
813!
814! The gridded data has a plaid distribution so the search is trivial.
815!
816 IF (rectangular) THEN
817 foundi=.false.
818 i_loop : DO i=lbi,ubi-1
819 IF ((xgrd(i ,1).le.xpos(mp,np)).and. &
820 & (xgrd(i+1,1).gt.xpos(mp,np))) THEN
821 imin=i
822 foundi=.true.
823 EXIT i_loop
824 END IF
825 END DO i_loop
826 foundj=.false.
827 j_loop : DO j=lbj,ubj-1
828 IF ((ygrd(1,j ).le.ypos(mp,np)).and. &
829 & (ygrd(1,j+1).gt.ypos(mp,np))) THEN
830 jmin=j
831 foundj=.true.
832 EXIT j_loop
833 END IF
834 END DO j_loop
835 found=foundi.and.foundj
836!
837! Check each position to find if it falls inside the whole domain.
838! Once it is stablished that it inside, find the exact cell to which
839! it belongs by successively dividing the domain by a half (binary
840! search).
841!
842 ELSE
843 found=try_range(ng, lbi, ubi, lbj, ubj, &
844 & xgrd, ygrd, &
845 & is, ie, js, je, &
846 & xpos(mp,np), ypos(mp,np))
847 IF (found) THEN
848 imin=is
849 imax=ie
850 jmin=js
851 jmax=je
852 DO while (((imax-imin).gt.1).or.((jmax-jmin).gt.1))
853 IF ((imax-imin).gt.1) THEN
854 i0=(imin+imax)/2
855 found=try_range(ng, lbi, ubi, lbj, ubj, &
856 & xgrd, ygrd, &
857 & imin, i0, jmin, jmax, &
858 & xpos(mp,np), ypos(mp,np))
859 IF (found) THEN
860 imax=i0
861 ELSE
862 imin=i0
863 END IF
864 END IF
865 IF ((jmax-jmin).gt.1) THEN
866 j0=(jmin+jmax)/2
867 found=try_range(ng, lbi, ubi, lbj, ubj, &
868 & xgrd, ygrd, &
869 & imin, imax, jmin, j0, &
870 & xpos(mp,np), ypos(mp,np))
871 IF (found) THEN
872 jmax=j0
873 ELSE
874 jmin=j0
875 END IF
876 END IF
877 END DO
878 found=(is.le.imin).and.(imin.le.ie).and. &
879 & (is.le.imax).and.(imax.le.ie).and. &
880 & (js.le.jmin).and.(jmin.le.je).and. &
881 & (js.le.jmax).and.(jmax.le.je)
882 END IF
883 END IF
884!
885! Knowing the correct cell, calculate the exact indices, accounting
886! for a possibly rotated grid. If spherical, convert all positions
887! to meters first.
888!
889 IF (found) THEN
890 IF (spherical) THEN
891 yfac=eradius*deg2rad
892 xfac=yfac*cos(ypos(mp,np)*deg2rad)
893 xpp=(xpos(mp,np)-xgrd(imin,jmin))*xfac
894 ypp=(ypos(mp,np)-ygrd(imin,jmin))*yfac
895 ELSE
896 xfac=1.0_r8
897 yfac=1.0_r8
898 xpp=xpos(mp,np)-xgrd(imin,jmin)
899 ypp=ypos(mp,np)-ygrd(imin,jmin)
900 END IF
901!
902! Use Law of Cosines to get cell parallelogram "shear" angle.
903!
904 diag2=((xgrd(imin+1,jmin)-xgrd(imin,jmin+1))*xfac)**2+ &
905 & ((ygrd(imin+1,jmin)-ygrd(imin,jmin+1))*yfac)**2
906 aa2=((xgrd(imin,jmin)-xgrd(imin+1,jmin))*xfac)**2+ &
907 & ((ygrd(imin,jmin)-ygrd(imin+1,jmin))*yfac)**2
908 bb2=((xgrd(imin,jmin)-xgrd(imin,jmin+1))*xfac)**2+ &
909 & ((ygrd(imin,jmin)-ygrd(imin,jmin+1))*yfac)**2
910 phi=asin((diag2-aa2-bb2)/(2.0_r8*sqrt(aa2*bb2)))
911!
912! Transform float position into curvilinear coordinates. Assume the
913! cell is rectanglar, for now.
914!
915 ang=angler(imin,jmin)
916 dx=xpp*cos(ang)+ypp*sin(ang)
917 dy=ypp*cos(ang)-xpp*sin(ang)
918!
919! Correct for parallelogram.
920!
921 dx=dx+dy*tan(phi)
922 dy=dy/cos(phi)
923!
924! Scale with cell side lengths to translate into cell indices.
925!
926 dx=min(max(0.0_r8,dx/sqrt(aa2)),1.0_r8)
927 dy=min(max(0.0_r8,dy/sqrt(bb2)),1.0_r8)
928 ipos(mp,np)=real(imin,r8)+dx
929 jpos(mp,np)=real(jmin,r8)+dy
930 END IF
931 END DO
932 END DO
933!
934 RETURN
try_range
logical function try_range(ng, lbi, ubi, lbj, ubj, xgrd, ygrd, imin, imax, jmin, jmax, xo, yo)
Definition interpolate.F:939

References mod_scalars::deg2rad, mod_scalars::eradius, mod_scalars::spherical, and try_range().

Referenced by grid_coords(), regrid_mod::regrid_nf90(), regrid_mod::regrid_pio(), and roms_interp_fractional().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ inside()

logical function __interpolate_f__::inside ( integer, intent(in) nb,
real(r8), dimension(:), intent(inout) xb,
real(r8), dimension(:), intent(inout) yb,
real(r8), intent(in) xo,
real(r8), intent(in) yo )
private

Definition at line 1013 of file interpolate.F.

1014!
1015!=======================================================================
1016! !
1017! Given the vectors Xb and Yb of size Nb, defining the coordinates !
1018! of a closed polygon, this function find if the point (Xo,Yo) is !
1019! inside the polygon. If the point (Xo,Yo) falls exactly on the !
1020! boundary of the polygon, it still considered inside. !
1021! !
1022! This algorithm does not rely on the setting of Xb(Nb)=Xb(1) and !
1023! Yb(Nb)=Yb(1). Instead, it assumes that the last closing segment !
1024! is (Xb(Nb),Yb(Nb)) --> (Xb(1),Yb(1)). !
1025! !
1026! Reference: !
1027! !
1028! Reid, C., 1969: A long way from Euclid. Oceanography EMR, !
1029! page 174. !
1030! !
1031! Algorithm: !
1032! !
1033! The decision whether the point is inside or outside the polygon !
1034! is done by counting the number of crossings from the ray (Xo,Yo) !
1035! to (Xo,-infinity), hereafter called meridian, by the boundary of !
1036! the polygon. In this counting procedure, a crossing is counted !
1037! as +2 if the crossing happens from "left to right" or -2 if from !
1038! "right to left". If the counting adds up to zero, then the point !
1039! is outside. Otherwise, it is either inside or on the boundary. !
1040! !
1041! This routine is a modified version of the Reid (1969) algorithm, !
1042! where all crossings were counted as positive and the decision is !
1043! made based on whether the number of crossings is even or odd. !
1044! This new algorithm may produce different results in cases where !
1045! Xo accidentally coinsides with one of the (Xb(k),k=1:Nb) points. !
1046! In this case, the crossing is counted here as +1 or -1 depending !
1047! of the sign of (Xb(k+1)-Xb(k)). Crossings are not counted if !
1048! Xo=Xb(k)=Xb(k+1). Therefore, if Xo=Xb(k0) and Yo>Yb(k0), and if !
1049! Xb(k0-1) < Xb(k0) < Xb(k0+1), the crossing is counted twice but !
1050! with weight +1 (for segments with k=k0-1 and k=k0). Similarly if !
1051! Xb(k0-1) > Xb(k0) > Xb(k0+1), the crossing is counted twice with !
1052! weight -1 each time. If, on the other hand, the meridian only !
1053! touches the boundary, that is, for example, Xb(k0-1) < Xb(k0)=Xo !
1054! and Xb(k0+1) < Xb(k0)=Xo, then the crossing is counted as +1 for !
1055! segment k=k0-1 and -1 for segment k=k0, resulting in no crossing. !
1056! !
1057! Note 1: (Explanation of the logical condition) !
1058! !
1059! Suppose that there exist two points (x1,y1)=(Xb(k),Yb(k)) and !
1060! (x2,y2)=(Xb(k+1),Yb(k+1)), such that, either (x1 < Xo < x2) or !
1061! (x1 > Xo > x2). Therefore, meridian x=Xo intersects the segment !
1062! (x1,y1) -> (x2,x2) and the ordinate of the point of intersection !
1063! is: !
1064! !
1065! y1*(x2-Xo) + y2*(Xo-x1) !
1066! y = ----------------------- !
1067! x2-x1 !
1068! !
1069! The mathematical statement that point (Xo,Yo) either coinsides !
1070! with the point of intersection or lies to the north (Yo>=y) from !
1071! it is, therefore, equivalent to the statement: !
1072! !
1073! Yo*(x2-x1) >= y1*(x2-Xo) + y2*(Xo-x1), if x2-x1 > 0 !
1074! or !
1075! Yo*(x2-x1) <= y1*(x2-Xo) + y2*(Xo-x1), if x2-x1 < 0 !
1076! !
1077! which, after noting that Yo*(x2-x1) = Yo*(x2-Xo + Xo-x1) may be !
1078! rewritten as: !
1079! !
1080! (Yo-y1)*(x2-Xo) + (Yo-y2)*(Xo-x1) >= 0, if x2-x1 > 0 !
1081! or !
1082! (Yo-y1)*(x2-Xo) + (Yo-y2)*(Xo-x1) <= 0, if x2-x1 < 0 !
1083! !
1084! and both versions can be merged into essentially the condition !
1085! that (Yo-y1)*(x2-Xo)+(Yo-y2)*(Xo-x1) has the same sign as x2-x1. !
1086! That is, the product of these two must be positive or zero. !
1087! !
1088!=======================================================================
1089!
1090! Imported variable declarations.
1091!
1092 integer, intent(in ) :: Nb
1093
1094 real(r8), intent(in ) :: Xo, Yo
1095
1096#ifdef ASSUMED_SHAPE
1097 real(r8), intent(inout) :: Xb(:), Yb(:)
1098#else
1099 real(r8), intent(inout) :: Xb(Nb+1), Yb(Nb+1)
1100#endif
1101!
1102! Local variable declarations.
1103!
1104 logical :: ItIsInside
1105!
1106 integer, parameter :: Nstep =128
1107
1108 integer :: crossings, i, inc, k, kk, nc
1109
1110 integer, dimension(Nstep) :: Sindex
1111!
1112 real(r8) :: dx1, dx2, dxy
1113!
1114!-----------------------------------------------------------------------
1115! Find intersections.
1116!-----------------------------------------------------------------------
1117!
1118! Set crossings counter and close the contour of the polygon.
1119!
1120 itisinside=.false.
1121 crossings=0
1122 xb(nb+1)=xb(1)
1123 yb(nb+1)=yb(1)
1124!
1125! The search is optimized. First select the indices of segments
1126! where Xb(k) is different from Xb(k+1) and Xo falls between them.
1127! Then, further investigate these segments in a separate loop.
1128! Doing it in two stages takes less time because the first loop is
1129! pipelined.
1130!
1131 DO kk=0,nb-1,nstep
1132 nc=0
1133 DO k=kk+1,min(kk+nstep,nb)
1134 IF (((xb(k+1)-xo)*(xo-xb(k)).ge.0.0_r8).and. &
1135 & (xb(k).ne.xb(k+1))) THEN
1136 nc=nc+1
1137 sindex(nc)=k
1138 END IF
1139 END DO
1140 DO i=1,nc
1141 k=sindex(i)
1142 IF (xb(k).ne.xb(k+1)) THEN
1143 dx1=xo-xb(k)
1144 dx2=xb(k+1)-xo
1145 dxy=dx2*(yo-yb(k))-dx1*(yb(k+1)-yo)
1146 inc=0
1147 IF ((xb(k).eq.xo).and.(yb(k).eq.yo)) THEN
1148 crossings=1
1149 itisinside=.true. ! terminate inner DO-loop
1150 EXIT
1151 ELSE IF (((dx1.eq.0.0_r8).and.(yo.ge.yb(k ))).or. &
1152 & ((dx2.eq.0.0_r8).and.(yo.ge.yb(k+1)))) THEN
1153 inc=1
1154 ELSE IF ((dx1*dx2.gt.0.0_r8).and. &
1155 & ((xb(k+1)-xb(k))*dxy.ge.0.0_r8)) THEN ! see note 1
1156 inc=2
1157 END IF
1158 IF (xb(k+1).gt.xb(k)) THEN
1159 crossings=crossings+inc
1160 ELSE
1161 crossings=crossings-inc
1162 END IF
1163 END IF
1164 END DO
1165 IF (itisinside) EXIT ! terminate outer DO-loop
1166 END DO
1167!
1168! Determine if point (Xo,Yo) is inside of closed polygon.
1169!
1170 IF (crossings.eq.0) THEN
1171 inside=.false.
1172 ELSE
1173 inside=.true.
1174 END IF
1175!
1176 RETURN
inside
logical function inside(nb, xb, yb, xo, yo)
Definition interpolate.F:1014

References inside().

Referenced by inside(), and try_range().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ linterp2d()

subroutine __interpolate_f__::linterp2d ( integer, intent(in) ng,
integer, intent(in) lbx,
integer, intent(in) ubx,
integer, intent(in) lby,
integer, intent(in) uby,
real(r8), dimension(lbx:ubx,lby:uby), intent(in) xinp,
real(r8), dimension(lbx:ubx,lby:uby), intent(in) yinp,
real(r8), dimension(lbx:ubx,lby:uby), intent(in) finp,
integer, intent(in) lbi,
integer, intent(in) ubi,
integer, intent(in) lbj,
integer, intent(in) ubj,
integer, intent(in) istr,
integer, intent(in) iend,
integer, intent(in) jstr,
integer, intent(in) jend,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) iout,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) jout,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) xout,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) yout,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) fout,
real(r8), intent(out), optional minval,
real(r8), intent(out), optional maxval )
private

Definition at line 167 of file interpolate.F.

174!
175!=======================================================================
176! !
177! Given any gridded 2D field, Finp, this routine linearly interpolate !
178! to locations (Xout,Yout). To facilitate the interpolation within !
179! any irregularly gridded 2D field, the fractional grid cell indices !
180! (Iout,Jout) with respect Finp are needed at input. Notice that the !
181! routine "hindices" can be used to compute these indices. !
182! !
183! On Input: !
184! !
185! ng Nested grid number. !
186! LBx I-dimension lower bound of gridded field, Finp. !
187! UBx I-dimension upper bound of gridded field, Finp. !
188! LBy J-dimension lower bound of gridded field, Finp. !
189! UBy J-dimension upper bound of gridded field, Finp. !
190! Xinp X-locations of gridded field, Finp. !
191! Yinp Y-locations of gridded field, Finp. !
192! Finp 2D field to interpolate from. !
193! LBi I-dimension Lower bound of data to interpolate, Fout. !
194! UBi I-dimension Upper bound of data to interpolate, Fout. !
195! LBj J-dimension Lower bound of data to interpolate, Fout. !
196! UBj J-dimension Upper bound of data to interpolate, Fout. !
197! Istr Starting data I-index to interpolate, Fout. !
198! Iend Ending data I-index to interpolate, Fout. !
199! Jstr Starting data J-index to interpolate, Fout. !
200! Jend Ending data J-index to interpolate, Fout. !
201! Iout I-fractional Xinp grid cell containing Xout. !
202! Jout J-fractional Yinp grid cell containing Yout. !
203! Xout X-locations to interpolate, Fout. !
204! Yout Y-locations to interpolate, Fout. !
205! !
206! On Output: !
207! !
208! Fout Interpolated 2D field. !
209! MinVal Interpolated field minimum value (Optional). !
210! MaxVal Interpolated field maximum value (Optional). !
211! !
212!=======================================================================
213!
214! Imported variable declarations.
215!
216 integer, intent(in) :: ng, LBx, UBx, LBy, UBy
217 integer, intent(in) :: LBi, UBi, LBj, UBj
218 integer, intent(in) :: Istr, Iend, Jstr, Jend
219!
220 real(r8), intent(in) :: Xinp(LBx:UBx,LBy:UBy)
221 real(r8), intent(in) :: Yinp(LBx:UBx,LBy:UBy)
222 real(r8), intent(in) :: Finp(LBx:UBx,LBy:UBy)
223
224 real(r8), intent(in) :: Iout(LBi:UBi,LBj:UBj)
225 real(r8), intent(in) :: Jout(LBi:UBi,LBj:UBj)
226 real(r8), intent(in) :: Xout(LBi:UBi,LBj:UBj)
227 real(r8), intent(in) :: Yout(LBi:UBi,LBj:UBj)
228
229 real(r8), intent(out) :: Fout(LBi:UBi,LBj:UBj)
230!
231 real(r8), intent(out), optional :: MinVal
232 real(r8), intent(out), optional :: MaxVal
233!
234! Local variable declarations.
235!
236 integer :: i, i1, i2, j, j1, j2
237!
238 real(r8) :: p1, p2, q1, q2
239 real(r8) :: Fmin, Fmax
240!
241!-----------------------------------------------------------------------
242! Linearly interpolate requested field
243!-----------------------------------------------------------------------
244!
245 fmin=1.0e+35_r8
246 fmax=-1.0e+35_r8
247 DO j=jstr,jend
248 DO i=istr,iend
249 i1=int(iout(i,j))
250 i2=i1+1
251 j1=int(jout(i,j))
252 j2=j1+1
253 IF (((lbx.le.i1).and.(i1.le.ubx)).and. &
254 & ((lby.le.j1).and.(j1.le.uby))) THEN
255 p2=real(i2-i1,r8)*(iout(i,j)-real(i1,r8))
256 q2=real(j2-j1,r8)*(jout(i,j)-real(j1,r8))
257 p1=1.0_r8-p2
258 q1=1.0_r8-q2
259 fout(i,j)=p1*q1*finp(i1,j1)+ &
260 & p2*q1*finp(i2,j1)+ &
261 & p2*q2*finp(i2,j2)+ &
262 & p1*q2*finp(i1,j2)
263 fmin=min(fmin,fout(i,j))
264 fmax=max(fmax,fout(i,j))
265 END IF
266 END DO
267 END DO
268!
269! Return minimum and maximum values.
270!
271 IF (PRESENT(minval)) THEN
272 minval=fmin
273 END IF
274 IF (PRESENT(maxval)) THEN
275 maxval=fmax
276 END IF
277!
278 RETURN

Referenced by regrid_mod::regrid_nf90(), regrid_mod::regrid_pio(), roms_horiz_interp_2d(), and roms_horiz_interp_3d().

Here is the caller graph for this function:

◆ roms_datum_column_interp()

subroutine __interpolate_f__::roms_datum_column_interp ( type (roms_interp_type), intent(inout) s,
real(r8), dimension(:,:,:), intent(in) fld_src,
real(r8), dimension(:,:), intent(out) datum,
integer, intent(in) method )
private

Definition at line 1453 of file interpolate.F.

1454!
1455!=======================================================================
1456! !
1457! Horizontally interpolates a 3D ROMS field to datum locations !
1458! level-by-level. The vertical interpolation will be carried out !
1459! elsewhere from the full column of values. This is what is !
1460! required by the JEDI-UFO operator. !
1461! !
1462! On Input: !
1463! !
1464! S Interpolation structure, TYPE(roms_interp_type), !
1465! initialized elsewhere. !
1466! fld_src 3D field source data (real). !
1467! method Interpolation method flag (integer) !
1468! BilinearMethod => bilinear interpolation !
1469! BicubicMethod => bicubic interpolation !
1470! !
1471! On Output: !
1472! !
1473! datum Interpolated values at datum locations (matrix): !
1474! datum(nlocs,nlevs) !
1475! !
1476!=======================================================================
1477!
1478! Imported variable declarations.
1479!
1480 TYPE (roms_interp_type), intent(inout) :: S
1481!
1482 real(r8), intent( in) :: fld_src(:,:,:)
1483 real(r8), intent(out) :: datum(:,:)
1484!
1485 integer, intent( in) :: method
1486!
1487! Local variable declarations.
1488!
1489 integer :: ic, iobs, i1, i2, j1, j2, k, nlevs, nlocs
1490!
1491 real(r8) :: p1, p2, q1, q2, wsum
1492
1493 real(r8), dimension(4) :: Hmat
1494!
1495 character (len=*), parameter :: MyFile = &
1496 & __FILE__//", roms_datum_interp_2d"
1497!
1498!-----------------------------------------------------------------------
1499! Interpolate from 2D state field to datum locations.
1500!-----------------------------------------------------------------------
1501!
1502 nlevs=SIZE(fld_src, dim=3)
1503 nlocs=SIZE(datum, dim=1)
1504!
1505 DO iobs=1,nlocs
1506 i1=int(s%x_dst(iobs,1))
1507 j1=int(s%y_dst(iobs,1))
1508 IF (((s%Istr_src.le.i1).and.(i1.lt.s%Iend_src)).and. &
1509 & ((s%Jstr_src.le.j1).and.(j1.lt.s%Jend_src))) THEN
1510 i2=i1+1
1511 j2=j1+1
1512 p2=real(i2-i1,r8)*(s%x_dst(iobs,1)-real(i1,r8))
1513 q2=real(j2-j1,r8)*(s%y_dst(iobs,1)-real(j1,r8))
1514 p1=1.0_r8-p2
1515 q1=1.0_r8-q2
1516 hmat(1)=p1*q1
1517 hmat(2)=p2*q1
1518 hmat(3)=p2*q2
1519 hmat(4)=p1*q2
1520#ifdef MASKING
1521 hmat(1)=hmat(1)*s%mask_src(i1,j1)
1522 hmat(2)=hmat(2)*s%mask_src(i2,j1)
1523 hmat(3)=hmat(3)*s%mask_src(i2,j2)
1524 hmat(4)=hmat(4)*s%mask_src(i1,j2)
1525 wsum=0.0_r8
1526 DO ic=1,4
1527 wsum=wsum+hmat(ic)
1528 END DO
1529 IF (wsum.gt.0.0_r8) THEN
1530 wsum=1.0_r8/wsum
1531 DO ic=1,4
1532 hmat(ic)=hmat(ic)*wsum
1533 END DO
1534 END IF
1535#endif
1536 DO k=1,nlevs
1537 datum(iobs,k)=hmat(1)*fld_src(i1,j1,k)+ &
1538 & hmat(2)*fld_src(i2,j1,k)+ &
1539 & hmat(3)*fld_src(i2,j2,k)+ &
1540 & hmat(4)*fld_src(i1,j2,k)
1541 END DO
1542 END IF
1543 END DO
1544!
1545 RETURN

◆ roms_datum_interp_2d()

subroutine __interpolate_f__::roms_datum_interp_2d ( type (roms_interp_type), intent(inout) s,
real(r8), dimension(:,:), intent(in) fld_src,
real(r8), dimension(:), intent(out) datum,
integer, intent(in) method,
real(r8), dimension(:), intent(inout), optional obsvetting )
private

Definition at line 1179 of file interpolate.F.

1181!
1182!=======================================================================
1183! !
1184! Interpolates a ROMS 2D field to observation locations. !
1185! !
1186! On Input: !
1187! !
1188! S Interpolation structure, TYPE(roms_interp_type), !
1189! initialized elsewhere !
1190! fld_src 2D source field data (real) !
1191! method Interpolation method flag (integer) !
1192! BilinearMethod => bilinear interpolation !
1193! BicubicMethod => bicubic interpolation !
1194! !
1195! On Output: !
1196! !
1197! datum Interpolated values at datum locations (vector). !
1198! ObsVetting Observation screenning flag, 0: reject or 1: accept !
1199! !
1200!=======================================================================
1201!
1202! Imported variable declarations.
1203!
1204 integer, intent(in ) :: method
1205!
1206 real(r8), intent(in ) :: fld_src(:,:)
1207 real(r8), intent(out ) :: datum(:)
1208
1209 real(r8), intent(inout), optional :: ObsVetting(:)
1210!
1211 TYPE (roms_interp_type), intent(inout) :: S
1212!
1213! Local variable declarations.
1214!
1215 integer :: ic, iobs, i1, i2, j1, j2, nlocs
1216!
1217 real(r8) :: p1, p2, q1, q2, wsum
1218
1219 real(r8), dimension(4) :: Hmat
1220!
1221 character (len=*), parameter :: MyFile = &
1222 & __FILE__//", roms_datum_interp_2d"
1223!
1224!-----------------------------------------------------------------------
1225! Interpolates a 2D ROMS state field to observation locations.
1226!-----------------------------------------------------------------------
1227!
1228 nlocs=SIZE(datum, dim=1)
1229!
1230 DO iobs=1,nlocs
1231 i1=int(s%x_dst(iobs,1))
1232 j1=int(s%y_dst(iobs,1))
1233 IF (((s%Istr_src.le.i1).and.(i1.lt.s%Iend_src)).and. &
1234 & ((s%Jstr_src.le.j1).and.(j1.lt.s%Jend_src))) THEN
1235 i2=i1+1
1236 j2=j1+1
1237 p2=real(i2-i1,r8)*(s%x_dst(iobs,1)-real(i1,r8))
1238 q2=real(j2-j1,r8)*(s%y_dst(iobs,1)-real(j1,r8))
1239 p1=1.0_r8-p2
1240 q1=1.0_r8-q2
1241 hmat(1)=p1*q1
1242 hmat(2)=p2*q1
1243 hmat(3)=p2*q2
1244 hmat(4)=p1*q2
1245#ifdef MASKING
1246 hmat(1)=hmat(1)*s%mask_src(i1,j1)
1247 hmat(2)=hmat(2)*s%mask_src(i2,j1)
1248 hmat(3)=hmat(3)*s%mask_src(i2,j2)
1249 hmat(4)=hmat(4)*s%mask_src(i1,j2)
1250 wsum=0.0_r8
1251 DO ic=1,4
1252 wsum=wsum+hmat(ic)
1253 END DO
1254 IF (wsum.gt.0.0_r8) THEN
1255 wsum=1.0_r8/wsum
1256 DO ic=1,4
1257 hmat(ic)=hmat(ic)*wsum
1258 END DO
1259 END IF
1260#endif
1261 datum(iobs)=hmat(1)*fld_src(i1,j1)+ &
1262 & hmat(2)*fld_src(i2,j1)+ &
1263 & hmat(3)*fld_src(i2,j2)+ &
1264 & hmat(4)*fld_src(i1,j2)
1265!
1266! Set observation screening flag, 0: reject or 1: accept.
1267!
1268 IF (PRESENT(obsvetting)) THEN
1269#ifdef MASKING
1270 IF (wsum.gt.0.0_r8) obsvetting(iobs)=1.0_r8
1271#else
1272 obsvetting(iobs)=1.0_r8
1273#endif
1274 END IF
1275 END IF
1276 END DO
1277!
1278 RETURN

◆ roms_datum_interp_3d()

subroutine __interpolate_f__::roms_datum_interp_3d ( type (roms_interp_type), intent(inout) s,
real(r8), dimension(:,:,:), intent(in) fld_src,
real(r8), dimension(:,:,:), intent(in) z_src,
real(r8), dimension(:), intent(in) z_locs,
real(r8), dimension(:), intent(out) datum,
integer, intent(in) method,
real(r8), dimension(:), intent(inout), optional obsvetting )
private

Definition at line 1281 of file interpolate.F.

1283!
1284!=======================================================================
1285! !
1286! Interpolates a 3D ROMS field to observation horizontal and vertical !
1287! locations. !
1288! !
1289! On Input: !
1290! !
1291! S Interpolation structure, TYPE(roms_interp_type), !
1292! initialized elsewhere !
1293! fld_src 3D source field data (real) !
1294! z_src 3D source field depth (real; meter; negative) !
1295! z_locs Datum depth locations (real vector). Negative values !
1296! are depths in meters while positive values are in !
1297! vertical fractional coordinates of 1 to N. !
1298! method Interpolation method flag (integer) !
1299! BilinearMethod => bilinear interpolation !
1300! BicubicMethod => bicubic interpolation !
1301! !
1302! On Output: !
1303! !
1304! datum Interpolated values at observation locations (vector) !
1305! ObsVetting Observation screenning flag, 0: reject or 1: accept !
1306! !
1307!=======================================================================
1308!
1309! Imported variable declarations.
1310!
1311 integer, intent( in) :: method
1312!
1313 real(r8), intent( in) :: fld_src(:,:,:)
1314 real(r8), intent( in) :: z_src(:,:,:)
1315 real(r8), intent( in) :: z_locs(:)
1316 real(r8), intent(out) :: datum(:)
1317
1318 real(r8), intent(inout), optional :: ObsVetting(:)
1319!
1320 TYPE (roms_interp_type), intent(inout) :: S
1321!
1322! Local variable declarations.
1323!
1324 integer :: ic, iobs, i1, i2, j1, j2, k, k1, k2, nlevs, nlocs
1325!
1326 real(r8) :: Zbot, Ztop, dz, p1, p2, q1, q2, r1, r2
1327 real(r8) :: w11, w12, w21, w22, wsum
1328
1329 real(r8), dimension(8) :: Hmat
1330!
1331 character (len=*), parameter :: MyFile = &
1332 & __FILE__//", roms_datum_interp_3d"
1333!
1334!-----------------------------------------------------------------------
1335! Interpolate from 3D ROMS state field at observation locations.
1336!-----------------------------------------------------------------------
1337!
1338 nlevs=SIZE(fld_src, dim=3)
1339 nlocs=SIZE(datum, dim=1)
1340!
1341 DO iobs=1,nlocs
1342 i1=int(s%x_dst(iobs,1))
1343 j1=int(s%y_dst(iobs,1))
1344 IF (((s%Istr_src.le.i1).and.(i1.lt.s%Iend_src)).and. &
1345 & ((s%Jstr_src.le.j1).and.(j1.lt.s%Jend_src))) THEN
1346 i2=i1+1
1347 j2=j1+1
1348 p2=real(i2-i1,r8)*(s%x_dst(iobs,1)-real(i1,r8))
1349 q2=real(j2-j1,r8)*(s%y_dst(iobs,1)-real(j1,r8))
1350 p1=1.0_r8-p2
1351 q1=1.0_r8-q2
1352 w11=p1*q1
1353 w21=p2*q1
1354 w22=p2*q2
1355 w12=p1*q2
1356 IF (z_locs(iobs).gt.0.0_r8) THEN
1357 k1=max(1,int(z_locs(iobs))) ! Positions in fractional
1358 k2=min(int(z_locs(iobs))+1,nlevs) ! levels
1359 r2=real(k2-k1,r8)*(z_locs(iobs)-real(k1,r8))
1360 r1=1.0_r8-r2
1361 ELSE
1362 ztop=z_src(i1,j1,nlevs)
1363 zbot=z_src(i1,j1,1 )
1364 IF (z_locs(iobs).ge.ztop) THEN
1365 r1=0.0_r8 ! If shallower, ignore.
1366 r2=0.0_r8
1367 IF (PRESENT(obsvetting)) obsvetting(iobs)=0.0_r8
1368 ELSE IF (zbot.ge.z_locs(iobs)) THEN
1369 r1=0.0_r8 ! If deeper, ignore.
1370 r2=0.0_r8
1371 IF (PRESENT(obsvetting)) obsvetting(iobs)=0.0_r8
1372 ELSE
1373 DO k=nlevs,2,-1 ! Otherwise, interpolate
1374 ztop=z_src(i1,j1,k ) ! to fractional level
1375 zbot=z_src(i1,j1,k-1)
1376 IF ((ztop.gt.z_locs(iobs)).and. &
1377 (z_locs(iobs).ge.zbot)) THEN
1378 k1=k-1
1379 k2=k
1380 END IF
1381 END DO
1382 dz=z_src(i1,j1,k2)-z_src(i1,j1,k1)
1383 r2=(z_locs(iobs)-z_src(i1,j1,k1))/dz
1384 r1=1.0_r8-r2
1385 END IF
1386 END IF
1387 IF ((r1+r2).gt.0.0_r8) THEN
1388 hmat(1)=w11*r1
1389 hmat(2)=w21*r1
1390 hmat(3)=w22*r1
1391 hmat(4)=w12*r1
1392 hmat(5)=w11*r2
1393 hmat(6)=w21*r2
1394 hmat(7)=w22*r2
1395 hmat(8)=w12*r2
1396#ifdef MASKING
1397 hmat(1)=hmat(1)*s%mask_src(i1,j1)
1398 hmat(2)=hmat(2)*s%mask_src(i2,j1)
1399 hmat(3)=hmat(3)*s%mask_src(i2,j2)
1400 hmat(4)=hmat(4)*s%mask_src(i1,j2)
1401 hmat(5)=hmat(5)*s%mask_src(i1,j1)
1402 hmat(6)=hmat(6)*s%mask_src(i2,j1)
1403 hmat(7)=hmat(7)*s%mask_src(i2,j2)
1404 hmat(8)=hmat(8)*s%mask_src(i1,j2)
1405 wsum=0.0_r8
1406 DO ic=1,8
1407 wsum=wsum+hmat(ic)
1408 END DO
1409 IF (wsum.gt.0.0_r8) THEN
1410 wsum=1.0_r8/wsum
1411 DO ic=1,8
1412 hmat(ic)=hmat(ic)*wsum
1413 END DO
1414 END IF
1415#endif
1416 datum(iobs)=hmat(1)*fld_src(i1,j1,k1)+ &
1417 & hmat(2)*fld_src(i2,j1,k1)+ &
1418 & hmat(3)*fld_src(i2,j2,k1)+ &
1419 & hmat(4)*fld_src(i1,j2,k1)+ &
1420 & hmat(5)*fld_src(i1,j1,k2)+ &
1421 & hmat(6)*fld_src(i2,j1,k2)+ &
1422 & hmat(7)*fld_src(i2,j2,k2)+ &
1423 & hmat(8)*fld_src(i1,j2,k2)
1424!
1425! Set observation screening flag, 0: reject or 1: accept.
1426!
1427 IF (PRESENT(obsvetting)) THEN
1428#ifdef MASKING
1429 IF (wsum.gt.0.0_r8) obsvetting(iobs)=1.0_r8
1430#else
1431 obsvetting(iobs)=1.0_r8
1432#endif
1433#ifndef ALLOW_BOTTOM_OBS
1434!
1435! Reject observations that lie in the lower bottom grid cell (k=1) to
1436! avoid clustering due shallowing of bathymetry during smoothing and
1437! coarse level half-thickness (-h < z_locs < z_src(:,:,1)) in deep
1438! water.
1439!
1440 IF ((z_locs(iobs).gt.0.0_r8).and. &
1441 (z_locs(iobs).le.1.0_r8)) THEN
1442 obsvetting(iobs)=0.0_r8
1443 END IF
1444#endif
1445 END IF
1446 END IF
1447 END IF
1448 END DO
1449!
1450 RETURN

◆ roms_horiz_interp_2d()

subroutine __interpolate_f__::roms_horiz_interp_2d ( type (roms_interp_type), intent(inout) s,
real(r8), dimension(:,:), intent(in) fld_src,
real(r8), dimension(:,:), intent(out) fld_dst,
integer, intent(in) method )
private

Definition at line 1548 of file interpolate.F.

1549!
1550!=======================================================================
1551! !
1552! Interpolate a ROMS 2D source field to a new destination locations. !
1553! !
1554! On Input: !
1555! !
1556! S Interpolation structure, TYPE(roms_interp_type), !
1557! initialized elsewhere. !
1558! fld_src 2D field source data (real). !
1559! method Interpolation method flag (integer) !
1560! BilinearMethod => bilinear interpolation !
1561! BicubicMethod => bicubic interpolation !
1562! !
1563! On Output: !
1564! !
1565! fld_dst 2D interpolated field data (real). !
1566! !
1567!=======================================================================
1568!
1569! Imported variable declarations.
1570!
1571 TYPE (roms_interp_type), intent(inout) :: S
1572!
1573 real(r8), intent(in ) :: fld_src(:,:)
1574 real(r8), intent(out) :: fld_dst(:,:)
1575!
1576 integer, intent(in ) :: method
1577!
1578! Local variable declarations.
1579!
1580 character (len=*), parameter :: MyFile = &
1581 & __FILE__//", roms_horiz_interp_2d"
1582!
1583!-----------------------------------------------------------------------
1584! Interpolate 2D field.
1585!-----------------------------------------------------------------------
1586!
1587 SELECT CASE (method)
1588 CASE (bilinearmethod)
1589 CALL linterp2d (s%ng, &
1590 & s%LBi_src, s%UBi_src, s%LBj_src, s%UBj_src, &
1591 & s%lon_src, s%lat_src, fld_src, &
1592 & s%LBi_dst, s%UBi_dst, s%LBj_dst, s%UBj_dst, &
1593 & s%Istr_dst, s%Iend_dst, &
1594 & s%Jstr_dst, s%Jend_dst, &
1595 & s%x_dst, s%y_dst, s%lon_dst, s%lat_dst, &
1596 & fld_dst, &
1597 & minval = s%min_dst, &
1598 & maxval = s%max_dst)
1599 CASE (bicubicmethod)
1600 CALL cinterp2d (s%ng, &
1601 & s%LBi_src, s%UBi_src, s%LBj_src, s%UBj_src, &
1602 & s%lon_src, s%lat_src, fld_src, &
1603 & s%LBi_dst, s%UBi_dst, s%LBj_dst, s%UBj_dst, &
1604 & s%Istr_dst, s%Iend_dst, &
1605 & s%Jstr_dst, s%Jend_dst, &
1606 & s%x_dst, s%y_dst, s%lon_dst, s%lat_dst, &
1607 & fld_dst, &
1608 & minval = s%min_dst, &
1609 & maxval = s%max_dst)
1610 CASE DEFAULT
1611 IF (master) WRITE(stdout,10) method
1612 exit_flag=5
1613 END SELECT
1614 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1615!
1616 10 FORMAT (' ROMS_HORIZ_INTERP_2D - Illegal interpolation method =', &
1617 & i0)
1618!
1619 RETURN
cinterp2d
subroutine cinterp2d(ng, lbx, ubx, lby, uby, xinp, yinp, finp, lbi, ubi, lbj, ubj, istr, iend, jstr, jend, iout, jout, xout, yout, fout, minval, maxval)
Definition interpolate.F:288
linterp2d
subroutine linterp2d(ng, lbx, ubx, lby, uby, xinp, yinp, finp, lbi, ubi, lbj, ubj, istr, iend, jstr, jend, iout, jout, xout, yout, fout, minval, maxval)
Definition interpolate.F:174

References cinterp2d(), mod_scalars::exit_flag, strings_mod::founderror(), linterp2d(), mod_parallel::master, mod_scalars::noerror, and mod_iounits::stdout.

Here is the call graph for this function:

◆ roms_horiz_interp_3d()

subroutine __interpolate_f__::roms_horiz_interp_3d ( type (roms_interp_type), intent(inout) s,
real(r8), dimension(:,:,:), intent(in) fld_src,
real(r8), dimension(:,:,:), intent(out) fld_dst,
integer, intent(in) method )
private

Definition at line 1622 of file interpolate.F.

1623!
1624!=======================================================================
1625! !
1626! Interpolate a ROMS 3D source field to a new destination locations !
1627! level-by-level. There is no vertical interpolation. !
1628! !
1629! On Input: !
1630! !
1631! S Interpolation structure, TYPE(roms_interp_type), !
1632! initialized elsewhere. !
1633! fld_src 3D field source data (real). !
1634! method Interpolation method flag (integer) !
1635! BilinearMethod => bilinear interpolation !
1636! BicubicMethod => bicubic interpolation !
1637! !
1638! On Output: !
1639! !
1640! fld_dst 3D interpolated field data (real). !
1641! !
1642!=======================================================================
1643!
1644! Imported variable declarations.
1645!
1646 TYPE (roms_interp_type), intent(inout) :: S
1647!
1648 real(r8), intent( in) :: fld_src(:,:,:)
1649 real(r8), intent(out) :: fld_dst(:,:,:)
1650!
1651 integer, intent( in) :: method
1652!
1653! Local variable declarations.
1654!
1655 integer :: k, nk
1656!
1657 character (len=*), parameter :: MyFile = &
1658 & __FILE__//", roms_horiz_interp_3d"
1659!
1660!-----------------------------------------------------------------------
1661! Interpolate 3D field level-by-level.
1662!-----------------------------------------------------------------------
1663!
1664 nk=SIZE(fld_src, dim=3)!
1665!
1666 SELECT CASE (method)
1667 CASE (bilinearmethod)
1668 DO k=1,nk
1669 CALL linterp2d (s%ng, &
1670 & s%LBi_src, s%UBi_src, s%LBj_src, s%UBj_src, &
1671 & s%lon_src, s%lat_src, fld_src(:,:,k), &
1672 & s%LBi_dst, s%UBi_dst, s%LBj_dst, s%UBj_dst, &
1673 & s%Istr_dst, s%Iend_dst, &
1674 & s%Jstr_dst, s%Jend_dst, &
1675 & s%x_dst, s%y_dst, s%lon_dst, s%lat_dst, &
1676 & fld_dst(:,:,k), &
1677 & minval = s%min_dst, &
1678 & maxval = s%max_dst)
1679 END DO
1680 CASE (bicubicmethod)
1681 DO k=1,nk
1682 CALL cinterp2d (s%ng, &
1683 & s%LBi_src, s%UBi_src, s%LBj_src, s%UBj_src, &
1684 & s%lon_src, s%lat_src, fld_src(:,:,k), &
1685 & s%LBi_dst, s%UBi_dst, s%LBj_dst, s%UBj_dst, &
1686 & s%Istr_dst, s%Iend_dst, &
1687 & s%Jstr_dst, s%Jend_dst, &
1688 & s%x_dst, s%y_dst, s%lon_dst, s%lat_dst, &
1689 & fld_dst(:,:,k), &
1690 & minval = s%min_dst, &
1691 & maxval = s%max_dst)
1692 END DO
1693 CASE DEFAULT
1694 IF (master) WRITE(stdout,10) method
1695 exit_flag=5
1696 END SELECT
1697 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1698!
1699 10 FORMAT (' ROMS_HORIZ_INTERP_3D - Illegal interpolation method =', &
1700 & i0)
1701!
1702 RETURN

References cinterp2d(), mod_scalars::exit_flag, strings_mod::founderror(), linterp2d(), mod_parallel::master, mod_scalars::noerror, and mod_iounits::stdout.

Here is the call graph for this function:

◆ roms_interp_create()

subroutine __interpolate_f__::roms_interp_create ( integer, intent(in) ng,
integer, intent(in) tile,
integer, intent(in) gtype,
type (roms_interp_type), intent(inout) s )
private

Definition at line 1705 of file interpolate.F.

1706!
1707!=======================================================================
1708! !
1709! Creates ROMS intepolation Object. The interpolation structure has !
1710! two components: source and destination. Here, the source component !
1711! is initialized. The destination component is initialized elsewhere !
1712! for flexibility. !
1713! !
1714!=======================================================================
1715!
1716! Imported variable declarations.
1717!
1718 integer, intent(in) :: ng, tile, gtype
1719!
1720 TYPE (roms_interp_type), intent(inout) :: S
1721!
1722! Local variable declarations.
1723!
1724 integer :: LBi, UBi, LBj, UBj
1725!
1726!-----------------------------------------------------------------------
1727! Creates source component of interpolation object.
1728!-----------------------------------------------------------------------
1729!
1730! If applicable, destroy interpolation object.
1731!
1732 IF (allocated(s%lon_src)) THEN
1733 CALL roms_interp_delete (s)
1734 END IF
1735!
1736! Allocate source component arrays.
1737!
1738 lbi=bounds(ng)%LBi(tile)
1739 ubi=bounds(ng)%UBi(tile)
1740 lbj=bounds(ng)%LBj(tile)
1741 ubj=bounds(ng)%UBj(tile)
1742!
1743 allocate ( s%lon_src(lbi:ubi,lbj:ubj) )
1744 allocate ( s%lat_src(lbi:ubi,lbj:ubj) )
1745 allocate ( s%angle_src(lbi:ubi,lbj:ubj) )
1746#ifdef MASKING
1747 allocate ( s%mask_src(lbi:ubi,lbj:ubj) )
1748#endif
1749!
1750! Associate and assign values to source component.
1751!
1752 s%ng = ng
1753 s%model = inlm
1754 s%spval = 0.0_r8
1755!
1756 SELECT CASE (gtype)
1757 CASE (r2dvar, r3dvar)
1758 s%lon_src = grid(ng)%lonr
1759 s%lat_src = grid(ng)%latr
1760 s%angle_src = grid(ng)%angler
1761#ifdef MASKING
1762 s%mask_src = grid(ng)%rmask
1763#endif
1764 s%Istr_src = bounds(ng)%IstrR(tile)
1765 s%Iend_src = bounds(ng)%IendR(tile)
1766 s%Jstr_src = bounds(ng)%JstrR(tile)
1767 s%Jend_src = bounds(ng)%JendR(tile)
1768 CASE (u2dvar, u3dvar)
1769 s%lon_src = grid(ng)%lonu
1770 s%lat_src = grid(ng)%latu
1771 s%angle_src = grid(ng)%angler
1772#ifdef MASKING
1773 s%mask_src = grid(ng)%umask
1774#endif
1775 s%Istr_src = bounds(ng)%Istr (tile)
1776 s%Iend_src = bounds(ng)%IendR(tile)
1777 s%Jstr_src = bounds(ng)%JstrR(tile)
1778 s%Jend_src = bounds(ng)%JendR(tile)
1779 CASE (v2dvar, v3dvar)
1780 s%lon_src = grid(ng)%lonv
1781 s%lat_src = grid(ng)%latv
1782 s%angle_src = grid(ng)%angler
1783#ifdef MASKING
1784 s%mask_src = grid(ng)%vmask
1785#endif
1786 s%Istr_src = bounds(ng)%IstrR(tile)
1787 s%Iend_src = bounds(ng)%IendR(tile)
1788 s%Jstr_src = bounds(ng)%Jstr (tile)
1789 s%Jend_src = bounds(ng)%JendR(tile)
1790 END SELECT
1791!
1792!
1793 s%LBi_src = lbi
1794 s%UBi_src = ubi
1795 s%LBj_src = lbj
1796 s%UBj_src = ubj
1797!
1798 RETURN
roms_interp_delete
subroutine roms_interp_delete(s)
Definition interpolate.F:1802

References mod_param::bounds, mod_grid::grid, mod_param::inlm, mod_param::r2dvar, mod_param::r3dvar, roms_interp_delete(), mod_param::u2dvar, mod_param::u3dvar, mod_param::v2dvar, and mod_param::v3dvar.

Here is the call graph for this function:

◆ roms_interp_delete()

subroutine __interpolate_f__::roms_interp_delete ( type (roms_interp_type), intent(inout) s)
private

Definition at line 1801 of file interpolate.F.

1802!
1803!=======================================================================
1804! !
1805! Deallocates pointer arrays in the intepolation structure, S. !
1806! !
1807!=======================================================================
1808!
1809! Imported variable declarations.
1810!
1811 TYPE (roms_interp_type), intent(inout) :: S
1812!
1813!-----------------------------------------------------------------------
1814! Deallocate interpolation structure.
1815!-----------------------------------------------------------------------
1816!
1817 IF (allocated(s%angle_src)) deallocate (s%angle_src)
1818 IF (allocated(s%lon_src)) deallocate (s%lon_src)
1819 IF (allocated(s%lat_src)) deallocate (s%lat_src)
1820 IF (allocated(s%mask_src)) deallocate (s%mask_src)
1821
1822 IF (allocated(s%lon_dst)) deallocate (s%lon_dst)
1823 IF (allocated(s%lat_dst)) deallocate (s%lat_dst)
1824 IF (allocated(s%mask_dst)) deallocate (s%mask_dst)
1825 IF (allocated(s%x_dst)) deallocate (s%x_dst)
1826 IF (allocated(s%y_dst)) deallocate (s%y_dst)
1827!
1828 RETURN

Referenced by roms_interp_create().

Here is the caller graph for this function:

◆ roms_interp_fractional()

subroutine __interpolate_f__::roms_interp_fractional ( type (roms_interp_type), intent(inout) s)
private

Definition at line 1831 of file interpolate.F.

1832!
1833!=======================================================================
1834! !
1835! Computes the horizontal fractional coordinates (S%x_dst, S%y_dst) !
1836! of the source cells containing the destination values. !
1837! !
1838!=======================================================================
1839!
1840! Imported variable declarations.
1841!
1842 TYPE (roms_interp_type), intent(inout) :: S
1843!
1844! Local variable declarations.
1845!
1846 integer :: Npts
1847!
1848!-----------------------------------------------------------------------
1849! Determine the horizontal fractional coordinates: S%x_dst and S%y_dst.
1850!-----------------------------------------------------------------------
1851!
1852 CALL hindices (s%ng, &
1853 & s%LBi_src, s%UBi_src, s%LBj_src, s%UBj_src, &
1854 & s%Istr_src , s%Iend_src, s%Jstr_src, s%Jend_src, &
1855 & s%angle_src, s%lon_src, s%lat_src, &
1856 & s%LBi_dst, s%UBi_dst, s%LBj_dst, s%UBj_dst, &
1857 & s%Istr_dst, s%Iend_dst, s%Jstr_dst, s%Jend_dst, &
1858 & s%lon_dst, s%lat_dst, s%x_dst, s%y_dst, &
1859 & s%spval, s%rectangular)
1860
1861#ifdef DISTRIBUTE
1862!
1863!-----------------------------------------------------------------------
1864! Collect fractional coordinates.
1865!-----------------------------------------------------------------------
1866!
1867 npts=SIZE(s%x_dst)
1868 CALL mp_assemble (s%ng, s%model, npts, s%spval, s%x_dst)
1869 CALL mp_assemble (s%ng, s%model, npts, s%spval, s%y_dst)
1870#endif
1871!
1872 RETURN
hindices
subroutine hindices(ng, lbi, ubi, lbj, ubj, is, ie, js, je, angler, xgrd, ygrd, lbm, ubm, lbn, ubn, ms, me, ns, ne, xpos, ypos, ipos, jpos, ijspv, rectangular)
Definition interpolate.F:729

References hindices().

Here is the call graph for this function:

◆ try_range()

logical function __interpolate_f__::try_range ( integer, intent(in) ng,
integer, intent(in) lbi,
integer, intent(in) ubi,
integer, intent(in) lbj,
integer, intent(in) ubj,
real(r8), dimension(lbi:,lbj:), intent(in) xgrd,
real(r8), dimension(lbi:,lbj:), intent(in) ygrd,
integer, intent(in) imin,
integer, intent(in) imax,
integer, intent(in) jmin,
integer, intent(in) jmax,
real(r8), intent(in) xo,
real(r8), intent(in) yo )
private

Definition at line 937 of file interpolate.F.

939!
940!=======================================================================
941! !
942! Given a grided domain with matrix coordinates Xgrd and Ygrd, this !
943! function finds if the point (Xo,Yo) is inside the box defined by !
944! the requested corners (Imin,Jmin) and (Imax,Jmax). It will return !
945! logical switch try_range=.TRUE. if (Xo,Yo) is inside, otherwise !
946! it will return false. !
947! !
948! Calls: inside !
949! !
950!=======================================================================
951!
952! Imported variable declarations.
953!
954 integer, intent(in) :: ng, LBi, UBi, LBj, UBj
955 integer, intent(in) :: Imin, Imax, Jmin, Jmax
956
957#ifdef ASSUMED_SHAPE
958 real(r8), intent(in) :: Xgrd(LBi:,LBj:)
959 real(r8), intent(in) :: Ygrd(LBi:,LBj:)
960#else
961 real(r8), intent(in) :: Xgrd(LBi:UBi,LBj:UBj)
962 real(r8), intent(in) :: Ygrd(LBi:UBi,LBj:UBj)
963#endif
964
965 real(r8), intent(in) :: Xo, Yo
966!
967! Local variable declarations.
968!
969 integer :: Nb, i, j, shft, ic
970
971 real(r8), dimension(2*(Jmax-Jmin+Imax-Imin)+1) :: Xb, Yb
972!
973!-----------------------------------------------------------------------
974! Define closed polygon.
975!-----------------------------------------------------------------------
976!
977! Note that the last point (Xb(Nb),Yb(Nb)) does not repeat first
978! point (Xb(1),Yb(1)). Instead, in function inside, it is implied
979! that the closing segment is (Xb(Nb),Yb(Nb))-->(Xb(1),Yb(1)). In
980! fact, function inside sets Xb(Nb+1)=Xb(1) and Yb(Nb+1)=Yb(1).
981!
982 nb=2*(jmax-jmin+imax-imin)
983 shft=1-imin
984 DO i=imin,imax-1
985 xb(i+shft)=xgrd(i,jmin)
986 yb(i+shft)=ygrd(i,jmin)
987 END DO
988 shft=1-jmin+imax-imin
989 DO j=jmin,jmax-1
990 xb(j+shft)=xgrd(imax,j)
991 yb(j+shft)=ygrd(imax,j)
992 END DO
993 shft=1+jmax-jmin+2*imax-imin
994 DO i=imax,imin+1,-1
995 xb(shft-i)=xgrd(i,jmax)
996 yb(shft-i)=ygrd(i,jmax)
997 END DO
998 shft=1+2*jmax-jmin+2*(imax-imin)
999 DO j=jmax,jmin+1,-1
1000 xb(shft-j)=xgrd(imin,j)
1001 yb(shft-j)=ygrd(imin,j)
1002 END DO
1003!
1004!-----------------------------------------------------------------------
1005! Check if point (Xo,Yo) is inside of the defined polygon.
1006!-----------------------------------------------------------------------
1007!
1008 try_range=inside(nb, xb, yb, xo, yo)
1009!
1010 RETURN

References inside(), and try_range().

Referenced by hindices(), and try_range().

Here is the call graph for this function:
Here is the caller graph for this function: