source: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/inter_barxy_m.F90 @ 5136

Last change on this file since 5136 was 5136, checked in by abarral, 8 weeks ago

Put comgeom.h, comgeom2.h into modules

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 12.0 KB
Line 
1! $Id$
2
3module inter_barxy_m
4
5  ! Authors: Robert SADOURNY, Phu LE VAN, Lionel GUEZ
6
7  IMPLICIT NONE
8
9  PRIVATE
10  public inter_barxy
11
12CONTAINS
13
14  SUBROUTINE inter_barxy(dlonid, dlatid, champ, rlonimod, rlatimod, champint)
15
16    USE lmdz_assert_eq, ONLY: assert_eq
17    USE lmdz_assert, ONLY: assert
18    USE lmdz_comgeom2, ONLY: aire, apoln, apols
19
20    INCLUDE "dimensions.h"
21    ! (for "iim", "jjm")
22
23    INCLUDE "paramet.h"
24    ! (for other included files)
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 modele (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 exprimees en degres)
115
116    USE lmdz_assert_eq, ONLY: assert_eq
117    USE lmdz_libmath, ONLY: minmax
118
119    IMPLICIT NONE
120
121    REAL, INTENT(IN) :: dlonid(:)
122    REAL, INTENT(IN) :: fdat(:)
123    REAL, INTENT(IN) :: rlonimod(:)
124
125    REAL inter_barx(size(rlonimod))
126
127    !    ...  Variables locales ...
128
129    INTEGER idatmax, imodmax
130    REAL xxid(size(dlonid) + 1), xxd(size(dlonid) + 1), fdd(size(dlonid) + 1)
131    REAL  fxd(size(dlonid) + 1), xchan(size(dlonid) + 1), fdchan(size(dlonid) + 1)
132    REAL  xxim(size(rlonimod))
133
134    REAL x0, xim0, dx, dxm
135    REAL chmin, chmax, pi
136
137    INTEGER imod, idat, i, ichang, id0, id1, nid, idatmax1
138
139    !-----------------------------------------------------
140
141    idatmax = assert_eq(size(dlonid), size(fdat), "inter_barx idatmax")
142    imodmax = size(rlonimod)
143
144    pi = 2. * ASIN(1.)
145
146    !   REDEFINITION DE L'ORIGINE DES ABSCISSES
147    !    A L'INTERFACE OUEST DE LA PREMIERE MAILLE DU MODELE 
148    DO imod = 1, imodmax
149      xxim(imod) = rlonimod(imod)
150    ENDDO
151
152    CALL minmax(imodmax, xxim, chmin, chmax)
153    IF(chmax<6.50)   THEN
154      DO imod = 1, imodmax
155        xxim(imod) = xxim(imod) * 180. / pi
156      ENDDO
157    ENDIF
158
159    xim0 = xxim(imodmax) - 360.
160
161    DO imod = 1, imodmax
162      xxim(imod) = xxim(imod) - xim0
163    ENDDO
164
165    idatmax1 = idatmax + 1
166
167    DO idat = 1, idatmax
168      xxd(idat) = dlonid(idat)
169    ENDDO
170
171    CALL minmax(idatmax, xxd, chmin, chmax)
172    IF(chmax<6.50)  THEN
173      DO idat = 1, idatmax
174        xxd(idat) = xxd(idat) * 180. / pi
175      ENDDO
176    ENDIF
177
178    DO idat = 1, idatmax
179      xxd(idat) = MOD(xxd(idat) - xim0, 360.)
180      fdd(idat) = fdat (idat)
181    ENDDO
182
183    i = 2
184    DO while (xxd(i) >= xxd(i - 1) .AND. i < idatmax)
185      i = i + 1
186    ENDDO
187    IF (xxd(i) < xxd(i - 1)) THEN
188      ichang = i
189      !  ***  reorganisation  des longitudes entre 0. et 360. degres ****
190      nid = idatmax - ichang + 1
191      DO i = 1, nid
192        xchan (i) = xxd(i + ichang - 1)
193        fdchan(i) = fdd(i + ichang - 1)
194      ENDDO
195      DO i = 1, ichang - 1
196        xchan (i + nid) = xxd(i)
197        fdchan(i + nid) = fdd(i)
198      ENDDO
199      DO i = 1, idatmax
200        xxd(i) = xchan(i)
201        fdd(i) = fdchan(i)
202      ENDDO
203    end IF
204
205    !    translation des champs de donnees par rapport
206    !    a la nouvelle origine, avec redondance de la
207    !       maille a cheval sur les bords
208
209    id0 = 0
210    id1 = 0
211
212    DO idat = 1, idatmax
213      IF (xxd(idatmax1 - idat)<360.) exit
214      id1 = id1 + 1
215    ENDDO
216
217    DO idat = 1, idatmax
218      IF (xxd(idat)>0.) exit
219      id0 = id0 + 1
220    END DO
221
222    IF(id1 /= 0) THEN
223      DO idat = 1, id1
224        xxid(idat) = xxd(idatmax - id1 + idat) - 360.
225        fxd (idat) = fdd(idatmax - id1 + idat)
226      END DO
227      DO idat = 1, idatmax - id1
228        xxid(idat + id1) = xxd(idat)
229        fxd (idat + id1) = fdd(idat)
230      END DO
231    end IF
232
233    IF(id0 /= 0) THEN
234      DO idat = 1, idatmax - id0
235        xxid(idat) = xxd(idat + id0)
236        fxd (idat) = fdd(idat + id0)
237      END DO
238
239      DO idat = 1, id0
240        xxid (idatmax - id0 + idat) = xxd(idat) + 360.
241        fxd  (idatmax - id0 + idat) = fdd(idat)
242      END DO
243    else
244      DO idat = 1, idatmax
245        xxid(idat) = xxd(idat)
246        fxd (idat) = fdd(idat)
247      ENDDO
248    end IF
249    xxid(idatmax1) = xxid(1) + 360.
250    fxd (idatmax1) = fxd(1)
251
252    !   initialisation du champ du modele
253
254    inter_barx(:) = 0.
255
256    ! iteration
257
258    x0 = xim0
259    dxm = 0.
260    imod = 1
261    idat = 1
262
263    do while (imod <= imodmax)
264      do while (xxim(imod)>xxid(idat))
265        dx = xxid(idat) - x0
266        dxm = dxm + dx
267        inter_barx(imod) = inter_barx(imod) + dx * fxd(idat)
268        x0 = xxid(idat)
269        idat = idat + 1
270      END DO
271      IF (xxim(imod)<xxid(idat)) THEN
272        dx = xxim(imod) - x0
273        dxm = dxm + dx
274        inter_barx(imod) = (inter_barx(imod) + dx * fxd(idat)) / dxm
275        x0 = xxim(imod)
276        dxm = 0.
277        imod = imod + 1
278      ELSE
279        dx = xxim(imod) - x0
280        dxm = dxm + dx
281        inter_barx(imod) = (inter_barx(imod) + dx * fxd(idat)) / dxm
282        x0 = xxim(imod)
283        dxm = 0.
284        imod = imod + 1
285        idat = idat + 1
286      END IF
287    END DO
288
289  END function inter_barx
290
291  !******************************
292
293  function inter_bary(yjdat, fdat, yjmod)
294
295    ! Interpolation barycentrique basée sur les aires.
296    ! Version unidimensionnelle, en latitude.
297    ! L'indice 1 correspond à l'interface maille 1 -- maille 2.
298
299    USE lmdz_assert, ONLY: assert
300
301    IMPLICIT NONE
302
303    REAL, INTENT(IN) :: yjdat(:)
304    ! (angles, ordonnées des interfaces des mailles des données, in
305    ! degrees, in increasing order)
306
307    REAL, INTENT(IN) :: fdat(:) ! champ de données
308
309    REAL, INTENT(IN) :: yjmod(:)
310    ! (ordonnées des interfaces des mailles du modèle)
311    ! (in degrees, in strictly increasing order)
312
313    REAL inter_bary(size(yjmod)) ! champ du modèle
314
315    ! Variables local to the procedure:
316
317    REAL y0, dy, dym
318    INTEGER jdat ! indice du champ de données
319    INTEGER jmod ! indice du champ du modèle
320
321    !------------------------------------
322
323    CALL assert(size(yjdat) == size(fdat), "inter_bary")
324
325    ! Initialisation des variables
326    inter_bary(:) = 0.
327    y0 = -90.
328    dym = 0.
329    jmod = 1
330    jdat = 1
331
332    do while (jmod <= size(yjmod))
333      do while (yjmod(jmod) > yjdat(jdat))
334        dy = yjdat(jdat) - y0
335        dym = dym + dy
336        inter_bary(jmod) = inter_bary(jmod) + dy * fdat(jdat)
337        y0 = yjdat(jdat)
338        jdat = jdat + 1
339      END DO
340      IF (yjmod(jmod) < yjdat(jdat)) THEN
341        dy = yjmod(jmod) - y0
342        dym = dym + dy
343        inter_bary(jmod) = (inter_bary(jmod) + dy * fdat(jdat)) / dym
344        y0 = yjmod(jmod)
345        dym = 0.
346        jmod = jmod + 1
347      ELSE
348        ! {yjmod(jmod) == yjdat(jdat)}
349        dy = yjmod(jmod) - y0
350        dym = dym + dy
351        inter_bary(jmod) = (inter_bary(jmod) + dy * fdat(jdat)) / dym
352        y0 = yjmod(jmod)
353        dym = 0.
354        jmod = jmod + 1
355        jdat = jdat + 1
356      END IF
357    END DO
358    ! Le test de fin suppose que l'interface 0 est commune aux deux
359    ! grilles "yjdat" et "yjmod".
360
361  END function inter_bary
362
363  !******************************
364
365  SUBROUTINE ord_coord(xi, xo, decrois)
366
367    ! This procedure receives an array of latitudes.
368    ! It converts them to degrees if they are in radians.
369    ! If the input latitudes are in decreasing order, the procedure
370    ! reverses their order.
371    ! Finally, the procedure adds 90° as the last value of the array.
372
373    USE lmdz_assert_eq, ONLY: assert_eq
374    USE comconst_mod, ONLY: pi
375
376    IMPLICIT NONE
377
378    REAL, INTENT(IN) :: xi(:)
379    ! (latitude, in degrees or radians, in increasing or decreasing order)
380    ! ("xi" should contain latitudes from pole to pole.
381    ! "xi" should contain the latitudes of the boundaries of grid
382    ! cells, not the centers of grid cells.
383    ! So the extreme values should not be 90° and -90°.)
384
385    REAL, INTENT(OUT) :: xo(:) ! angles in degrees
386    LOGICAL, INTENT(OUT) :: decrois
387
388    ! Variables  local to the procedure:
389    INTEGER nmax, i
390
391    !--------------------
392
393    nmax = assert_eq(size(xi), size(xo) - 1, "ord_coord")
394
395    ! Check monotonicity:
396    decrois = xi(2) < xi(1)
397    DO i = 3, nmax
398      IF (decrois .neqv. xi(i) < xi(i - 1)) stop &
399              '"ord_coord":  latitudes are not monotonic'
400    ENDDO
401
402    IF (abs(xi(1)) < pi) THEN
403      ! "xi" contains latitudes in radians
404      xo(:nmax) = xi(:) * 180. / pi
405    else
406      ! "xi" contains latitudes in degrees
407      xo(:nmax) = xi(:)
408    end IF
409
410    IF (ABS(abs(xo(1)) - 90) < 0.001 .OR. ABS(abs(xo(nmax)) - 90) < 0.001) THEN
411      print *, "ord_coord"
412      PRINT *, '"xi" should contain the latitudes of the boundaries of ' &
413              // 'grid cells, not the centers of grid cells.'
414      STOP
415    ENDIF
416
417    IF (decrois) xo(:nmax) = xo(nmax:1:- 1)
418    xo(nmax + 1) = 90.
419
420  END SUBROUTINE ord_coord
421
422  !***********************************
423
424  function ord_coordm(xi)
425
426    ! This procedure converts to degrees, if necessary, and inverts the
427    ! order.
428
429    USE comconst_mod, ONLY: pi
430
431    IMPLICIT NONE
432
433    REAL, INTENT(IN) :: xi(:) ! angle, in rad or degrees
434    REAL ord_coordm(size(xi)) ! angle, in degrees
435
436    !-----------------------------
437
438    IF (xi(1) < 6.5) THEN
439      ! "xi" is in rad
440      ord_coordm(:) = xi(size(xi):1:-1) * 180. / pi
441    else
442      ! "xi" is in degrees
443      ord_coordm(:) = xi(size(xi):1:-1)
444    ENDIF
445
446  END function ord_coordm
447
448END MODULE inter_barxy_m
Note: See TracBrowser for help on using the repository browser.