| [878] | 1 | subroutine bilinearbig(nX,nY,x_arr,y_arr,f2d_arr,x_in,y_in,f,ind) |
|---|
| [873] | 2 | |
|---|
| 3 | ! Necessary for interpolation of continuum data |
|---|
| 4 | ! optimized by A. Spiga 01/2013 |
|---|
| 5 | |
|---|
| 6 | implicit none |
|---|
| 7 | |
|---|
| [878] | 8 | integer nX,nY,i,j,ind,b |
|---|
| [873] | 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 | |
|---|
| 17 | integer strlen |
|---|
| 18 | character*100 label |
|---|
| 19 | label='subroutine bilinear' |
|---|
| 20 | |
|---|
| [918] | 21 | |
|---|
| [873] | 22 | x=x_in |
|---|
| 23 | y=y_in |
|---|
| 24 | |
|---|
| 25 | !! AS: important to optimize here because the array is quite large |
|---|
| 26 | !! ... and actually calculations only need to be done once |
|---|
| [903] | 27 | !! IF ind=-9999 we have not calculated yet |
|---|
| [878] | 28 | if ( ind == -9999) then |
|---|
| [873] | 29 | !1st check we're within the wavenumber range |
|---|
| 30 | if ((x.lt.x_arr(2)).or.(x.gt.x_arr(nX-2))) then |
|---|
| [878] | 31 | ind=-1 |
|---|
| [918] | 32 | return |
|---|
| [873] | 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) |
|---|
| [878] | 40 | ind=i-1 |
|---|
| [873] | 41 | end do |
|---|
| 42 | endif |
|---|
| [903] | 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 |
|---|
| [873] | 48 | f=0.0D+0 |
|---|
| 49 | return |
|---|
| [903] | 50 | !! Or we already determined ind -- so we just proceed |
|---|
| [873] | 51 | else |
|---|
| [878] | 52 | x1=x_arr(ind) |
|---|
| 53 | x2=x_arr(ind+1) |
|---|
| [873] | 54 | endif |
|---|
| 55 | |
|---|
| 56 | ! ... and for y within the temperature range |
|---|
| 57 | if ((y.lt.y_arr(1)).or.(y.gt.y_arr(nY))) then |
|---|
| 58 | write(*,*) 'Warning from bilinearH2H2:' |
|---|
| 59 | write(*,*) 'Outside continuum temperature range!' |
|---|
| 60 | if(y.lt.y_arr(1))then |
|---|
| 61 | y=y_arr(1)+0.01 |
|---|
| 62 | endif |
|---|
| 63 | if(y.gt.y_arr(nY))then |
|---|
| 64 | y=y_arr(nY)-0.01 |
|---|
| 65 | endif |
|---|
| 66 | else |
|---|
| 67 | j=1 |
|---|
| 68 | y2=y_arr(j) |
|---|
| 69 | do while ( y2 .le. y ) |
|---|
| 70 | y1=y2 |
|---|
| 71 | j=j+1 |
|---|
| 72 | y2=y_arr(j) |
|---|
| 73 | b=j-1 |
|---|
| 74 | end do |
|---|
| 75 | endif |
|---|
| [918] | 76 | |
|---|
| [878] | 77 | f11=f2d_arr(ind,b) |
|---|
| 78 | f21=f2d_arr(ind+1,b) |
|---|
| 79 | f12=f2d_arr(ind,b+1) |
|---|
| 80 | f22=f2d_arr(ind+1,b+1) |
|---|
| [873] | 81 | |
|---|
| 82 | call bilinear(f,f11,f21,f12,f22,x,x1,x2,y,y1,y2) |
|---|
| 83 | |
|---|
| 84 | return |
|---|
| 85 | end subroutine bilinearbig |
|---|
