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

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