trunk/dyn3d/inter_barxy.f

module inter_barxy_m

  ! From inter_barxy.F, version 1.1.1.1 2004/05/19 12:53:07

  implicit none

  private
  public inter_barxy

contains

  SUBROUTINE inter_barxy(dlonid, dlatid, champ, rlonimod, rlatimod, champint)

    ! Author: P. Le Van

    use nr_util, only: assert_eq, assert
    use dimens_m, only: iim, jjm
    use comgeom, only: aire_2d, apoln, apols

    REAL, intent(in):: dlonid(:)
    ! (longitude from input file, in rad, from -pi to pi)

    REAL, intent(in):: dlatid(:), champ(:, :), rlonimod(:)

    REAL, intent(in):: rlatimod(:)
    ! (latitude angle, in degrees or rad, in strictly decreasing order)

    real, intent(out):: champint(:, :)
    ! Si taille de la seconde dim = jjm + 1, on veut interpoler sur les
    ! jjm+1 latitudes rlatu du modele (latitudes des scalaires et de U)
    ! Si taille de la seconde dim = jjm, on veut interpoler sur les
    ! jjm latitudes rlatv du modèle (latitudes de V) 

    ! Variables local to the procedure:

    REAL champy(iim, size(champ, 2))
    integer j, i, jnterfd, jmods

    REAL yjmod(size(champint, 2))
    ! (angle, in degrees, in strictly increasing order)

    REAL   yjdat(size(dlatid) + 1) ! angle, in degrees, in increasing order
    LOGICAL decrois ! "dlatid" is in decreasing order

    !-----------------------------------

    jnterfd = assert_eq(size(champ, 2) - 1, size(dlatid), &
         "inter_barxy jnterfd")
    jmods = size(champint, 2)
    call assert(size(champ, 1) == size(dlonid), "inter_barxy size(champ, 1)")
    call assert((/size(rlonimod), size(champint, 1)/) == iim, &
         "inter_barxy iim")
    call assert(any(jmods == (/jjm, jjm + 1/)), 'inter_barxy jmods')
    call assert(size(rlatimod) == jjm, "inter_barxy size(rlatimod)")

    ! Check decreasing order for "rlatimod":
    DO i = 2, jjm
       IF (rlatimod(i) >= rlatimod(i-1)) then
          print *, '"inter_barxy": "rlatimod" should be strictly decreasing'
          stop 1
       end IF
    ENDDO

    yjmod(:jjm) = ord_coordm(rlatimod)
    IF (jmods == jjm + 1) THEN
       IF (90. - yjmod(jjm) < 0.01) then
          print *, '"inter_barxy": with jmods = jjm + 1, ' &
               // 'yjmod(jjm) should be < 90.'
          stop 1
       end IF
    ELSE
       ! jmods = jjm
       IF (ABS(yjmod(jjm) - 90.) > 0.01) then
          print *, '"inter_barxy": with jmods = jjm, yjmod(jjm) should be 90.'
          stop 1
       end IF
    ENDIF

    if (jmods == jjm + 1) yjmod(jjm + 1) = 90.

    DO j = 1, jnterfd + 1
       champy(:, j) = inter_barx(dlonid, champ(:, j), rlonimod)
    ENDDO

    CALL ord_coord(dlatid, yjdat, decrois) 
    IF (decrois) champy(:, :) = champy(:, jnterfd + 1:1:-1)
    DO i = 1, iim
       champint(i, :) = inter_bary(yjdat, champy(i, :), yjmod)
    ENDDO
    champint(:, :) = champint(:, jmods:1:-1)

    IF (jmods == jjm + 1) THEN
       ! Valeurs uniques aux poles
       champint(:, 1) = SUM(aire_2d(:iim,  1) * champint(:, 1)) / apoln
       champint(:, jjm + 1) = SUM(aire_2d(:iim, jjm + 1) &
            * champint(:, jjm + 1)) / apols
    ENDIF

  END SUBROUTINE inter_barxy

  !******************************

  function inter_barx(dlonid, fdat, rlonimod) 

    ! From dyn3d/inter_barx.F, v 1.1.1.1 2004/05/19 12:53:06

    ! Auteurs : Robert Sadourny, P. Le Van

    !        INTERPOLATION BARYCENTRIQUE BASEE SUR LES AIRES
    !            VERSION UNIDIMENSIONNELLE  ,   EN  LONGITUDE .

    !     idat : indice du champ de donnees, de 1 a idatmax
    !     imod : indice du champ du modele,  de 1 a  imodmax
    !     fdat(idat) : champ de donnees (entrees)
    !     inter_barx(imod) : champ du modele (sorties)
    !     dlonid(idat): abscisses des interfaces des mailles donnees
    !     rlonimod(imod): abscisses des interfaces des mailles modele
    !      ( L'indice 1 correspond a l'interface mailLE 1 / maille 2)
    !      ( Les abscisses sont exprimées en degres)

    use nr_util, only: assert_eq

    IMPLICIT NONE

    REAL, intent(in):: dlonid(:)
    real, intent(in):: fdat(:)
    real, intent(in):: rlonimod(:)

    real inter_barx(size(rlonimod))

    !    ...  Variables locales ... 

    INTEGER idatmax, imodmax
    REAL xxid(size(dlonid)+1), xxd(size(dlonid)+1), fdd(size(dlonid)+1)
    REAL  fxd(size(dlonid)+1), xchan(size(dlonid)+1), fdchan(size(dlonid)+1) 
    REAL  xxim(size(rlonimod))

    REAL x0, xim0, dx, dxm
    REAL chmin, chmax, pi

    INTEGER imod, idat, i, ichang, id0, id1, nid, idatmax1

    !-----------------------------------------------------

    idatmax = assert_eq(size(dlonid), size(fdat), "inter_barx idatmax")
    imodmax = size(rlonimod)

    pi = 2. * ASIN(1.)

    !   REDEFINITION DE L'ORIGINE DES ABSCISSES
    !    A L'INTERFACE OUEST DE LA PREMIERE MAILLE DU MODELE  
    DO imod = 1, imodmax
       xxim(imod) = rlonimod(imod)
    ENDDO

    CALL minmax( imodmax, xxim, chmin, chmax)
    IF( chmax.LT.6.50 )   THEN
       DO imod = 1, imodmax
          xxim(imod) = xxim(imod) * 180./pi
       ENDDO
    ENDIF

    xim0 = xxim(imodmax) - 360.

    DO imod = 1, imodmax
       xxim(imod) = xxim(imod) - xim0
    ENDDO

    idatmax1 = idatmax +1

    DO idat = 1, idatmax
       xxd(idat) = dlonid(idat)
    ENDDO

    CALL minmax( idatmax, xxd, chmin, chmax)
    IF( chmax.LT.6.50 )  THEN
       DO idat = 1, idatmax
          xxd(idat) = xxd(idat) * 180./pi
       ENDDO
    ENDIF

    DO idat = 1, idatmax
       xxd(idat) = AMOD( xxd(idat) - xim0, 360. )
       fdd(idat) = fdat (idat)
    ENDDO

    i = 2
    DO while (xxd(i) >= xxd(i-1) .and. i < idatmax)
       i = i + 1
    ENDDO
    IF (xxd(i) < xxd(i-1)) THEN
       ichang = i
       !  ***  reorganisation  des longitudes entre 0. et 360. degres ****
       nid = idatmax - ichang +1
       DO i = 1, nid
          xchan (i) = xxd(i+ichang -1 )
          fdchan(i) = fdd(i+ichang -1 )
       ENDDO
       DO i=1, ichang -1
          xchan (i+ nid) = xxd(i)
          fdchan(i+nid) = fdd(i) 
       ENDDO
       DO i =1, idatmax
          xxd(i) = xchan(i)
          fdd(i) = fdchan(i)
       ENDDO
    end IF

    !    translation des champs de donnees par rapport
    !    a la nouvelle origine, avec redondance de la
    !       maille a cheval sur les bords

    id0 = 0
    id1 = 0

    DO idat = 1, idatmax
       IF ( xxd( idatmax1- idat ).LT.360.) exit
       id1 = id1 + 1
    ENDDO

    DO idat = 1, idatmax
       IF (xxd(idat).GT.0.) exit
       id0 = id0 + 1
    END DO

    IF( id1 /= 0 ) then
       DO idat = 1, id1
          xxid(idat) = xxd(idatmax - id1 + idat) - 360.
          fxd (idat) = fdd(idatmax - id1 + idat)     
       END DO
       DO idat = 1, idatmax - id1
          xxid(idat + id1) = xxd(idat)
          fxd (idat + id1) = fdd(idat)
       END DO
    end IF

    IF(id0 /= 0) then
       DO idat = 1, idatmax - id0
          xxid(idat) = xxd(idat + id0)
          fxd (idat) = fdd(idat + id0)
       END DO

       DO idat = 1, id0
          xxid (idatmax - id0 + idat) =  xxd(idat) + 360.
          fxd  (idatmax - id0 + idat) =  fdd(idat)   
       END DO
    else 
       DO idat = 1, idatmax
          xxid(idat)  = xxd(idat)
          fxd (idat)  = fdd(idat)
       ENDDO
    end IF
    xxid(idatmax1) = xxid(1) + 360.
    fxd (idatmax1) = fxd(1)

    !   initialisation du champ du modele

    inter_barx(:) = 0.

    ! iteration

    x0   = xim0
    dxm  = 0.
    imod = 1
    idat = 1

    do while (imod <= imodmax)
       do while (xxim(imod).GT.xxid(idat))
          dx   = xxid(idat) - x0
          dxm  = dxm + dx
          inter_barx(imod) = inter_barx(imod) + dx * fxd(idat)
          x0   = xxid(idat)
          idat = idat + 1
       end do
       IF (xxim(imod).LT.xxid(idat)) THEN
          dx   = xxim(imod) - x0
          dxm  = dxm + dx
          inter_barx(imod) = (inter_barx(imod) + dx * fxd(idat)) / dxm
          x0   = xxim(imod)
          dxm  = 0.
          imod = imod + 1
       ELSE
          dx   = xxim(imod) - x0
          dxm  = dxm + dx
          inter_barx(imod) = (inter_barx(imod) + dx * fxd(idat)) / dxm
          x0   = xxim(imod)
          dxm  = 0.
          imod = imod + 1
          idat = idat + 1
       END IF
    end do

  END function inter_barx

  !******************************

  function inter_bary(yjdat, fdat, yjmod)

    ! From dyn3d/inter_bary.F, version 1.1.1.1 2004/05/19 12:53:06
    ! Authors: R. Sadourny, P. Le Van

    ! Interpolation barycentrique basée sur les aires.
    ! Version unidimensionnelle, en latitude.
    ! L'indice 1 correspond à l'interface maille 1 -- maille 2.

    use nr_util, only: assert

    IMPLICIT NONE

    REAL, intent(in):: yjdat(:)
    ! (angles, ordonnées des interfaces des mailles des données, in
    ! degrees, in increasing order)

    REAL, intent(in):: fdat(:) ! champ de données

    REAL, intent(in):: yjmod(:)
    ! (ordonnées des interfaces des mailles du modèle)
    ! (in degrees, in strictly increasing order)

    REAL inter_bary(size(yjmod)) ! champ du modèle

    ! Variables local to the procedure:

    REAL y0, dy, dym 
    INTEGER jdat ! indice du champ de données
    integer jmod ! indice du champ du modèle

    !------------------------------------

    call assert(size(yjdat) == size(fdat), "inter_bary")

    ! Initialisation des variables
    inter_bary(:) = 0.
    y0    = -90.
    dym   = 0.
    jmod  = 1
    jdat  = 1

    do while (jmod <= size(yjmod))
       do while (yjmod(jmod) > yjdat(jdat))
          dy         = yjdat(jdat) - y0
          dym        = dym + dy
          inter_bary(jmod) = inter_bary(jmod) + dy * fdat(jdat)
          y0         = yjdat(jdat)
          jdat       = jdat + 1
       end do
       IF (yjmod(jmod) < yjdat(jdat)) THEN
          dy         = yjmod(jmod) - y0
          dym        = dym + dy
          inter_bary(jmod) = (inter_bary(jmod) + dy * fdat(jdat)) / dym
          y0         = yjmod(jmod)
          dym        = 0.
          jmod       = jmod + 1
       ELSE
          ! {yjmod(jmod) == yjdat(jdat)}
          dy         = yjmod(jmod) - y0
          dym        = dym + dy
          inter_bary(jmod) = (inter_bary(jmod) + dy * fdat(jdat)) / dym
          y0         = yjmod(jmod)
          dym        = 0.
          jmod       = jmod + 1
          jdat       = jdat + 1
       END IF
    end do
    ! Le test de fin suppose que l'interface 0 est commune aux deux
    ! grilles "yjdat" et "yjmod".

  END function inter_bary

  !******************************

  SUBROUTINE ord_coord(xi, xo, decrois)

    ! From dyn3d/ord_coord.F, version 1.1.1.1 2004/05/19 12:53:06
    ! Author : P. Le Van

    ! This procedure receives an array of latitudes.
    ! It converts them to degrees if they are in radians.
    ! If the input latitudes are in decreasing order, the procedure
    ! reverses their order.
    ! Finally, the procedure adds 90° as the last value of the array.

    use nr_util, only: assert_eq, pi

    IMPLICIT NONE

    REAL, intent(in):: xi(:)
    ! (latitude, in degrees or radians, in increasing or decreasing order)
    ! ("xi" should contain latitudes from pole to pole.
    ! "xi" should contain the latitudes of the boundaries of grid
    ! cells, not the centers of grid cells.
    ! So the extreme values should not be 90° and -90°.)

    REAL, intent(out):: xo(:) ! angles in degrees
    LOGICAL, intent(out):: decrois

    ! Variables  local to the procedure:
    INTEGER nmax, i

    !--------------------

    nmax = assert_eq(size(xi), size(xo) - 1, "ord_coord")

    ! Check monotonicity:
    decrois = xi(2) < xi(1)
    DO i = 3, nmax
       IF (decrois .neqv. xi(i) < xi(i-1)) then
          print *, '"ord_coord":  latitudes are not monotonic'
          stop 1
       end IF
    ENDDO

    IF (abs(xi(1)) < pi) then
       ! "xi" contains latitudes in radians
       xo(:nmax) = xi(:) * 180. / pi
    else
       ! "xi" contains latitudes in degrees
       xo(:nmax) = xi(:)
    end IF

    IF (ABS(abs(xo(1)) - 90) < 0.001 .or. ABS(abs(xo(nmax)) - 90) < 0.001) THEN
       print *, "ord_coord"
       PRINT *, '"xi" should contain the latitudes of the boundaries of ' &
            // 'grid cells, not the centers of grid cells.'
       STOP 1
    ENDIF

    IF (decrois) xo(:nmax) = xo(nmax:1:- 1)
    xo(nmax + 1) = 90.

  END SUBROUTINE ord_coord

  !***********************************

  function ord_coordm(xi)

    ! From dyn3d/ord_coordm.F, version 1.1.1.1 2004/05/19 12:53:06
    ! Author : P. Le Van

    ! This procedure converts to degrees, if necessary, and inverts the
    ! order.

    use nr_util, only: pi

    IMPLICIT NONE

    REAL, intent(in):: xi(:) ! angle, in rad or degrees
    REAL ord_coordm(size(xi)) ! angle, in degrees

    !-----------------------------

    IF (xi(1) < 6.5) THEN
       ! "xi" is in rad
       ord_coordm(:) = xi(size(xi):1:-1) * 180. / pi
    else
       ! "xi" is in degrees
       ord_coordm(:) = xi(size(xi):1:-1)
    ENDIF

  END function ord_coordm

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