| 1 | subroutine bilinearbig(nX,nY,x_arr,y_arr,f2d_arr,x_in,y_in,f,ind) |
|---|
| 2 | |
|---|
| 3 | ! Necessary for interpolation of continuum data |
|---|
| 4 | ! optimized by A. Spiga 01/2013 |
|---|
| 5 | |
|---|
| 6 | implicit none |
|---|
| 7 | |
|---|
| 8 | integer nX,nY,i,j,ind,b |
|---|
| 9 | |
|---|
| 10 | real*8 x_in,y_in,x1,x2,y1,y2 |
|---|
| 11 | real*8 f,f11,f12,f21,f22,fA,fB |
|---|
| 12 | real*8 x_arr(nX) |
|---|
| 13 | real*8 y_arr(nY) |
|---|
| 14 | real*8 f2d_arr(nX,nY) |
|---|
| 15 | real*8,save :: x,y |
|---|
| 16 | !$OMP THREADPRIVATE(x,y) |
|---|
| 17 | |
|---|
| 18 | integer strlen |
|---|
| 19 | character*100 label |
|---|
| 20 | label='subroutine bilinear' |
|---|
| 21 | |
|---|
| 22 | |
|---|
| 23 | x=x_in |
|---|
| 24 | y=y_in |
|---|
| 25 | |
|---|
| 26 | !! AS: important to optimize here because the array is quite large |
|---|
| 27 | !! ... and actually calculations only need to be done once |
|---|
| 28 | !! IF ind=-9999 we have not calculated yet |
|---|
| 29 | if ( ind == -9999) then |
|---|
| 30 | !1st check we're within the wavenumber range |
|---|
| 31 | if ((x.lt.x_arr(2)).or.(x.gt.x_arr(nX-2))) then |
|---|
| 32 | ind=-1 |
|---|
| 33 | else |
|---|
| 34 | i=1 |
|---|
| 35 | x2=x_arr(i) |
|---|
| 36 | do while ( x2 .le. x ) |
|---|
| 37 | x1=x2 |
|---|
| 38 | i=i+1 |
|---|
| 39 | x2=x_arr(i) |
|---|
| 40 | ind=i-1 |
|---|
| 41 | end do |
|---|
| 42 | endif |
|---|
| 43 | endif |
|---|
| 44 | |
|---|
| 45 | !! Either we already saw we are out of wavenumber range |
|---|
| 46 | !! ... and we just have to set f=0 and exit |
|---|
| 47 | if ( ind == -1) then |
|---|
| 48 | f=0.0D+0 |
|---|
| 49 | return |
|---|
| 50 | !! Or we already determined ind -- so we just proceed |
|---|
| 51 | else |
|---|
| 52 | x1=x_arr(ind) |
|---|
| 53 | x2=x_arr(ind+1) |
|---|
| 54 | endif |
|---|
| 55 | |
|---|
| 56 | ! ... and for y within the temperature range |
|---|
| 57 | if ((y.le.y_arr(1)).or.(y.ge.y_arr(nY))) then |
|---|
| 58 | !print*,y_arr(1),y_arr(nY) |
|---|
| 59 | !write(*,*) 'Warning from bilinearbig routine:' |
|---|
| 60 | !write(*,*) 'Outside continuum temperature range!' |
|---|
| 61 | if(y.le.y_arr(1))then |
|---|
| 62 | y=y_arr(1)+0.01 |
|---|
| 63 | b=1 |
|---|
| 64 | y1=y_arr(b) |
|---|
| 65 | y2=y_arr(b+1) |
|---|
| 66 | endif |
|---|
| 67 | if(y.ge.y_arr(nY))then |
|---|
| 68 | y=y_arr(nY)-0.01 |
|---|
| 69 | b=nY-1 |
|---|
| 70 | y1=y_arr(b) |
|---|
| 71 | y2=y_arr(b+1) |
|---|
| 72 | endif |
|---|
| 73 | else |
|---|
| 74 | j=1 |
|---|
| 75 | y2=y_arr(j) |
|---|
| 76 | do while ( y2 .lt. y ) |
|---|
| 77 | y1=y2 |
|---|
| 78 | j=j+1 |
|---|
| 79 | y2=y_arr(j) |
|---|
| 80 | b=j-1 |
|---|
| 81 | end do |
|---|
| 82 | endif |
|---|
| 83 | |
|---|
| 84 | f11=f2d_arr(ind,b) |
|---|
| 85 | f21=f2d_arr(ind+1,b) |
|---|
| 86 | f12=f2d_arr(ind,b+1) |
|---|
| 87 | f22=f2d_arr(ind+1,b+1) |
|---|
| 88 | |
|---|
| 89 | call bilinear(f,f11,f21,f12,f22,x,x1,x2,y,y1,y2) |
|---|
| 90 | |
|---|
| 91 | return |
|---|
| 92 | end subroutine bilinearbig |
|---|
| 93 | |
|---|
| 94 | !------------------------------------------------------------------------- |
|---|
| 95 | subroutine bilinear(f,f11,f21,f12,f22,x,x1,x2,y,y1,y2) |
|---|
| 96 | ! Used for interpolation of continuum data |
|---|
| 97 | |
|---|
| 98 | implicit none |
|---|
| 99 | |
|---|
| 100 | real*8 x,y,x1,x2,y1,y2 |
|---|
| 101 | real*8 f,f11,f12,f21,f22,fA,fB |
|---|
| 102 | |
|---|
| 103 | ! 1st in x-direction |
|---|
| 104 | fA=f11*(x2-x)/(x2-x1)+f21*(x-x1)/(x2-x1) |
|---|
| 105 | fB=f12*(x2-x)/(x2-x1)+f22*(x-x1)/(x2-x1) |
|---|
| 106 | |
|---|
| 107 | ! then in y-direction |
|---|
| 108 | f=fA*(y2-y)/(y2-y1)+fB*(y-y1)/(y2-y1) |
|---|
| 109 | |
|---|
| 110 | return |
|---|
| 111 | end subroutine bilinear |
|---|
| 112 | |
|---|
