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

Last change on this file since 5209 was 5160, checked in by abarral, 7 weeks ago

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