libf/dyn3d/grid_atob.f90

module grid_atob

  ! 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 (le 1 avril 1994) (voir aussi A. Harzallah et L. Fairhead)

    ! Méthode naïve pour transformer un champ d'une grille fine à une
    ! grille grossière. Je considère que les nouveaux points occupent
    ! une zone adjacente qui comprend un ou plusieurs anciens points

    ! Aucune pondération n'est considérée (voir grille_p)

    !           (c)
    !        ----d-----
    !        | . . . .|
    !        |        |
    !     (b)a . * . .b(a)
    !        |        |
    !        | . . . .|
    !        ----c-----
    !           (d)
    !=======================================================================
    ! INPUT:
    !        imdep, jmdep: dimensions X et Y pour depart
    !        xdata, ydata: coordonnees X et Y pour depart
    !        entree: champ d'entree a transformer
    ! OUTPUT:
    !        imar, jmar: dimensions X et Y d'arrivee
    !        x, y: coordonnees X et Y d'arrivee
    !        grille_m: champ de sortie deja transforme
    !=======================================================================

    use numer_rec, only: assert_eq

    REAL, intent(in):: xdata(:),ydata(:) 
    REAL, intent(in):: entree(:, :)
    REAL, intent(in):: x(:), y(:)

    real grille_m(size(x), size(y))

    ! Variables local to the procedure:
    INTEGER imdep, jmdep, imar, jmar
    INTEGER i, j, ii, jj
    REAL a(2200),b(2200),c(1100),d(1100)
    REAL number(2200,1100)
    REAL distans(2200*1100)
    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)

    IF (imar.GT.2200 .OR. jmar.GT.1100) THEN
       PRINT*, 'imar ou jmar trop grand', imar, jmar
       STOP 1
    ENDIF

    ! 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

    ! Si aucun ancien point ne tombe sur une zone, c'est un probleme

    DO i = 1, imar
       DO j = 1, jmar
          IF (number(i,j) .GT. 0.001) THEN
             grille_m(i,j) = grille_m(i,j) / number(i,j)
          ELSE
             PRINT*, 'probleme,i,j=', i,j
             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).LT.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
             PRINT*, "solution:", ij_proche, i_proche, j_proche
             grille_m(i,j) = entree(i_proche,j_proche)
          ENDIF
       ENDDO
    ENDDO

  END function grille_m

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

  SUBROUTINE dist_sphe(rf_lon,rf_lat,rlon,rlat,im,jm,distance)

    ! Auteur: Laurent Li (le 30 decembre 1996)

    ! Ce programme calcule la distance minimale (selon le grand cercle)
    ! entre deux points sur la terre

    INTEGER, intent(in):: im, jm ! dimensions
    REAL, intent(in):: rf_lon ! longitude du point de reference (degres)
    REAL, intent(in):: rf_lat ! latitude du point de reference (degres)
    REAL, intent(in):: rlon(im), rlat(jm) ! longitude et latitude des points

    REAL, intent(out):: distance(im,jm) ! distances en metre

    REAL rlon1, rlat1
    REAL rlon2, rlat2
    REAL dist
    REAL pa, pb, p, pi

    REAL radius
    PARAMETER (radius=6371229.)
    integer i, j

    pi = 4.0 * ATAN(1.0)

    DO j = 1, jm
       DO  i = 1, im
      
          rlon1=rf_lon
          rlat1=rf_lat
          rlon2=rlon(i)
          rlat2=rlat(j)
          pa = pi/2.0 - rlat1*pi/180.0 ! dist. entre pole n et point a
          pb = pi/2.0 - rlat2*pi/180.0 ! dist. entre pole n et point b
          p = (rlon1-rlon2)*pi/180.0 ! angle entre a et b (leurs meridiens)
      
          dist = ACOS( COS(pa)*COS(pb) + SIN(pa)*SIN(pb)*COS(p))
          dist = radius * dist
          distance(i,j) = dist
      
       end DO
    end DO

  END SUBROUTINE dist_sphe

