source: LMDZ6/branches/Amaury_dev/libf/misc/regr_conserv_m.F90 @ 5114

Last change on this file since 5114 was 5113, checked in by abarral, 4 months ago

Rename modules in misc from *_mod > lmdz_*
Put cbrt.f90, ch*.f90, pch*.f90 in new lmdz_libmath_pch.f90
File size: 15.5 KB
Line 
1MODULE regr_conserv_m
2
3  USE lmdz_assert_eq,   ONLY: assert_eq
4  USE lmdz_assert,      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)/2.)
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,:)/2.)
145    IF(ix==2) vt(:,it) = vt(:,it)+idt*(b-a)*(vs(:,is)+co*slope(:,is)/2.)
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,:,:)/2.)
205    IF(ix==2) vt(:,it,:) = vt(:,it,:)+idt*(b-a)*(vs(:,is,:)+co*slope(:,is,:)/2.)
206    IF(ix==3) vt(:,:,it) = vt(:,:,it)+idt*(b-a)*(vs(:,:,is)+co*slope(:,:,is)/2.)
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,:,:,:)/2.)
267    IF(ix==2) vt(:,it,:,:) = vt(:,it,:,:)+idt*(b-a)*(vs(:,is,:,:)+co*slope(:,is,:,:)/2.)
268    IF(ix==3) vt(:,:,it,:) = vt(:,:,it,:)+idt*(b-a)*(vs(:,:,is,:)+co*slope(:,:,is,:)/2.)
269    IF(ix==4) vt(:,:,:,it) = vt(:,:,:,it)+idt*(b-a)*(vs(:,:,:,is)+co*slope(:,:,:,is)/2.)
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,:,:,:,:)/2.)
331    IF(ix==2) vt(:,it,:,:,:) = vt(:,it,:,:,:)+idt*(b-a)*(vs(:,is,:,:,:)+co*slope(:,is,:,:,:)/2.)
332    IF(ix==3) vt(:,:,it,:,:) = vt(:,:,it,:,:)+idt*(b-a)*(vs(:,:,is,:,:)+co*slope(:,:,is,:,:)/2.)
333    IF(ix==4) vt(:,:,:,it,:) = vt(:,:,:,it,:)+idt*(b-a)*(vs(:,:,:,is,:)+co*slope(:,:,:,is,:)/2.)
334    IF(ix==5) vt(:,:,:,:,it) = vt(:,:,:,:,it)+idt*(b-a)*(vs(:,:,:,:,is)+co*slope(:,:,:,:,is)/2.)
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.