source: LMDZ6/branches/Amaury_dev/libf/misc/lmdz_regr_conserv.f90 @ 5128

Last change on this file since 5128 was 5119, checked in by abarral, 2 months ago

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