end module grid_atob
1	module grid_atob
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	!=======================================================================
12	! Z. X. Li (le 1 avril 1994) (voir aussi A. Harzallah et L. Fairhead)
13
14	! Méthode naïve pour transformer un champ d'une grille fine à une
15	! grille grossière. Je considère que les nouveaux points occupent
16	! une zone adjacente qui comprend un ou plusieurs anciens points
17
18	! Aucune pondération n'est considérée (voir grille_p)
19
20	! (c)
21	! ----d-----
22	! \| . . . .\|
23	! \| \|
24	! (b)a . * . .b(a)
25	! \| \|
26	! \| . . . .\|
27	! ----c-----
28	! (d)
29	!=======================================================================
30	! INPUT:
31	! imdep, jmdep: dimensions X et Y pour depart
32	! xdata, ydata: coordonnees X et Y pour depart
33	! entree: champ d'entree a transformer
34	! OUTPUT:
35	! imar, jmar: dimensions X et Y d'arrivee
36	! x, y: coordonnees X et Y d'arrivee
37	! grille_m: champ de sortie deja transforme
38	!=======================================================================
39
40	use numer_rec, only: assert_eq
41
42	REAL, intent(in):: xdata(:),ydata(:)
43	REAL, intent(in):: entree(:, :)
44	REAL, intent(in):: x(:), y(:)
45
46	real grille_m(size(x), size(y))
47
48	! Variables local to the procedure:
49	INTEGER imdep, jmdep, imar, jmar
50	INTEGER i, j, ii, jj
51	REAL a(2200),b(2200),c(1100),d(1100)
52	REAL number(2200,1100)
53	REAL distans(2200*1100)
54	INTEGER i_proche, j_proche, ij_proche
55	REAL zzmin
56
57	!-------------------------
58
59	print *, "Call sequence information: grille_m"
60
61	imdep = assert_eq(size(xdata), size(entree, 1), "grille_m")
62	jmdep = assert_eq(size(ydata), size(entree, 2), "grille_m")
63	imar = size(x)
64	jmar = size(y)
65
66	IF (imar.GT.2200 .OR. jmar.GT.1100) THEN
67	PRINT*, 'imar ou jmar trop grand', imar, jmar
68	STOP 1
69	ENDIF
70
71	! Calculer les limites des zones des nouveaux points
72
73	a(1) = x(1) - (x(2)-x(1))/2.0
74	b(1) = (x(1)+x(2))/2.0
75	DO i = 2, imar-1
76	a(i) = b(i-1)
77	b(i) = (x(i)+x(i+1))/2.0
78	ENDDO
79	a(imar) = b(imar-1)
80	b(imar) = x(imar) + (x(imar)-x(imar-1))/2.0
81
82	c(1) = y(1) - (y(2)-y(1))/2.0
83	d(1) = (y(1)+y(2))/2.0
84	DO j = 2, jmar-1
85	c(j) = d(j-1)
86	d(j) = (y(j)+y(j+1))/2.0
87	ENDDO
88	c(jmar) = d(jmar-1)
89	d(jmar) = y(jmar) + (y(jmar)-y(jmar-1))/2.0
90
91	DO i = 1, imar
92	DO j = 1, jmar
93	number(i,j) = 0.0
94	grille_m(i,j) = 0.0
95	ENDDO
96	ENDDO
97
98	! Determiner la zone sur laquelle chaque ancien point se trouve
99
100	DO ii = 1, imar
101	DO jj = 1, jmar
102	DO i = 1, imdep
103	IF( ( xdata(i)-a(ii) >= 1.e-5.AND.xdata(i)-b(ii) <= 1.e-5 ).OR. &
104	(xdata(i)-a(ii) <= 1.e-5.AND.xdata(i)-b(ii) >= 1.e-5 ) ) &
105	THEN
106	DO j = 1, jmdep
107	IF((ydata(j)-c(jj) >= 1.e-5.AND.ydata(j)-d(jj) <= 1.e-5 ) &
108	.OR. (ydata(j)-c(jj) <= 1.e-5 .AND. &
109	ydata(j)-d(jj) >= 1.e-5)) THEN
110	number(ii,jj) = number(ii,jj) + 1.0
111	grille_m(ii,jj) = grille_m(ii,jj) + entree(i,j)
112	ENDIF
113	ENDDO
114	ENDIF
115	ENDDO
116	ENDDO
117	ENDDO
118
119	! Si aucun ancien point ne tombe sur une zone, c'est un probleme
120
121	DO i = 1, imar
122	DO j = 1, jmar
123	IF (number(i,j) .GT. 0.001) THEN
124	grille_m(i,j) = grille_m(i,j) / number(i,j)
125	ELSE
126	PRINT*, 'probleme,i,j=', i,j
127	CALL dist_sphe(x(i),y(j),xdata,ydata,imdep,jmdep,distans)
128	ij_proche = 1
129	zzmin = distans(ij_proche)
130	DO ii = 2, imdep*jmdep
131	IF (distans(ii).LT.zzmin) THEN
132	zzmin = distans(ii)
133	ij_proche = ii
134	ENDIF
135	ENDDO
136	j_proche = (ij_proche-1)/imdep + 1
137	i_proche = ij_proche - (j_proche-1)*imdep
138	PRINT*, "solution:", ij_proche, i_proche, j_proche
139	grille_m(i,j) = entree(i_proche,j_proche)
140	ENDIF
141	ENDDO
142	ENDDO
143
144	END function grille_m
145
146	!**************************************************
147
148	SUBROUTINE dist_sphe(rf_lon,rf_lat,rlon,rlat,im,jm,distance)
149
150	! Auteur: Laurent Li (le 30 decembre 1996)
151
152	! Ce programme calcule la distance minimale (selon le grand cercle)
153	! entre deux points sur la terre
154
155	INTEGER, intent(in):: im, jm ! dimensions
156	REAL, intent(in):: rf_lon ! longitude du point de reference (degres)
157	REAL, intent(in):: rf_lat ! latitude du point de reference (degres)
158	REAL, intent(in):: rlon(im), rlat(jm) ! longitude et latitude des points
159
160	REAL, intent(out):: distance(im,jm) ! distances en metre
161
162	REAL rlon1, rlat1
163	REAL rlon2, rlat2
164	REAL dist
165	REAL pa, pb, p, pi
166
167	REAL radius
168	PARAMETER (radius=6371229.)
169	integer i, j
170
171	pi = 4.0 * ATAN(1.0)
172
173	DO j = 1, jm
174	DO i = 1, im
175
176	rlon1=rf_lon
177	rlat1=rf_lat
178	rlon2=rlon(i)
179	rlat2=rlat(j)
180	pa = pi/2.0 - rlat1*pi/180.0 ! dist. entre pole n et point a
181	pb = pi/2.0 - rlat2*pi/180.0 ! dist. entre pole n et point b
182	p = (rlon1-rlon2)*pi/180.0 ! angle entre a et b (leurs meridiens)
183
184	dist = ACOS( COS(pa)COS(pb) + SIN(pa)SIN(pb)*COS(p))
185	dist = radius * dist
186	distance(i,j) = dist
187
188	end DO
189	end DO
190
191	END SUBROUTINE dist_sphe
192
193	end module grid_atob