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

Last change on this file since 5139 was 5119, checked in by abarral, 4 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.