| 1 | MODULE mod_bicub |
|---|
| 2 | CONTAINS |
|---|
| 3 | SUBROUTINE bicub ( pout, px, py, kndx, kndy, pax, pay, kpts |
|---|
| 4 | $ , pin, ki, kj) |
|---|
| 5 | C**** |
|---|
| 6 | C ***************************** |
|---|
| 7 | C * OASIS ROUTINE - LEVEL 3 * |
|---|
| 8 | C * ------------- ------- * |
|---|
| 9 | C ***************************** |
|---|
| 10 | C |
|---|
| 11 | C**** *bicub* - Bicubic interpolation |
|---|
| 12 | C |
|---|
| 13 | C Purpose: Proceed to a bicubic interpolation |
|---|
| 14 | C ------- |
|---|
| 15 | C Warning : this interpolation is correct ONLY for a regular source grid |
|---|
| 16 | C |
|---|
| 17 | C |
|---|
| 18 | C** Interface: |
|---|
| 19 | C --------- |
|---|
| 20 | C *CALL* *bicub* ( zo, px, py, kndx, kndy, pax, pay, kpts |
|---|
| 21 | C $ , z, ki, kj) |
|---|
| 22 | C |
|---|
| 23 | C** Method |
|---|
| 24 | C ------ |
|---|
| 25 | C |
|---|
| 26 | C * * * * |
|---|
| 27 | C |
|---|
| 28 | C * * * * |
|---|
| 29 | C # ==> pt (x,y) |
|---|
| 30 | C * (=) * * ==> = pt (kndx, kndy) |
|---|
| 31 | C |
|---|
| 32 | C * * * * |
|---|
| 33 | C |
|---|
| 34 | C Input: |
|---|
| 35 | C ----- |
|---|
| 36 | C px : longitudes of target grid |
|---|
| 37 | C py : latitudes of target grid |
|---|
| 38 | C kndx : index of source point in source longitude |
|---|
| 39 | C kndy : index of source point in source latitude |
|---|
| 40 | C pax : longitudes of source grid |
|---|
| 41 | C pay : latitudes of source grid |
|---|
| 42 | C kpts : dimension of target grid |
|---|
| 43 | C pin : input field on source grid |
|---|
| 44 | C ki, kj : dimension of source grid |
|---|
| 45 | C |
|---|
| 46 | C Output: |
|---|
| 47 | C ------ |
|---|
| 48 | C pout : interpolated field on target grid |
|---|
| 49 | C |
|---|
| 50 | C Workspace: |
|---|
| 51 | C --------- |
|---|
| 52 | C Local variables |
|---|
| 53 | C zy1, zy2, zy3, zy4, i, j, zdx, zdy ,z1, z2, z3, z4 |
|---|
| 54 | C Statement function |
|---|
| 55 | C cubic |
|---|
| 56 | C |
|---|
| 57 | C Externals: |
|---|
| 58 | C --------- |
|---|
| 59 | C None |
|---|
| 60 | C |
|---|
| 61 | C Reference: |
|---|
| 62 | C --------- |
|---|
| 63 | C See OASIS manual (1995) |
|---|
| 64 | C |
|---|
| 65 | C History: |
|---|
| 66 | C ------- |
|---|
| 67 | C Version Programmer Date Description |
|---|
| 68 | C ------- ---------- ---- ----------- |
|---|
| 69 | C 2.0 O. Marti 96/07/15 Created |
|---|
| 70 | C |
|---|
| 71 | C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
|---|
| 72 | C |
|---|
| 73 | C * |
|---|
| 74 | IMPLICIT NONE |
|---|
| 75 | C |
|---|
| 76 | C* ---------------------------- Argument declarations ------------------- |
|---|
| 77 | C |
|---|
| 78 | INTEGER, INTENT ( in) :: kpts, ki, kj |
|---|
| 79 | REAL, DIMENSION ( kpts), INTENT ( in) :: px, py |
|---|
| 80 | REAL, DIMENSION ( -1: ki + 2), INTENT ( in) :: pax |
|---|
| 81 | REAL, DIMENSION ( -1: kj + 2), INTENT ( in) :: pay |
|---|
| 82 | REAL, DIMENSION ( -1: ki + 2, -1: kj + 2), INTENT ( in) :: pin |
|---|
| 83 | INTEGER, DIMENSION ( kpts), INTENT ( in) :: kndx, kndy |
|---|
| 84 | REAL, DIMENSION ( kpts), INTENT ( out) :: pout |
|---|
| 85 | C |
|---|
| 86 | C* ---------------------------- Local declarations ---------------------- |
|---|
| 87 | C |
|---|
| 88 | REAL :: zy1, zy2, zy3, zy4 |
|---|
| 89 | INTEGER :: jn, ji, jj |
|---|
| 90 | C |
|---|
| 91 | C* ---------------------------- Statement fucntions---------------------- |
|---|
| 92 | C |
|---|
| 93 | REAL :: cubic, zdx, zdy , z1, z2, z3, z4 |
|---|
| 94 | C |
|---|
| 95 | cubic ( z1, z2, z3, z4, zdx) = (((( z4 - z1)*0.1666666666666 |
|---|
| 96 | $ + 0.5 * ( z2 - z3) |
|---|
| 97 | $ ) * zdx + 0.5 * ( z1 + z3) - z2) * zdx + z3 |
|---|
| 98 | $ - 0.1666666666666 * z4 - 0.5 * z2 - 0. |
|---|
| 99 | $ 3333333333333 * z1) * zdx + z2 |
|---|
| 100 | C |
|---|
| 101 | C* ---------------------------- Poema verses ---------------------------- |
|---|
| 102 | C |
|---|
| 103 | C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
|---|
| 104 | C |
|---|
| 105 | C* 1. Interpolation |
|---|
| 106 | C ------------- |
|---|
| 107 | C |
|---|
| 108 | DO jn = 1, kpts |
|---|
| 109 | ji = kndx ( jn) |
|---|
| 110 | jj = kndy ( jn) |
|---|
| 111 | zdx = ( px ( jn) - pax ( ji)) / ( pax ( ji + 1) - pax ( ji)) |
|---|
| 112 | zdy = ( py ( jn) - pay ( jj)) / ( pay ( jj + 1) - pay ( jj)) |
|---|
| 113 | zy1 = cubic ( pin ( ji-1, jj-1), pin ( ji , jj-1) |
|---|
| 114 | $ , pin ( ji+1, jj-1), pin ( ji+2, jj-1), zdx) |
|---|
| 115 | zy2 = cubic ( pin ( ji-1, jj ), pin ( ji , jj ) |
|---|
| 116 | $ , pin ( ji+1, jj ), pin ( ji+2, jj ), zdx) |
|---|
| 117 | zy3 = cubic ( pin ( ji-1, jj+1), pin ( ji , jj+1) |
|---|
| 118 | $ , pin ( ji+1, jj+1), pin ( ji+2, jj+1), zdx) |
|---|
| 119 | zy4 = cubic ( pin ( ji-1, jj+2), pin ( ji , jj+2) |
|---|
| 120 | $ , pin ( ji+1, jj+2), pin ( ji+2, jj+2), zdx) |
|---|
| 121 | pout ( jn) = cubic ( zy1, zy2, zy3, zy4, zdy) |
|---|
| 122 | END DO |
|---|
| 123 | C |
|---|
| 124 | C |
|---|
| 125 | C* 3. End of routine |
|---|
| 126 | C -------------- |
|---|
| 127 | C |
|---|
| 128 | RETURN |
|---|
| 129 | END SUBROUTINE bicub |
|---|
| 130 | END MODULE mod_bicub |
|---|
