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

Last change on this file since 5159 was 5119, checked in by abarral, 11 months ago

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