source: LMDZ6/branches/Amaury_dev/libf/misc/lmdz_regr_lint.f90 @ 5117

Last change on this file since 5117 was 5117, checked in by abarral, 2 months ago

rename modules properly lmdz_*
move some unused files to obsolete/
(lint) uppercase fortran keywords
File size: 8.9 KB
Line 
1MODULE lmdz_regr_lint
2
3  USE lmdz_assert_eq, ONLY: assert_eq
4  USE lmdz_assert, ONLY: assert
5  USE lmdz_interpolation, ONLY: hunt
6
7  IMPLICIT NONE
8
9  ! Purpose: Each procedure regrids by linear interpolation along dimension "ix"
10  !          the input field "vs" to the output field "vt".
11  ! Remark:
12  !   * "vs" and "vt" have the same dimensions except Nr. ix (ns for vs, nt for vt)
13
14  ! Argument                         Type         Description
15  !-------------------------------------------------------------------------------
16  !  INTEGER, INTENT(IN)  :: ix      Scalar       dimension regridded <=rank(vs)
17  !  REAL,    INTENT(IN)  :: vs(*)   Rank>=1      source grid field values
18  !  REAL,    INTENT(IN)  :: xs(:)   Vector(ns)   centers of source grid, asc. order
19  !  REAL,    INTENT(IN)  :: xt(:)   Vector(nt)   centers of target grid, asc. order
20  !  REAL,    INTENT(OUT) :: vt(*)   Rank>=1      regridded field
21
22  INTERFACE regr_lint
23    ! The procedures differ only from the rank of the input/output fields.
24    MODULE PROCEDURE regr1_lint, regr2_lint, regr3_lint, &
25            regr4_lint, regr5_lint
26  END INTERFACE
27
28  PRIVATE
29  PUBLIC :: regr_lint
30
31CONTAINS
32
33
34  !-------------------------------------------------------------------------------
35
36  SUBROUTINE regr1_lint(ix, vs, xs, xt, vt)
37
38    !-------------------------------------------------------------------------------
39    ! Arguments:
40    INTEGER, INTENT(IN) :: ix
41    REAL, INTENT(IN) :: vs(:)
42    REAL, INTENT(IN) :: xs(:)
43    REAL, INTENT(IN) :: xt(:)
44    REAL, INTENT(OUT) :: vt(:)
45    !-------------------------------------------------------------------------------
46    ! Local variables:
47    INTEGER :: is, it, ns, nt, isb ! "is" bound between 1 and "ns - 1"
48    REAL :: r
49    !-------------------------------------------------------------------------------
50    CALL check_size(ix, SHAPE(vs), SHAPE(vt), SIZE(xs), SIZE(xt), ns, nt)
51    is = -1 ! go immediately to bisection on first CALL to "hunt"
52    DO it = 1, SIZE(xt)
53      CALL hunt(xs, xt(it), is); isb = MIN(MAX(is, 1), ns - 1)
54      r = (xs(isb + 1) - xt(it)) / (xs(isb + 1) - xs(isb))
55      vt(it) = r * vs(isb) + (1. - r) * vs(isb + 1)
56    END DO
57
58  END SUBROUTINE regr1_lint
59
60  !-------------------------------------------------------------------------------
61
62
63  !-------------------------------------------------------------------------------
64
65  SUBROUTINE regr2_lint(ix, vs, xs, xt, vt)
66
67    !-------------------------------------------------------------------------------
68    ! Arguments:
69    INTEGER, INTENT(IN) :: ix
70    REAL, INTENT(IN) :: vs(:, :)
71    REAL, INTENT(IN) :: xs(:)
72    REAL, INTENT(IN) :: xt(:)
73    REAL, INTENT(OUT) :: vt(:, :)
74    !-------------------------------------------------------------------------------
75    ! Local variables:
76    INTEGER :: is, it, ns, nt, isb ! "is" bound between 1 and "ns - 1"
77    REAL :: r
78    !-------------------------------------------------------------------------------
79    CALL check_size(ix, SHAPE(vs), SHAPE(vt), SIZE(xs), SIZE(xt), ns, nt)
80    is = -1 ! go immediately to bisection on first CALL to "hunt"
81    DO it = 1, SIZE(xt)
82      CALL hunt(xs, xt(it), is); isb = MIN(MAX(is, 1), ns - 1)
83      r = (xs(isb + 1) - xt(it)) / (xs(isb + 1) - xs(isb))
84      IF(ix==1) vt(it, :) = r * vs(isb, :) + (1. - r) * vs(isb + 1, :)
85      IF(ix==2) vt(:, it) = r * vs(:, isb) + (1. - r) * vs(:, isb + 1)
86    END DO
87
88  END SUBROUTINE regr2_lint
89
90  !-------------------------------------------------------------------------------
91
92
93  !-------------------------------------------------------------------------------
94
95  SUBROUTINE regr3_lint(ix, vs, xs, xt, vt)
96
97    !-------------------------------------------------------------------------------
98    ! Arguments:
99    INTEGER, INTENT(IN) :: ix
100    REAL, INTENT(IN) :: vs(:, :, :)
101    REAL, INTENT(IN) :: xs(:)
102    REAL, INTENT(IN) :: xt(:)
103    REAL, INTENT(OUT) :: vt(:, :, :)
104    !-------------------------------------------------------------------------------
105    ! Local variables:
106    INTEGER :: is, it, ns, nt, isb ! "is" bound between 1 and "ns - 1"
107    REAL :: r
108    !-------------------------------------------------------------------------------
109    CALL check_size(ix, SHAPE(vs), SHAPE(vt), SIZE(xs), SIZE(xt), ns, nt)
110    is = -1 ! go immediately to bisection on first CALL to "hunt"
111    DO it = 1, SIZE(xt)
112      CALL hunt(xs, xt(it), is); isb = MIN(MAX(is, 1), ns - 1)
113      r = (xs(isb + 1) - xt(it)) / (xs(isb + 1) - xs(isb))
114      IF(ix==1) vt(it, :, :) = r * vs(isb, :, :) + (1. - r) * vs(isb + 1, :, :)
115      IF(ix==2) vt(:, it, :) = r * vs(:, isb, :) + (1. - r) * vs(:, isb + 1, :)
116      IF(ix==3) vt(:, :, it) = r * vs(:, :, isb) + (1. - r) * vs(:, :, isb + 1)
117    END DO
118
119  END SUBROUTINE regr3_lint
120
121  !-------------------------------------------------------------------------------
122
123
124  !-------------------------------------------------------------------------------
125
126  SUBROUTINE regr4_lint(ix, vs, xs, xt, vt)
127
128    !-------------------------------------------------------------------------------
129    ! Arguments:
130    INTEGER, INTENT(IN) :: ix
131    REAL, INTENT(IN) :: vs(:, :, :, :)
132    REAL, INTENT(IN) :: xs(:)
133    REAL, INTENT(IN) :: xt(:)
134    REAL, INTENT(OUT) :: vt(:, :, :, :)
135    !-------------------------------------------------------------------------------
136    ! Local variables:
137    INTEGER :: is, it, ns, nt, isb ! "is" bound between 1 and "ns - 1"
138    REAL :: r
139    !-------------------------------------------------------------------------------
140    CALL check_size(ix, SHAPE(vs), SHAPE(vt), SIZE(xs), SIZE(xt), ns, nt)
141    is = -1 ! go immediately to bisection on first CALL to "hunt"
142    DO it = 1, SIZE(xt)
143      CALL hunt(xs, xt(it), is); isb = MIN(MAX(is, 1), ns - 1)
144      r = (xs(isb + 1) - xt(it)) / (xs(isb + 1) - xs(isb))
145      IF(ix==1) vt(it, :, :, :) = r * vs(isb, :, :, :) + (1. - r) * vs(isb + 1, :, :, :)
146      IF(ix==2) vt(:, it, :, :) = r * vs(:, isb, :, :) + (1. - r) * vs(:, isb + 1, :, :)
147      IF(ix==3) vt(:, :, it, :) = r * vs(:, :, isb, :) + (1. - r) * vs(:, :, isb + 1, :)
148      IF(ix==4) vt(:, :, :, it) = r * vs(:, :, :, isb) + (1. - r) * vs(:, :, :, isb + 1)
149    END DO
150
151  END SUBROUTINE regr4_lint
152
153  !-------------------------------------------------------------------------------
154
155
156  !-------------------------------------------------------------------------------
157
158  SUBROUTINE regr5_lint(ix, vs, xs, xt, vt)
159
160    !-------------------------------------------------------------------------------
161    ! Arguments:
162    INTEGER, INTENT(IN) :: ix
163    REAL, INTENT(IN) :: vs(:, :, :, :, :)
164    REAL, INTENT(IN) :: xs(:)
165    REAL, INTENT(IN) :: xt(:)
166    REAL, INTENT(OUT) :: vt(:, :, :, :, :)
167    !-------------------------------------------------------------------------------
168    ! Local variables:
169    INTEGER :: is, it, ns, nt, isb ! "is" bound between 1 and "ns - 1"
170    REAL :: r
171    !-------------------------------------------------------------------------------
172    CALL check_size(ix, SHAPE(vs), SHAPE(vt), SIZE(xs), SIZE(xt), ns, nt)
173    is = -1 ! go immediately to bisection on first CALL to "hunt"
174    DO it = 1, SIZE(xt)
175      CALL hunt(xs, xt(it), is); isb = MIN(MAX(is, 1), ns - 1)
176      r = (xs(isb + 1) - xt(it)) / (xs(isb + 1) - xs(isb))
177      IF(ix==1) vt(it, :, :, :, :) = r * vs(isb, :, :, :, :) + (1. - r) * vs(isb + 1, :, :, :, :)
178      IF(ix==2) vt(:, it, :, :, :) = r * vs(:, isb, :, :, :) + (1. - r) * vs(:, isb + 1, :, :, :)
179      IF(ix==3) vt(:, :, it, :, :) = r * vs(:, :, isb, :, :) + (1. - r) * vs(:, :, isb + 1, :, :)
180      IF(ix==4) vt(:, :, :, it, :) = r * vs(:, :, :, isb, :) + (1. - r) * vs(:, :, :, isb + 1, :)
181      IF(ix==5) vt(:, :, :, :, it) = r * vs(:, :, :, :, isb) + (1. - r) * vs(:, :, :, :, isb + 1)
182    END DO
183
184  END SUBROUTINE regr5_lint
185
186  !-------------------------------------------------------------------------------
187
188
189  !-------------------------------------------------------------------------------
190
191  SUBROUTINE check_size(ix, svs, svt, nxs, nxt, ns, nt)
192
193    !-------------------------------------------------------------------------------
194    ! Arguments:
195    INTEGER, INTENT(IN) :: ix, svs(:), svt(:), nxs, nxt
196    INTEGER, INTENT(OUT) :: ns, nt
197    !-------------------------------------------------------------------------------
198    ! Local variables:
199    INTEGER :: rk, is
200    CHARACTER(LEN = 80) :: sub, msg
201    !-------------------------------------------------------------------------------
202    rk = SIZE(svs)
203    WRITE(sub, '(a,2i0,a)')"regr", rk, ix, "_lint"
204    CALL assert(ix>=1.AND.ix<=rk, TRIM(sub) // ": ix exceeds fields rank")
205    DO is = 1, rk; IF(is==ix) CYCLE
206    WRITE(msg, '(a,i1)')TRIM(sub) // " n", is
207    CALL assert(svs(is)==svt(is), msg)
208    END DO
209    ns = assert_eq(svs(ix), nxs, TRIM(sub) // " ns")
210    nt = assert_eq(svt(ix), nxt, TRIM(sub) // " nt")
211
212  END SUBROUTINE check_size
213
214  !-------------------------------------------------------------------------------
215
216END MODULE lmdz_regr_lint
217
218!-------------------------------------------------------------------------------
219
Note: See TracBrowser for help on using the repository browser.