/[lmdze]/trunk/dyn3d/grid_atob.f
ViewVC logotype

Annotation of /trunk/dyn3d/grid_atob.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 32 - (hide annotations)
Tue Apr 6 17:52:58 2010 UTC (14 years, 2 months ago) by guez
Original Path: trunk/libf/dyn3d/grid_atob.f90
File size: 5740 byte(s) 
Split "stringop.f90" into single-procedure files. Gathered files in directory
"IOIPSL/Stringop".

Split "flincom.f90" into "flincom.f90" and "flinget.f90". Removed
unused procedures from module "flincom". Removed unused argument
"filename" of procedure "flinopen_nozoom".

Removed unused files.

Split "grid_change.f90" into "grid_change.f90" and
"gr_phy_write_3d.f90".

Removed unused procedures from modules "calendar", "ioipslmpp",
"grid_atob", "gath_cpl" and "getincom". Removed unused procedures in
files "ppm3d.f" and "thermcell.f".

Split "mathelp.f90" into "mathelp.f90" and "mathop.f90".

Removed unused variable "dpres" of module "comvert".

Use argument "itau" instead of local variables "iadvtr" and "first" to
control algorithm in procedure "fluxstokenc".

Removed unused arguments of procedure "integrd".

Removed useless computations at the end of procedure "leapfrog".

Merged common block "matrfil" into module "parafilt".
1 guez 3 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 guez 15 function grille_m(xdata, ydata, entree, x, y)
10 guez 3
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 guez 13 use numer_rec, only: assert_eq
41 guez 3
42     REAL, intent(in):: xdata(:),ydata(:)
43     REAL, intent(in):: entree(:, :)
44     REAL, intent(in):: x(:), y(:)
45    
46 guez 15 real grille_m(size(x), size(y))
47 guez 3
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 guez 32 !**************************************************
147 guez 3
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
  ViewVC Help
Powered by ViewVC 1.1.21