trunk/dyn3d/grille_m.f

module grille_m_m

  ! From grid_atob.F, v 1.1.1.1, 2004/05/19 12:53:05

  IMPLICIT none

contains

  function grille_m(xdata, ydata, entree, x, y)

    ! Z. X. Li (1er avril 1994) (voir aussi A. Harzallah et L. Fairhead)

    ! M\'ethode na\"ive pour transformer un champ d'une grille fine \`a une
    ! grille grossi\`ere. Je consid\`ere que les nouveaux points occupent
    ! une zone adjacente qui comprend un ou plusieurs anciens points.

    ! Aucune pond\'eration n'est consid\'er\'ee (voir
    ! grille_p). Cf. grille_m.txt.

    use  dist_sphe_m, only: dist_sphe
    use nr_util, only: assert_eq

    ! Coordonn\'ees :
    REAL, intent(in):: xdata(:) ! (imdep)
    REAL, intent(in):: ydata(:) ! (jmdep)

    REAL, intent(in):: entree(:, :) ! (imdep, jmdep) champ \`a transformer

    ! Coordonn\'ees :
    REAL, intent(in):: x(:) ! (imar)
    REAL, intent(in):: y(:) ! (jmar)

    real grille_m(size(x), size(y)) ! (imar, jmar) champ transform\'e

    ! Local:
    INTEGER imdep, jmdep, imar, jmar
    INTEGER i, j, ii, jj
    REAL a(size(x)), b(size(x)) ! (imar)
    real c(size(y)), d(size(y)) ! (jmar)
    REAL number(size(x), size(y)) ! (imar, jmar) 
    REAL distans(size(xdata) * size(ydata)) ! (imdep * jmdep)
    INTEGER i_proche, j_proche, ij_proche
    REAL zzmin

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

    imdep = assert_eq(size(xdata), size(entree, 1), "grille_m")
    jmdep = assert_eq(size(ydata), size(entree, 2), "grille_m")
    imar = size(x)
    jmar = size(y)

    ! Calculer les limites des zones des nouveaux points

    a(1) = x(1) - (x(2)-x(1))/2.0
    b(1) = (x(1)+x(2))/2.0
    DO i = 2, imar-1
       a(i) = b(i-1)
       b(i) = (x(i)+x(i+1))/2.0
    ENDDO
    a(imar) = b(imar-1)
    b(imar) = x(imar) + (x(imar)-x(imar-1))/2.0

    c(1) = y(1) - (y(2)-y(1))/2.0
    d(1) = (y(1)+y(2))/2.0
    DO j = 2, jmar-1
       c(j) = d(j-1)
       d(j) = (y(j)+y(j+1))/2.0
    ENDDO
    c(jmar) = d(jmar-1)
    d(jmar) = y(jmar) + (y(jmar)-y(jmar-1))/2.0

    DO i = 1, imar
       DO j = 1, jmar
          number(i, j) = 0.0
          grille_m(i, j) = 0.0
       ENDDO
    ENDDO

    ! Determiner la zone sur laquelle chaque ancien point se trouve

    DO ii = 1, imar
       DO jj = 1, jmar
          DO i = 1, imdep
             IF((xdata(i)-a(ii) >= 1.e-5.AND.xdata(i)-b(ii) <= 1.e-5).OR. &
                  (xdata(i)-a(ii) <= 1.e-5.AND.xdata(i)-b(ii) >= 1.e-5)) &
                  THEN
                DO j = 1, jmdep
                   IF((ydata(j)-c(jj) >= 1.e-5.AND.ydata(j)-d(jj) <= 1.e-5) &
                        .OR. (ydata(j)-c(jj) <= 1.e-5 .AND. &
                        ydata(j)-d(jj) >= 1.e-5)) THEN
                      number(ii, jj) = number(ii, jj) + 1.0
                      grille_m(ii, jj) = grille_m(ii, jj) + entree(i, j)
                   ENDIF
                ENDDO
             ENDIF
          ENDDO
       ENDDO
    ENDDO

    DO i = 1, imar
       DO j = 1, jmar
          IF (number(i, j) > 0.001) THEN
             grille_m(i, j) = grille_m(i, j) / number(i, j)
          ELSE
             ! Si aucun ancien point ne tombe sur une zone, c'est un probl\`eme
             CALL dist_sphe(x(i), y(j), xdata, ydata, imdep, jmdep, distans)
             ij_proche = 1
             zzmin = distans(ij_proche)
             DO ii = 2, imdep*jmdep
                IF (distans(ii) < zzmin) THEN
                   zzmin = distans(ii)
                   ij_proche = ii
                ENDIF
             ENDDO
             j_proche = (ij_proche-1)/imdep + 1
             i_proche = ij_proche - (j_proche-1)*imdep
             grille_m(i, j) = entree(i_proche, j_proche)
          ENDIF
       ENDDO
    ENDDO

    if (any(number <= 0.001)) print *, "problem in grille_m"

  END function grille_m

end module grille_m_m
1	guez	212	module grille_m_m
2	guez	3
3	guez	212	! From grid_atob.F, v 1.1.1.1, 2004/05/19 12:53:05
4	guez	3
5			IMPLICIT none
6
7			contains
8
9	guez	15	function grille_m(xdata, ydata, entree, x, y)
10	guez	3
11	guez	212	! Z. X. Li (1er avril 1994) (voir aussi A. Harzallah et L. Fairhead)
12	guez	3
13	guez	212	! M\'ethode na\"ive pour transformer un champ d'une grille fine \`a une
14			! grille grossi\`ere. Je consid\`ere que les nouveaux points occupent
15			! une zone adjacente qui comprend un ou plusieurs anciens points.
16	guez	3
17	guez	212	! Aucune pond\'eration n'est consid\'er\'ee (voir
18			! grille_p). Cf. grille_m.txt.
19	guez	3
20	guez	212	use dist_sphe_m, only: dist_sphe
21	guez	36	use nr_util, only: assert_eq
22	guez	3
23	guez	212	! Coordonn\'ees :
24			REAL, intent(in):: xdata(:) ! (imdep)
25			REAL, intent(in):: ydata(:) ! (jmdep)
26	guez	3
27	guez	212	REAL, intent(in):: entree(:, :) ! (imdep, jmdep) champ \`a transformer
28	guez	3
29	guez	212	! Coordonn\'ees :
30			REAL, intent(in):: x(:) ! (imar)
31			REAL, intent(in):: y(:) ! (jmar)
32
33			real grille_m(size(x), size(y)) ! (imar, jmar) champ transform\'e
34
35			! Local:
36	guez	3	INTEGER imdep, jmdep, imar, jmar
37			INTEGER i, j, ii, jj
38	guez	212	REAL a(size(x)), b(size(x)) ! (imar)
39			real c(size(y)), d(size(y)) ! (jmar)
40			REAL number(size(x), size(y)) ! (imar, jmar)
41			REAL distans(size(xdata) * size(ydata)) ! (imdep * jmdep)
42	guez	3	INTEGER i_proche, j_proche, ij_proche
43			REAL zzmin
44
45			!-------------------------
46
47			imdep = assert_eq(size(xdata), size(entree, 1), "grille_m")
48			jmdep = assert_eq(size(ydata), size(entree, 2), "grille_m")
49			imar = size(x)
50			jmar = size(y)
51
52			! Calculer les limites des zones des nouveaux points
53
54			a(1) = x(1) - (x(2)-x(1))/2.0
55			b(1) = (x(1)+x(2))/2.0
56			DO i = 2, imar-1
57			a(i) = b(i-1)
58			b(i) = (x(i)+x(i+1))/2.0
59			ENDDO
60			a(imar) = b(imar-1)
61			b(imar) = x(imar) + (x(imar)-x(imar-1))/2.0
62
63			c(1) = y(1) - (y(2)-y(1))/2.0
64			d(1) = (y(1)+y(2))/2.0
65			DO j = 2, jmar-1
66			c(j) = d(j-1)
67			d(j) = (y(j)+y(j+1))/2.0
68			ENDDO
69			c(jmar) = d(jmar-1)
70			d(jmar) = y(jmar) + (y(jmar)-y(jmar-1))/2.0
71
72			DO i = 1, imar
73			DO j = 1, jmar
74	guez	212	number(i, j) = 0.0
75			grille_m(i, j) = 0.0
76	guez	3	ENDDO
77			ENDDO
78
79			! Determiner la zone sur laquelle chaque ancien point se trouve
80
81			DO ii = 1, imar
82			DO jj = 1, jmar
83			DO i = 1, imdep
84	guez	212	IF((xdata(i)-a(ii) >= 1.e-5.AND.xdata(i)-b(ii) <= 1.e-5).OR. &
85			(xdata(i)-a(ii) <= 1.e-5.AND.xdata(i)-b(ii) >= 1.e-5)) &
86	guez	3	THEN
87			DO j = 1, jmdep
88	guez	212	IF((ydata(j)-c(jj) >= 1.e-5.AND.ydata(j)-d(jj) <= 1.e-5) &
89	guez	3	.OR. (ydata(j)-c(jj) <= 1.e-5 .AND. &
90			ydata(j)-d(jj) >= 1.e-5)) THEN
91	guez	212	number(ii, jj) = number(ii, jj) + 1.0
92			grille_m(ii, jj) = grille_m(ii, jj) + entree(i, j)
93	guez	3	ENDIF
94			ENDDO
95			ENDIF
96			ENDDO
97			ENDDO
98			ENDDO
99
100			DO i = 1, imar
101			DO j = 1, jmar
102	guez	212	IF (number(i, j) > 0.001) THEN
103			grille_m(i, j) = grille_m(i, j) / number(i, j)
104	guez	3	ELSE
105	guez	212	! Si aucun ancien point ne tombe sur une zone, c'est un probl\`eme
106			CALL dist_sphe(x(i), y(j), xdata, ydata, imdep, jmdep, distans)
107	guez	3	ij_proche = 1
108			zzmin = distans(ij_proche)
109			DO ii = 2, imdep*jmdep
110	guez	212	IF (distans(ii) < zzmin) THEN
111	guez	3	zzmin = distans(ii)
112			ij_proche = ii
113			ENDIF
114			ENDDO
115			j_proche = (ij_proche-1)/imdep + 1
116			i_proche = ij_proche - (j_proche-1)*imdep
117	guez	212	grille_m(i, j) = entree(i_proche, j_proche)
118	guez	3	ENDIF
119			ENDDO
120			ENDDO
121
122	guez	212	if (any(number <= 0.001)) print *, "problem in grille_m"
123
124	guez	3	END function grille_m
125
126	guez	212	end module grille_m_m