source: LMDZ5/trunk/libf/inter_barxy_m.F90 @ 1630

Last change on this file since 1630 was 1630, checked in by Laurent Fairhead, 12 years ago

Importation initiale du répertoire dyn3dmem


Initial import of dyn3dmem directory

File size: 12.2 KB
RevLine 
[1630]1module inter_barxy_m
2
3  ! Authors: Robert SADOURNY, Phu LE VAN, Lionel GUEZ
4
5  implicit none
6
7  private
8  public inter_barxy
9
10contains
11
12  SUBROUTINE inter_barxy(dlonid, dlatid, champ, rlonimod, rlatimod, champint)
13
14    use assert_eq_m, only: assert_eq
15    use assert_m, only: assert
16
17    include "dimensions.h"
18    ! (for "iim", "jjm")
19
20    include "paramet.h"
21    ! (for other included files)
22
23    include "comgeom2.h"
24    ! (for "aire", "apoln", "apols")
25
26    REAL, intent(in):: dlonid(:)
27    ! (longitude from input file, in rad, from -pi to pi)
28
29    REAL, intent(in):: dlatid(:), champ(:, :), rlonimod(:)
30
31    REAL, intent(in):: rlatimod(:)
32    ! (latitude angle, in degrees or rad, in strictly decreasing order)
33
34    real, intent(out):: champint(:, :)
35    ! Si taille de la seconde dim = jjm + 1, on veut interpoler sur les
36    ! jjm+1 latitudes rlatu du modele (latitudes des scalaires et de U)
37    ! Si taille de la seconde dim = jjm, on veut interpoler sur les
38    ! jjm latitudes rlatv du modèle (latitudes de V)
39
40    ! Variables local to the procedure:
41
42    REAL champy(iim, size(champ, 2))
43    integer j, i, jnterfd, jmods
44
45    REAL yjmod(size(champint, 2))
46    ! (angle, in degrees, in strictly increasing order)
47
48    REAL   yjdat(size(dlatid) + 1) ! angle, in degrees, in increasing order
49    LOGICAL decrois ! "dlatid" is in decreasing order
50
51    !-----------------------------------
52
53    jnterfd = assert_eq(size(champ, 2) - 1, size(dlatid), &
54         "inter_barxy jnterfd")
55    jmods = size(champint, 2)
56    call assert(size(champ, 1) == size(dlonid), "inter_barxy size(champ, 1)")
57    call assert((/size(rlonimod), size(champint, 1)/) == iim, &
58         "inter_barxy iim")
59    call assert(any(jmods == (/jjm, jjm + 1/)), 'inter_barxy jmods')
60    call assert(size(rlatimod) == jjm, "inter_barxy size(rlatimod)")
61
62    ! Check decreasing order for "rlatimod":
63    DO i = 2, jjm
64       IF (rlatimod(i) >= rlatimod(i-1)) stop &
65            '"inter_barxy": "rlatimod" should be strictly decreasing'
66    ENDDO
67
68    yjmod(:jjm) = ord_coordm(rlatimod)
69    IF (jmods == jjm + 1) THEN
70       IF (90. - yjmod(jjm) < 0.01) stop &
71            '"inter_barxy": with jmods = jjm + 1, yjmod(jjm) should be < 90.'
72    ELSE
73       ! jmods = jjm
74       IF (ABS(yjmod(jjm) - 90.) > 0.01) stop &
75            '"inter_barxy": with jmods = jjm, yjmod(jjm) should be 90.'
76    ENDIF
77
78    if (jmods == jjm + 1) yjmod(jjm + 1) = 90.
79
80    DO j = 1, jnterfd + 1
81       champy(:, j) = inter_barx(dlonid, champ(:, j), rlonimod)
82    ENDDO
83
84    CALL ord_coord(dlatid, yjdat, decrois)
85    IF (decrois) champy(:, :) = champy(:, jnterfd + 1:1:-1)
86    DO i = 1, iim
87       champint(i, :) = inter_bary(yjdat, champy(i, :), yjmod)
88    ENDDO
89    champint(:, :) = champint(:, jmods:1:-1)
90
91    IF (jmods == jjm + 1) THEN
92       ! Valeurs uniques aux poles
93       champint(:, 1) = SUM(aire(:iim,  1) * champint(:, 1)) / apoln
94       champint(:, jjm + 1) = SUM(aire(:iim, jjm + 1) &
95            * champint(:, jjm + 1)) / apols
96    ENDIF
97
98  END SUBROUTINE inter_barxy
99
100  !******************************
101
102  function inter_barx(dlonid, fdat, rlonimod)
103
104    !        INTERPOLATION BARYCENTRIQUE BASEE SUR LES AIRES
105    !            VERSION UNIDIMENSIONNELLE  ,   EN  LONGITUDE .
106
107    !     idat : indice du champ de donnees, de 1 a idatmax
108    !     imod : indice du champ du modele,  de 1 a  imodmax
109    !     fdat(idat) : champ de donnees (entrees)
110    !     inter_barx(imod) : champ du modele (sorties)
111    !     dlonid(idat): abscisses des interfaces des mailles donnees
112    !     rlonimod(imod): abscisses des interfaces des mailles modele
113    !      ( L'indice 1 correspond a l'interface mailLE 1 / maille 2)
114    !      ( Les abscisses sont exprimées en degres)
115
116    use assert_eq_m, only: assert_eq
117
118    IMPLICIT NONE
119
120    REAL, intent(in):: dlonid(:)
121    real, intent(in):: fdat(:)
122    real, intent(in):: rlonimod(:)
123
124    real inter_barx(size(rlonimod))
125
126    !    ...  Variables locales ...
127
128    INTEGER idatmax, imodmax
129    REAL xxid(size(dlonid)+1), xxd(size(dlonid)+1), fdd(size(dlonid)+1)
130    REAL  fxd(size(dlonid)+1), xchan(size(dlonid)+1), fdchan(size(dlonid)+1)
131    REAL  xxim(size(rlonimod))
132
133    REAL x0, xim0, dx, dxm
134    REAL chmin, chmax, pi
135
136    INTEGER imod, idat, i, ichang, id0, id1, nid, idatmax1
137
138    !-----------------------------------------------------
139
140    idatmax = assert_eq(size(dlonid), size(fdat), "inter_barx idatmax")
141    imodmax = size(rlonimod)
142
143    pi = 2. * ASIN(1.)
144
145    !   REDEFINITION DE L'ORIGINE DES ABSCISSES
146    !    A L'INTERFACE OUEST DE LA PREMIERE MAILLE DU MODELE 
147    DO imod = 1, imodmax
148       xxim(imod) = rlonimod(imod)
149    ENDDO
150
151    CALL minmax( imodmax, xxim, chmin, chmax)
152    IF( chmax.LT.6.50 )   THEN
153       DO imod = 1, imodmax
154          xxim(imod) = xxim(imod) * 180./pi
155       ENDDO
156    ENDIF
157
158    xim0 = xxim(imodmax) - 360.
159
160    DO imod = 1, imodmax
161       xxim(imod) = xxim(imod) - xim0
162    ENDDO
163
164    idatmax1 = idatmax +1
165
166    DO idat = 1, idatmax
167       xxd(idat) = dlonid(idat)
168    ENDDO
169
170    CALL minmax( idatmax, xxd, chmin, chmax)
171    IF( chmax.LT.6.50 )  THEN
172       DO idat = 1, idatmax
173          xxd(idat) = xxd(idat) * 180./pi
174       ENDDO
175    ENDIF
176
177    DO idat = 1, idatmax
178       xxd(idat) = AMOD( xxd(idat) - xim0, 360. )
179       fdd(idat) = fdat (idat)
180    ENDDO
181
182    i = 2
183    DO while (xxd(i) >= xxd(i-1) .and. i < idatmax)
184       i = i + 1
185    ENDDO
186    IF (xxd(i) < xxd(i-1)) THEN
187       ichang = i
188       !  ***  reorganisation  des longitudes entre 0. et 360. degres ****
189       nid = idatmax - ichang +1
190       DO i = 1, nid
191          xchan (i) = xxd(i+ichang -1 )
192          fdchan(i) = fdd(i+ichang -1 )
193       ENDDO
194       DO i=1, ichang -1
195          xchan (i+ nid) = xxd(i)
196          fdchan(i+nid) = fdd(i)
197       ENDDO
198       DO i =1, idatmax
199          xxd(i) = xchan(i)
200          fdd(i) = fdchan(i)
201       ENDDO
202    end IF
203
204    !    translation des champs de donnees par rapport
205    !    a la nouvelle origine, avec redondance de la
206    !       maille a cheval sur les bords
207
208    id0 = 0
209    id1 = 0
210
211    DO idat = 1, idatmax
212       IF ( xxd( idatmax1- idat ).LT.360.) exit
213       id1 = id1 + 1
214    ENDDO
215
216    DO idat = 1, idatmax
217       IF (xxd(idat).GT.0.) exit
218       id0 = id0 + 1
219    END DO
220
221    IF( id1 /= 0 ) then
222       DO idat = 1, id1
223          xxid(idat) = xxd(idatmax - id1 + idat) - 360.
224          fxd (idat) = fdd(idatmax - id1 + idat)     
225       END DO
226       DO idat = 1, idatmax - id1
227          xxid(idat + id1) = xxd(idat)
228          fxd (idat + id1) = fdd(idat)
229       END DO
230    end IF
231
232    IF(id0 /= 0) then
233       DO idat = 1, idatmax - id0
234          xxid(idat) = xxd(idat + id0)
235          fxd (idat) = fdd(idat + id0)
236       END DO
237
238       DO idat = 1, id0
239          xxid (idatmax - id0 + idat) =  xxd(idat) + 360.
240          fxd  (idatmax - id0 + idat) =  fdd(idat)   
241       END DO
242    else
243       DO idat = 1, idatmax
244          xxid(idat)  = xxd(idat)
245          fxd (idat)  = fdd(idat)
246       ENDDO
247    end IF
248    xxid(idatmax1) = xxid(1) + 360.
249    fxd (idatmax1) = fxd(1)
250
251    !   initialisation du champ du modele
252
253    inter_barx(:) = 0.
254
255    ! iteration
256
257    x0   = xim0
258    dxm  = 0.
259    imod = 1
260    idat = 1
261
262    do while (imod <= imodmax)
263       do while (xxim(imod).GT.xxid(idat))
264          dx   = xxid(idat) - x0
265          dxm  = dxm + dx
266          inter_barx(imod) = inter_barx(imod) + dx * fxd(idat)
267          x0   = xxid(idat)
268          idat = idat + 1
269       end do
270       IF (xxim(imod).LT.xxid(idat)) THEN
271          dx   = xxim(imod) - x0
272          dxm  = dxm + dx
273          inter_barx(imod) = (inter_barx(imod) + dx * fxd(idat)) / dxm
274          x0   = xxim(imod)
275          dxm  = 0.
276          imod = imod + 1
277       ELSE
278          dx   = xxim(imod) - x0
279          dxm  = dxm + dx
280          inter_barx(imod) = (inter_barx(imod) + dx * fxd(idat)) / dxm
281          x0   = xxim(imod)
282          dxm  = 0.
283          imod = imod + 1
284          idat = idat + 1
285       END IF
286    end do
287
288  END function inter_barx
289
290  !******************************
291
292  function inter_bary(yjdat, fdat, yjmod)
293
294    ! Interpolation barycentrique basée sur les aires.
295    ! Version unidimensionnelle, en latitude.
296    ! L'indice 1 correspond à l'interface maille 1 -- maille 2.
297
298    use assert_m, only: assert
299
300    IMPLICIT NONE
301
302    REAL, intent(in):: yjdat(:)
303    ! (angles, ordonnées des interfaces des mailles des données, in
304    ! degrees, in increasing order)
305
306    REAL, intent(in):: fdat(:) ! champ de données
307
308    REAL, intent(in):: yjmod(:)
309    ! (ordonnées des interfaces des mailles du modèle)
310    ! (in degrees, in strictly increasing order)
311
312    REAL inter_bary(size(yjmod)) ! champ du modèle
313
314    ! Variables local to the procedure:
315
316    REAL y0, dy, dym
317    INTEGER jdat ! indice du champ de données
318    integer jmod ! indice du champ du modèle
319
320    !------------------------------------
321
322    call assert(size(yjdat) == size(fdat), "inter_bary")
323
324    ! Initialisation des variables
325    inter_bary(:) = 0.
326    y0    = -90.
327    dym   = 0.
328    jmod  = 1
329    jdat  = 1
330
331    do while (jmod <= size(yjmod))
332       do while (yjmod(jmod) > yjdat(jdat))
333          dy         = yjdat(jdat) - y0
334          dym        = dym + dy
335          inter_bary(jmod) = inter_bary(jmod) + dy * fdat(jdat)
336          y0         = yjdat(jdat)
337          jdat       = jdat + 1
338       end do
339       IF (yjmod(jmod) < yjdat(jdat)) THEN
340          dy         = yjmod(jmod) - y0
341          dym        = dym + dy
342          inter_bary(jmod) = (inter_bary(jmod) + dy * fdat(jdat)) / dym
343          y0         = yjmod(jmod)
344          dym        = 0.
345          jmod       = jmod + 1
346       ELSE
347          ! {yjmod(jmod) == yjdat(jdat)}
348          dy         = yjmod(jmod) - y0
349          dym        = dym + dy
350          inter_bary(jmod) = (inter_bary(jmod) + dy * fdat(jdat)) / dym
351          y0         = yjmod(jmod)
352          dym        = 0.
353          jmod       = jmod + 1
354          jdat       = jdat + 1
355       END IF
356    end do
357    ! Le test de fin suppose que l'interface 0 est commune aux deux
358    ! grilles "yjdat" et "yjmod".
359
360  END function inter_bary
361
362  !******************************
363
364  SUBROUTINE ord_coord(xi, xo, decrois)
365
366    ! This procedure receives an array of latitudes.
367    ! It converts them to degrees if they are in radians.
368    ! If the input latitudes are in decreasing order, the procedure
369    ! reverses their order.
370    ! Finally, the procedure adds 90° as the last value of the array.
371
372    use assert_eq_m, only: assert_eq
373
374    IMPLICIT NONE
375
376    include "comconst.h"
377    ! (for "pi")
378
379    REAL, intent(in):: xi(:)
380    ! (latitude, in degrees or radians, in increasing or decreasing order)
381    ! ("xi" should contain latitudes from pole to pole.
382    ! "xi" should contain the latitudes of the boundaries of grid
383    ! cells, not the centers of grid cells.
384    ! So the extreme values should not be 90° and -90°.)
385
386    REAL, intent(out):: xo(:) ! angles in degrees
387    LOGICAL, intent(out):: decrois
388
389    ! Variables  local to the procedure:
390    INTEGER nmax, i
391
392    !--------------------
393
394    nmax = assert_eq(size(xi), size(xo) - 1, "ord_coord")
395
396    ! Check monotonicity:
397    decrois = xi(2) < xi(1)
398    DO i = 3, nmax
399       IF (decrois .neqv. xi(i) < xi(i-1)) stop &
400            '"ord_coord":  latitudes are not monotonic'
401    ENDDO
402
403    IF (abs(xi(1)) < pi) then
404       ! "xi" contains latitudes in radians
405       xo(:nmax) = xi(:) * 180. / pi
406    else
407       ! "xi" contains latitudes in degrees
408       xo(:nmax) = xi(:)
409    end IF
410
411    IF (ABS(abs(xo(1)) - 90) < 0.001 .or. ABS(abs(xo(nmax)) - 90) < 0.001) THEN
412       print *, "ord_coord"
413       PRINT *, '"xi" should contain the latitudes of the boundaries of ' &
414            // 'grid cells, not the centers of grid cells.'
415       STOP
416    ENDIF
417
418    IF (decrois) xo(:nmax) = xo(nmax:1:- 1)
419    xo(nmax + 1) = 90.
420
421  END SUBROUTINE ord_coord
422
423  !***********************************
424
425  function ord_coordm(xi)
426
427    ! This procedure converts to degrees, if necessary, and inverts the
428    ! order.
429
430    IMPLICIT NONE
431
432    include "comconst.h"
433    ! (for "pi")
434
435    REAL, intent(in):: xi(:) ! angle, in rad or degrees
436    REAL ord_coordm(size(xi)) ! angle, in degrees
437
438    !-----------------------------
439
440    IF (xi(1) < 6.5) THEN
441       ! "xi" is in rad
442       ord_coordm(:) = xi(size(xi):1:-1) * 180. / pi
443    else
444       ! "xi" is in degrees
445       ord_coordm(:) = xi(size(xi):1:-1)
446    ENDIF
447
448  END function ord_coordm
449
450end module inter_barxy_m
Note: See TracBrowser for help on using the repository browser.