source: trunk/LMDZ.COMMON/libf/dyn3d_common/inter_barxy_m.F90 @ 3552

Last change on this file since 3552 was 1422, checked in by milmd, 10 years ago

In GENERIC, MARS and COMMON models replace some include files by modules (usefull for decoupling physics with dynamics).
File size: 12.2 KB
Line 
1!
2! $Id$
3!
4module inter_barxy_m
5
6  ! Authors: Robert SADOURNY, Phu LE VAN, Lionel GUEZ
7
8  implicit none
9
10  private
11  public inter_barxy
12
13contains
14
15  SUBROUTINE inter_barxy(dlonid, dlatid, champ, rlonimod, rlatimod, champint)
16
17    use assert_eq_m, only: assert_eq
18    use assert_m, only: assert
19
20    include "dimensions.h"
21    ! (for "iim", "jjm")
22
23    include "paramet.h"
24    ! (for other included files)
25
26    include "comgeom2.h"
27    ! (for "aire", "apoln", "apols")
28
29    REAL, intent(in):: dlonid(:)
30    ! (longitude from input file, in rad, from -pi to pi)
31
32    REAL, intent(in):: dlatid(:), champ(:, :), rlonimod(:)
33
34    REAL, intent(in):: rlatimod(:)
35    ! (latitude angle, in degrees or rad, in strictly decreasing order)
36
37    real, intent(out):: champint(:, :)
38    ! Si taille de la seconde dim = jjm + 1, on veut interpoler sur les
39    ! jjm+1 latitudes rlatu du modele (latitudes des scalaires et de U)
40    ! Si taille de la seconde dim = jjm, on veut interpoler sur les
41    ! jjm latitudes rlatv du modele (latitudes de V)
42
43    ! Variables local to the procedure:
44
45    REAL champy(iim, size(champ, 2))
46    integer j, i, jnterfd, jmods
47
48    REAL yjmod(size(champint, 2))
49    ! (angle, in degrees, in strictly increasing order)
50
51    REAL   yjdat(size(dlatid) + 1) ! angle, in degrees, in increasing order
52    LOGICAL decrois ! "dlatid" is in decreasing order
53
54    !-----------------------------------
55
56    jnterfd = assert_eq(size(champ, 2) - 1, size(dlatid), &
57         "inter_barxy jnterfd")
58    jmods = size(champint, 2)
59    call assert(size(champ, 1) == size(dlonid), "inter_barxy size(champ, 1)")
60    call assert((/size(rlonimod), size(champint, 1)/) == iim, &
61         "inter_barxy iim")
62    call assert(any(jmods == (/jjm, jjm + 1/)), 'inter_barxy jmods')
63    call assert(size(rlatimod) == jjm, "inter_barxy size(rlatimod)")
64
65    ! Check decreasing order for "rlatimod":
66    DO i = 2, jjm
67       IF (rlatimod(i) >= rlatimod(i-1)) stop &
68            '"inter_barxy": "rlatimod" should be strictly decreasing'
69    ENDDO
70
71    yjmod(:jjm) = ord_coordm(rlatimod)
72    IF (jmods == jjm + 1) THEN
73       IF (90. - yjmod(jjm) < 0.01) stop &
74            '"inter_barxy": with jmods = jjm + 1, yjmod(jjm) should be < 90.'
75    ELSE
76       ! jmods = jjm
77       IF (ABS(yjmod(jjm) - 90.) > 0.01) stop &
78            '"inter_barxy": with jmods = jjm, yjmod(jjm) should be 90.'
79    ENDIF
80
81    if (jmods == jjm + 1) yjmod(jjm + 1) = 90.
82
83    DO j = 1, jnterfd + 1
84       champy(:, j) = inter_barx(dlonid, champ(:, j), rlonimod)
85    ENDDO
86
87    CALL ord_coord(dlatid, yjdat, decrois)
88    IF (decrois) champy(:, :) = champy(:, jnterfd + 1:1:-1)
89    DO i = 1, iim
90       champint(i, :) = inter_bary(yjdat, champy(i, :), yjmod)
91    ENDDO
92    champint(:, :) = champint(:, jmods:1:-1)
93
94    IF (jmods == jjm + 1) THEN
95       ! Valeurs uniques aux poles
96       champint(:, 1) = SUM(aire(:iim,  1) * champint(:, 1)) / apoln
97       champint(:, jjm + 1) = SUM(aire(:iim, jjm + 1) &
98            * champint(:, jjm + 1)) / apols
99    ENDIF
100
101  END SUBROUTINE inter_barxy
102
103  !******************************
104
105  function inter_barx(dlonid, fdat, rlonimod)
106
107    !        INTERPOLATION BARYCENTRIQUE BASEE SUR LES AIRES
108    !            VERSION UNIDIMENSIONNELLE  ,   EN  LONGITUDE .
109
110    !     idat : indice du champ de donnees, de 1 a idatmax
111    !     imod : indice du champ du modele,  de 1 a  imodmax
112    !     fdat(idat) : champ de donnees (entrees)
113    !     inter_barx(imod) : champ du modele (sorties)
114    !     dlonid(idat): abscisses des interfaces des mailles donnees
115    !     rlonimod(imod): abscisses des interfaces des mailles modele
116    !      ( L'indice 1 correspond a l'interface mailLE 1 / maille 2)
117    !      ( Les abscisses sont exprimees en degres)
118
119    use assert_eq_m, only: assert_eq
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.LT.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.LT.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 ).LT.360.) exit
216       id1 = id1 + 1
217    ENDDO
218
219    DO idat = 1, idatmax
220       IF (xxd(idat).GT.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).GT.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).LT.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 assert_m, 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 assert_eq_m, only: assert_eq
376    USE comconst_mod, ONLY: pi
377
378    IMPLICIT NONE
379
380    ! (for "pi")
381
382    REAL, intent(in):: xi(:)
383    ! (latitude, in degrees or radians, in increasing or decreasing order)
384    ! ("xi" should contain latitudes from pole to pole.
385    ! "xi" should contain the latitudes of the boundaries of grid
386    ! cells, not the centers of grid cells.
387    ! So the extreme values should not be 90° and -90°.)
388
389    REAL, intent(out):: xo(:) ! angles in degrees
390    LOGICAL, intent(out):: decrois
391
392    ! Variables  local to the procedure:
393    INTEGER nmax, i
394
395    !--------------------
396
397    nmax = assert_eq(size(xi), size(xo) - 1, "ord_coord")
398
399    ! Check monotonicity:
400    decrois = xi(2) < xi(1)
401    DO i = 3, nmax
402       IF (decrois .neqv. xi(i) < xi(i-1)) stop &
403            '"ord_coord":  latitudes are not monotonic'
404    ENDDO
405
406    IF (abs(xi(1)) < pi) then
407       ! "xi" contains latitudes in radians
408       xo(:nmax) = xi(:) * 180. / pi
409    else
410       ! "xi" contains latitudes in degrees
411       xo(:nmax) = xi(:)
412    end IF
413
414    IF (ABS(abs(xo(1)) - 90) < 0.001 .or. ABS(abs(xo(nmax)) - 90) < 0.001) THEN
415       print *, "ord_coord"
416       PRINT *, '"xi" should contain the latitudes of the boundaries of ' &
417            // 'grid cells, not the centers of grid cells.'
418       STOP
419    ENDIF
420
421    IF (decrois) xo(:nmax) = xo(nmax:1:- 1)
422    xo(nmax + 1) = 90.
423
424  END SUBROUTINE ord_coord
425
426  !***********************************
427
428  function ord_coordm(xi)
429
430    ! This procedure converts to degrees, if necessary, and inverts the
431    ! order.
432
433    USE comconst_mod, ONLY: pi
434
435    IMPLICIT NONE
436
437    ! (for "pi")
438
439    REAL, intent(in):: xi(:) ! angle, in rad or degrees
440    REAL ord_coordm(size(xi)) ! angle, in degrees
441
442    !-----------------------------
443
444    IF (xi(1) < 6.5) THEN
445       ! "xi" is in rad
446       ord_coordm(:) = xi(size(xi):1:-1) * 180. / pi
447    else
448       ! "xi" is in degrees
449       ord_coordm(:) = xi(size(xi):1:-1)
450    ENDIF
451
452  END function ord_coordm
453
454end module inter_barxy_m
Note: See TracBrowser for help on using the repository browser.