source: LMDZ6/branches/contrails/libf/misc/regr_lint_m.f90 @ 5426

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

.f90 <-> .F90 depending on cpp key use
File size: 8.1 KB
Line 
1MODULE regr_lint_m
2
3  USE assert_eq_m,   ONLY: assert_eq
4  USE assert_m,      ONLY: assert
5  USE 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!
36SUBROUTINE 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
58END SUBROUTINE regr1_lint
59!
60!-------------------------------------------------------------------------------
61
62
63!-------------------------------------------------------------------------------
64!
65SUBROUTINE 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
88END SUBROUTINE regr2_lint
89!
90!-------------------------------------------------------------------------------
91
92
93!-------------------------------------------------------------------------------
94!
95SUBROUTINE 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
119END SUBROUTINE regr3_lint
120!
121!-------------------------------------------------------------------------------
122
123
124!-------------------------------------------------------------------------------
125!
126SUBROUTINE 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
151END SUBROUTINE regr4_lint
152!
153!-------------------------------------------------------------------------------
154
155
156!-------------------------------------------------------------------------------
157!
158SUBROUTINE 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
184END SUBROUTINE regr5_lint
185!
186!-------------------------------------------------------------------------------
187
188
189!-------------------------------------------------------------------------------
190!
191SUBROUTINE 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
212END SUBROUTINE check_size
213!
214!-------------------------------------------------------------------------------
215
216END MODULE regr_lint_m
217!
218!-------------------------------------------------------------------------------
219
Note: See TracBrowser for help on using the repository browser.