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

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

    print *, "Call sequence information: grille_m"

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