source: LMDZ5/branches/testing/libf/misc/regr_conserv_m.F90 @ 5326

Last change on this file since 5326 was 2839, checked in by Laurent Fairhead, 8 years ago

Merged trunk changes r2785:2838 into testing branch
File size: 15.4 KB
Line 
1MODULE regr_conserv_m
2
3  USE assert_eq_m,   ONLY: assert_eq
4  USE assert_m,      ONLY: assert
5  USE interpolation, ONLY: locate
6
7  IMPLICIT NONE
8
9! Purpose: Each procedure regrids a piecewise linear function (not necessarily
10!          continuous) by averaging it. This is a conservative regridding.
11!          The regridding operation is done along dimension "ix" of the input
12!          field "vs". Input are positions of cell edges.
13! Remarks:
14!   * The target grid should be included in the source grid:
15!     no extrapolation is allowed.
16!   * Field on target grid "vt" has same rank, slopes and averages as "vs".
17!   * If optional argument "slope" is not given, 0 is assumed for slopes.
18!     Then the regridding is first order instead of second.
19!   * "vs" and "vt" have the same dimensions except Nr. ix (ns for vs, nt for vt)
20
21! Argument                         Type         Description
22!-------------------------------------------------------------------------------
23!  INTEGER, INTENT(IN)  :: ix      Scalar       dimension regridded <=rank(vs)
24!  REAL,    INTENT(IN)  :: vs(*)   Rank>=1      averages in source grid cells
25!  REAL,    INTENT(IN)  :: xs(:)   Vector(ns+1) edges of source grid, asc. order
26!  REAL,    INTENT(IN)  :: xt(:)   Vector(nt+1) edges of target grid, asc. order
27!  REAL,    INTENT(OUT) :: vt(*)   Rank>=1      regridded field
28! [REAL,    INTENT(IN)  :: slope]  Rank>=1      slopes in source grid cells
29
30  INTERFACE regr_conserv
31    ! The procedures differ only from the rank of the input/output fields.
32    MODULE PROCEDURE regr1_conserv, regr2_conserv, regr3_conserv, &
33                     regr4_conserv, regr5_conserv
34  END INTERFACE
35
36  PRIVATE
37  PUBLIC :: regr_conserv
38
39CONTAINS
40
41!-------------------------------------------------------------------------------
42!
43SUBROUTINE regr1_conserv(ix, vs, xs, xt, vt, slope)
44!
45!-------------------------------------------------------------------------------
46! Arguments:
47  INTEGER,        INTENT(IN)  :: ix
48  REAL,           INTENT(IN)  :: vs(:)
49  REAL,           INTENT(IN)  :: xs(:), xt(:)
50  REAL,           INTENT(OUT) :: vt(:)
51  REAL, OPTIONAL, INTENT(IN)  :: slope(:)
52!-------------------------------------------------------------------------------
53! Local variables:
54  INTEGER :: is, it, ns, nt
55  REAL    :: co, idt
56  LOGICAL :: lslope
57!-------------------------------------------------------------------------------
58  lslope=PRESENT(slope)
59  CALL check_size(ix,SHAPE(vs),SHAPE(vt),xs,xt,ns,nt)
60  is=locate(xs,xt(1))                    !--- 1<= is <= ns (no extrapolation)
61  vt(:)=0.
62  DO it=1, nt; idt=1./(xt(it+1)-xt(it))
63    IF(xt(it+1)<=xs(is+1)) THEN          !--- xs(is)<=xt(it)<=xt(it+1)<=xs(is+1)
64      CALL mean_lin(xt(it),xt(it+1))
65    ELSE
66      CALL mean_lin(xt(it),xs(is+1)); is=is+1
67      DO WHILE(xs(is+1)<xt(it+1))        !--- 1<=is<=ns-1
68      CALL mean_lin(xs(is),xs(is+1)); is=is+1
69      END DO
70      CALL mean_lin(xs(is),xt(it+1))     !--- 1<=is<=ns
71    END IF
72    IF(xs(is+1)==xt(it+1)) is=is+1       !--- 1<=is<=ns .OR. it==nt
73  END DO
74
75CONTAINS
76
77!-------------------------------------------------------------------------------
78SUBROUTINE mean_lin(a,b)  ! mean [a,b] of the linear function in [xs(is),xs(is+1)]
79!-------------------------------------------------------------------------------
80! Arguments:
81  REAL, INTENT(IN)    :: a, b
82!-------------------------------------------------------------------------------
83  IF(lslope.AND.(a/=xs(is).OR.b/=xs(is+1))) THEN; co=0.
84    IF(a==xt(it  )) co=co+xt(it  )-xs(is  )
85    IF(b==xt(it+1)) co=co+xt(it+1)-xs(is+1)
86    vt(it) = vt(it)+idt*(b-a)*(vs(is)+co*slope(is))
87  ELSE
88    vt(it) = vt(it)+idt*(b-a)* vs(is)
89  END IF
90
91END SUBROUTINE mean_lin
92!-------------------------------------------------------------------------------
93
94END SUBROUTINE regr1_conserv
95!
96!-------------------------------------------------------------------------------
97
98
99!-------------------------------------------------------------------------------
100!
101SUBROUTINE regr2_conserv(ix, vs, xs, xt, vt, slope)
102!
103!-------------------------------------------------------------------------------
104! Arguments:
105  INTEGER,        INTENT(IN)  :: ix
106  REAL,           INTENT(IN)  :: vs(:,:)
107  REAL,           INTENT(IN)  :: xs(:), xt(:)
108  REAL,           INTENT(OUT) :: vt(:,:)
109  REAL, OPTIONAL, INTENT(IN)  :: slope(:,:)
110!-------------------------------------------------------------------------------
111! Local variables:
112  INTEGER :: is, it, ns, nt
113  REAL    :: co, idt
114  LOGICAL :: lslope
115!-------------------------------------------------------------------------------
116  lslope=PRESENT(slope)
117  CALL check_size(ix,SHAPE(vs),SHAPE(vt),xs,xt,ns,nt)
118  is=locate(xs,xt(1)) ! 1 <= is <= ns, because we forbid extrapolation
119  vt(:,:)=0.
120  DO it=1, nt; idt=1./(xt(it+1)-xt(it))
121    IF(xt(it+1)<=xs(is+1)) THEN          !--- xs(is)<=xt(it)<=xt(it+1)<=xs(is+1)
122      CALL mean_lin(xt(it),xt(it+1))
123    ELSE
124      CALL mean_lin(xt(it),xs(is+1)); is=is+1
125      DO WHILE(xs(is+1)<xt(it+1))        !--- 1<=is<=ns-1
126      CALL mean_lin(xs(is),xs(is+1)); is=is+1
127      END DO
128      CALL mean_lin(xs(is),xt(it+1))     !--- 1<=is<=ns
129    END IF
130    IF(xs(is+1)==xt(it+1)) is=is+1       !--- 1<=is<=ns .OR. it==nt
131  END DO
132
133CONTAINS
134
135!-------------------------------------------------------------------------------
136SUBROUTINE mean_lin(a,b)  ! mean [a,b] of the linear function in [xs(is),xs(is+1)]
137!-------------------------------------------------------------------------------
138! Arguments:
139  REAL, INTENT(IN)    :: a, b
140!-------------------------------------------------------------------------------
141  IF(lslope.AND.(a/=xs(is).OR.b/=xs(is+1))) THEN; co=0.
142    IF(a==xt(it  )) co=co+xt(it  )-xs(is  )
143    IF(b==xt(it+1)) co=co+xt(it+1)-xs(is+1)
144    IF(ix==1) vt(it,:) = vt(it,:)+idt*(b-a)*(vs(is,:)+co*slope(is,:))
145    IF(ix==2) vt(:,it) = vt(:,it)+idt*(b-a)*(vs(:,is)+co*slope(:,is))
146  ELSE
147    IF(ix==1) vt(it,:) = vt(it,:)+idt*(b-a)* vs(is,:)
148    IF(ix==2) vt(:,it) = vt(:,it)+idt*(b-a)* vs(:,is)
149  END IF
150
151END SUBROUTINE mean_lin
152!-------------------------------------------------------------------------------
153
154END SUBROUTINE regr2_conserv
155!
156!-------------------------------------------------------------------------------
157
158
159!-------------------------------------------------------------------------------
160!
161SUBROUTINE regr3_conserv(ix, vs, xs, xt, vt, slope)
162!
163!-------------------------------------------------------------------------------
164! Arguments:
165  INTEGER,        INTENT(IN)  :: ix
166  REAL,           INTENT(IN)  :: vs(:,:,:)
167  REAL,           INTENT(IN)  :: xs(:), xt(:)
168  REAL,           INTENT(OUT) :: vt(:,:,:)
169  REAL, OPTIONAL, INTENT(IN)  :: slope(:,:,:)
170!-------------------------------------------------------------------------------
171! Local variables:
172  INTEGER :: is, it, ns, nt
173  REAL    :: co, idt
174  LOGICAL :: lslope
175!-------------------------------------------------------------------------------
176  lslope=PRESENT(slope)
177  CALL check_size(ix,SHAPE(vs),SHAPE(vt),xs,xt,ns,nt)
178  is=locate(xs,xt(1)) ! 1 <= is <= ns, because we forbid extrapolation
179  vt(:,:,:)=0.
180  DO it=1, nt; idt=1./(xt(it+1)-xt(it))
181    IF(xt(it+1)<=xs(is+1)) THEN          !--- xs(is)<=xt(it)<=xt(it+1)<=xs(is+1)
182      CALL mean_lin(xt(it),xt(it+1))
183    ELSE
184      CALL mean_lin(xt(it),xs(is+1)); is=is+1
185      DO WHILE(xs(is+1)<xt(it+1))        !--- 1<=is<=ns-1
186      CALL mean_lin(xs(is),xs(is+1)); is=is+1
187      END DO
188      CALL mean_lin(xs(is),xt(it+1))     !--- 1<=is<=ns
189    END IF
190    IF(xs(is+1)==xt(it+1)) is=is+1       !--- 1<=is<=ns .OR. it==nt
191  END DO
192
193CONTAINS
194
195!-------------------------------------------------------------------------------
196SUBROUTINE mean_lin(a,b)  ! mean [a,b] of the linear function in [xs(is),xs(is+1)]
197!-------------------------------------------------------------------------------
198! Arguments:
199  REAL, INTENT(IN)    :: a, b
200!-------------------------------------------------------------------------------
201  IF(lslope.AND.(a/=xs(is).OR.b/=xs(is+1))) THEN; co=0.
202    IF(a==xt(it  )) co=co+xt(it  )-xs(is  )
203    IF(b==xt(it+1)) co=co+xt(it+1)-xs(is+1)
204    IF(ix==1) vt(it,:,:) = vt(it,:,:)+idt*(b-a)*(vs(is,:,:)+co*slope(is,:,:))
205    IF(ix==2) vt(:,it,:) = vt(:,it,:)+idt*(b-a)*(vs(:,is,:)+co*slope(:,is,:))
206    IF(ix==3) vt(:,:,it) = vt(:,:,it)+idt*(b-a)*(vs(:,:,is)+co*slope(:,:,is))
207  ELSE
208    IF(ix==1) vt(it,:,:) = vt(it,:,:)+idt*(b-a)* vs(is,:,:)
209    IF(ix==2) vt(:,it,:) = vt(:,it,:)+idt*(b-a)* vs(:,is,:)
210    IF(ix==3) vt(:,:,it) = vt(:,:,it)+idt*(b-a)* vs(:,:,is)
211  END IF
212
213END SUBROUTINE mean_lin
214!-------------------------------------------------------------------------------
215
216END SUBROUTINE regr3_conserv
217!
218!-------------------------------------------------------------------------------
219
220
221!-------------------------------------------------------------------------------
222!
223SUBROUTINE regr4_conserv(ix, vs, xs, xt, vt, slope)
224!
225!-------------------------------------------------------------------------------
226! Arguments:
227  INTEGER,        INTENT(IN)  :: ix
228  REAL,           INTENT(IN)  :: vs(:,:,:,:)
229  REAL,           INTENT(IN)  :: xs(:), xt(:)
230  REAL,           INTENT(OUT) :: vt(:,:,:,:)
231  REAL, OPTIONAL, INTENT(IN)  :: slope(:,:,:,:)
232!-------------------------------------------------------------------------------
233! Local variables:
234  INTEGER :: is, it, ns, nt
235  REAL    :: co, idt
236  LOGICAL :: lslope
237!-------------------------------------------------------------------------------
238  lslope=PRESENT(slope)
239  CALL check_size(ix,SHAPE(vs),SHAPE(vt),xs,xt,ns,nt)
240  is=locate(xs,xt(1)) ! 1 <= is <= ns, because we forbid extrapolation
241  vt(:,:,:,:)=0.
242  DO it=1, nt; idt=1./(xt(it+1)-xt(it))
243    IF(xt(it+1)<=xs(is+1)) THEN          !--- xs(is)<=xt(it)<=xt(it+1)<=xs(is+1)
244      CALL mean_lin(xt(it),xt(it+1))
245    ELSE
246      CALL mean_lin(xt(it),xs(is+1)); is=is+1
247      DO WHILE(xs(is+1)<xt(it+1))        !--- 1<=is<=ns-1
248      CALL mean_lin(xs(is),xs(is+1)); is=is+1
249      END DO
250      CALL mean_lin(xs(is),xt(it+1))     !--- 1<=is<=ns
251    END IF
252    IF(xs(is+1)==xt(it+1)) is=is+1       !--- 1<=is<=ns .OR. it==nt
253  END DO
254
255CONTAINS
256
257!-------------------------------------------------------------------------------
258SUBROUTINE mean_lin(a,b)  ! mean [a,b] of the linear function in [xs(is),xs(is+1)]
259!-------------------------------------------------------------------------------
260! Arguments:
261  REAL, INTENT(IN) :: a, b
262!-------------------------------------------------------------------------------
263  IF(lslope.AND.(a/=xs(is).OR.b/=xs(is+1))) THEN; co=0.
264    IF(a==xt(it  )) co=co+xt(it  )-xs(is  )
265    IF(b==xt(it+1)) co=co+xt(it+1)-xs(is+1)
266    IF(ix==1) vt(it,:,:,:) = vt(it,:,:,:)+idt*(b-a)*(vs(is,:,:,:)+co*slope(is,:,:,:))
267    IF(ix==2) vt(:,it,:,:) = vt(:,it,:,:)+idt*(b-a)*(vs(:,is,:,:)+co*slope(:,is,:,:))
268    IF(ix==3) vt(:,:,it,:) = vt(:,:,it,:)+idt*(b-a)*(vs(:,:,is,:)+co*slope(:,:,is,:))
269    IF(ix==4) vt(:,:,:,it) = vt(:,:,:,it)+idt*(b-a)*(vs(:,:,:,is)+co*slope(:,:,:,is))
270  ELSE
271    IF(ix==1) vt(it,:,:,:) = vt(it,:,:,:)+idt*(b-a)* vs(is,:,:,:)
272    IF(ix==2) vt(:,it,:,:) = vt(:,it,:,:)+idt*(b-a)* vs(:,is,:,:)
273    IF(ix==3) vt(:,:,it,:) = vt(:,:,it,:)+idt*(b-a)* vs(:,:,is,:)
274    IF(ix==4) vt(:,:,:,it) = vt(:,:,:,it)+idt*(b-a)* vs(:,:,:,is)
275  END IF
276
277END SUBROUTINE mean_lin
278!-------------------------------------------------------------------------------
279
280END SUBROUTINE regr4_conserv
281!
282!-------------------------------------------------------------------------------
283
284
285!-------------------------------------------------------------------------------
286!
287SUBROUTINE regr5_conserv(ix, vs, xs, xt, vt, slope)
288!
289!-------------------------------------------------------------------------------
290! Arguments:
291  INTEGER,        INTENT(IN)  :: ix
292  REAL,           INTENT(IN)  :: vs(:,:,:,:,:)
293  REAL,           INTENT(IN)  :: xs(:), xt(:)
294  REAL,           INTENT(OUT) :: vt(:,:,:,:,:)
295  REAL, OPTIONAL, INTENT(IN)  :: slope(:,:,:,:,:)
296!-------------------------------------------------------------------------------
297! Local variables:
298  INTEGER :: is, it, ns, nt
299  REAL    :: co, idt
300  LOGICAL :: lslope
301!-------------------------------------------------------------------------------
302  lslope=PRESENT(slope)
303  CALL check_size(ix,SHAPE(vs),SHAPE(vt),xs,xt,ns,nt)
304  is=locate(xs,xt(1)) ! 1 <= is <= ns, because we forbid extrapolation
305  vt(:,:,:,:,:)=0.
306  DO it=1, nt; idt=1./(xt(it+1)-xt(it))
307    IF(xt(it+1)<=xs(is+1)) THEN          !--- xs(is)<=xt(it)<=xt(it+1)<=xs(is+1)
308      CALL mean_lin(xt(it),xt(it+1))
309    ELSE
310      CALL mean_lin(xt(it),xs(is+1)); is=is+1
311      DO WHILE(xs(is+1)<xt(it+1))        !--- 1<=is<=ns-1
312      CALL mean_lin(xs(is),xs(is+1)); is=is+1
313      END DO
314      CALL mean_lin(xs(is),xt(it+1))     !--- 1<=is<=ns
315    END IF
316    IF(xs(is+1)==xt(it+1)) is=is+1     !--- 1<=is<=ns .OR. it==nt
317  END DO
318
319CONTAINS
320
321!-------------------------------------------------------------------------------
322SUBROUTINE mean_lin(a,b)  ! mean [a,b] of the linear function in [xs(is),xs(is+1)]
323!-------------------------------------------------------------------------------
324! Arguments:
325  REAL, INTENT(IN) :: a, b
326!-------------------------------------------------------------------------------
327  IF(lslope.AND.(a/=xs(is).OR.b/=xs(is+1))) THEN; co=0.
328    IF(a==xt(it  )) co=co+xt(it  )-xs(is  )
329    IF(b==xt(it+1)) co=co+xt(it+1)-xs(is+1)
330    IF(ix==1) vt(it,:,:,:,:) = vt(it,:,:,:,:)+idt*(b-a)*(vs(is,:,:,:,:)+co*slope(is,:,:,:,:))
331    IF(ix==2) vt(:,it,:,:,:) = vt(:,it,:,:,:)+idt*(b-a)*(vs(:,is,:,:,:)+co*slope(:,is,:,:,:))
332    IF(ix==3) vt(:,:,it,:,:) = vt(:,:,it,:,:)+idt*(b-a)*(vs(:,:,is,:,:)+co*slope(:,:,is,:,:))
333    IF(ix==4) vt(:,:,:,it,:) = vt(:,:,:,it,:)+idt*(b-a)*(vs(:,:,:,is,:)+co*slope(:,:,:,is,:))
334    IF(ix==5) vt(:,:,:,:,it) = vt(:,:,:,:,it)+idt*(b-a)*(vs(:,:,:,:,is)+co*slope(:,:,:,:,is))
335  ELSE
336    IF(ix==1) vt(it,:,:,:,:) = vt(it,:,:,:,:)+idt*(b-a)* vs(is,:,:,:,:)
337    IF(ix==2) vt(:,it,:,:,:) = vt(:,it,:,:,:)+idt*(b-a)* vs(:,is,:,:,:)
338    IF(ix==3) vt(:,:,it,:,:) = vt(:,:,it,:,:)+idt*(b-a)* vs(:,:,is,:,:)
339    IF(ix==4) vt(:,:,:,it,:) = vt(:,:,:,it,:)+idt*(b-a)* vs(:,:,:,is,:)
340    IF(ix==5) vt(:,:,:,:,it) = vt(:,:,:,:,it)+idt*(b-a)* vs(:,:,:,:,is)
341  END IF
342
343END SUBROUTINE mean_lin
344!-------------------------------------------------------------------------------
345
346END SUBROUTINE regr5_conserv
347!
348!-------------------------------------------------------------------------------
349
350
351!-------------------------------------------------------------------------------
352!
353SUBROUTINE check_size(ix,svs,svt,xs,xt,ns,nt)
354!
355!-------------------------------------------------------------------------------
356! Arguments:
357  INTEGER, INTENT(IN)  :: ix, svs(:), svt(:)
358  REAL,    INTENT(IN)  :: xs(:), xt(:)
359  INTEGER, INTENT(OUT) :: ns,    nt
360!-------------------------------------------------------------------------------
361! Local variables:
362  INTEGER :: rk, is, ii, n
363  CHARACTER(LEN=80) :: sub, msg
364!-------------------------------------------------------------------------------
365  WRITE(sub,'(a,2i0,a)')"regr",SIZE(svs),ix,"_conserv"
366  rk=assert_eq(SIZE(svs),SIZE(svt),TRIM(sub)//': inconsistent ranks')
367  WRITE(msg,'(a,2(i0,a))')TRIM(sub)//": ix (",ix,") exceeds field rank (",rk,")"
368  CALL assert(ix>=1.AND.ix<=rk,msg)
369  ii=0
370  DO is=1,rk; IF(is==ix) CYCLE
371    WRITE(msg,'(a,i0)')TRIM(sub)//" n",is
372    n=assert_eq(svs(is),svt(is),msg)
373    ii=ii+1
374  END DO
375  ns=assert_eq(svs(ix),SIZE(xs)-1,TRIM(sub)//" ns")
376  nt=assert_eq(svt(ix),SIZE(xt)-1,TRIM(sub)//" nt")
377
378  !--- Quick check on sort order:
379  CALL assert(xs(1)<xs(2), TRIM(sub)//": xs bad order")
380  CALL assert(xt(1)<xt(2), TRIM(sub)//": xt bad order")
381  CALL assert(xs(1)<=xt(1).AND.xt(nt+1)<=xs(ns+1), TRIM(sub)//": extrapolation")
382
383END SUBROUTINE check_size
384!
385!-------------------------------------------------------------------------------
386
387
388END MODULE regr_conserv_m
Note: See TracBrowser for help on using the repository browser.