1 | MODULE projection |
---|
2 | !!----------------------------------------------------------- |
---|
3 | !! |
---|
4 | !! to make a polar stereographic projection |
---|
5 | !! in the area of the north pole |
---|
6 | !! |
---|
7 | !! Created by Brice Lemaire on 05/2010. |
---|
8 | !! |
---|
9 | !!----------------------------------------------------------- |
---|
10 | USE readwrite |
---|
11 | ! |
---|
12 | IMPLICIT NONE |
---|
13 | PUBLIC |
---|
14 | ! |
---|
15 | REAL*8,DIMENSION(4) :: dmixlam, dmixphi |
---|
16 | REAL*8 :: dray = 6357 !rayon polaire de la Terre (km) |
---|
17 | REAL*8 :: PI = ACOS(-1.0) |
---|
18 | REAL*8 :: dlong0 = 0. !longitude du centre de projection (origine) |
---|
19 | REAL*8 :: dlat0 = 90. !latitude '' |
---|
20 | INTEGER :: idisc = 0 !to check the +/- 180 discontinuity |
---|
21 | INTEGER :: m |
---|
22 | ! |
---|
23 | CONTAINS |
---|
24 | !******************************************************** |
---|
25 | ! SUBROUTINE stereo_projection * |
---|
26 | ! * |
---|
27 | ! * |
---|
28 | ! CALLED by cfg_tools.F90 * |
---|
29 | !******************************************************** |
---|
30 | SUBROUTINE stereo_projection(kji,kjj,kjk,knorth_pole,kway) |
---|
31 | ! |
---|
32 | INTEGER,INTENT(IN) :: kji,kjj,kjk |
---|
33 | LOGICAL,INTENT(IN) :: knorth_pole |
---|
34 | INTEGER,INTENT(IN) :: kway |
---|
35 | INTEGER :: klon, klat |
---|
36 | REAL*8,DIMENSION(4) :: dlx, dly |
---|
37 | REAL*8,DIMENSION(4) :: dllam, dlphi, dlk |
---|
38 | REAL*8 :: dl_lat0, dl_long0 |
---|
39 | ! |
---|
40 | dl_lat0 = dlat0 * PI/180. |
---|
41 | dl_long0 = dlong0 * PI/180. |
---|
42 | ! |
---|
43 | !Either we interpolate along longitude or along latitude |
---|
44 | IF(kway.EQ.1)THEN |
---|
45 | klon = 1 |
---|
46 | klat = 0 |
---|
47 | ELSEIF(kway.EQ.2) THEN |
---|
48 | klon = 0 |
---|
49 | klat = 1 |
---|
50 | ENDIF |
---|
51 | ! |
---|
52 | ! Check the discontinuity +/- 180 |
---|
53 | IF(kway.EQ.1)THEN |
---|
54 | IF(smixgrd%glam(kji,kjj).LT.180..AND.smixgrd%glam(kji+3*nn_rhox,kjj).GT.180.) THEN |
---|
55 | IF(smixgrd%glam(kji+2*nn_rhox,kjj).LT.180.) THEN |
---|
56 | idisc = 1 |
---|
57 | ELSEIF(smixgrd%glam(kji+1*nn_rhox,kjj).GT.180.) THEN |
---|
58 | idisc = 2 |
---|
59 | ELSEIF(smixgrd%glam(kji+2*nn_rhox,kjj).GT.180.AND.smixgrd%glam(kji+1*nn_rhox,kjj).LT.180.)THEN |
---|
60 | idisc = 3 |
---|
61 | ENDIF |
---|
62 | ELSE |
---|
63 | idisc = 0 |
---|
64 | ENDIF |
---|
65 | ENDIF |
---|
66 | ! |
---|
67 | dllam(1) = smixgrd%glam(kji,kjj) * PI/180. |
---|
68 | dlphi(1) = smixgrd%gphi(kji,kjj) * PI/180. |
---|
69 | dlk(1) = (2*dray) / (1 + (SIN(dl_lat0)*SIN(dlphi(1))) + (COS(dl_lat0)*COS(dlphi(1))*COS(dllam(1) - dl_long0))) |
---|
70 | ! |
---|
71 | dllam(2) = smixgrd%glam(kji+1*nn_rhox*klon,kjj+1*nn_rhoy*klat) * PI/180. |
---|
72 | dlphi(2) = smixgrd%gphi(kji+1*nn_rhox*klon,kjj+1*nn_rhoy*klat) * PI/180. |
---|
73 | dlk(2) = (2*dray) / (1 + (SIN(dl_lat0)*SIN(dlphi(2))) + (COS(dl_lat0)*COS(dlphi(2))*COS(dllam(2) - dl_long0))) |
---|
74 | ! |
---|
75 | dllam(3) = smixgrd%glam(kji+2*nn_rhox*klon,kjj+2*nn_rhoy*klat) * PI/180. |
---|
76 | dlphi(3) = smixgrd%gphi(kji+2*nn_rhox*klon,kjj+2*nn_rhoy*klat) * PI/180. |
---|
77 | dlk(3) = (2*dray) / (1 + (SIN(dl_lat0)*SIN(dlphi(3))) + (COS(dl_lat0)*COS(dlphi(3))*COS(dllam(3) - dl_long0))) |
---|
78 | ! |
---|
79 | dllam(4) = smixgrd%glam(kji+3*nn_rhox*klon,kjj+3*nn_rhoy*klat) * PI/180. |
---|
80 | dlphi(4) = smixgrd%gphi(kji+3*nn_rhox*klon,kjj+3*nn_rhoy*klat) * PI/180. |
---|
81 | dlk(4) = (2*dray) / (1 + (SIN(dl_lat0)*SIN(dlphi(4))) + (COS(dl_lat0)*COS(dlphi(4))*COS(dllam(4) - dl_long0))) |
---|
82 | ! |
---|
83 | dlx(1) = dlk(1) * COS(dlphi(1)) * SIN(dllam(1) - dl_long0) |
---|
84 | dly(1) = dlk(1) * ((COS(dl_lat0) * SIN(dlphi(1))) - (SIN(dl_lat0) * COS(dlphi(1)) * COS(dllam(1) - dl_long0))) |
---|
85 | ! |
---|
86 | dlx(2) = dlk(2) * COS(dlphi(2)) * SIN(dllam(2) - dl_long0) |
---|
87 | dly(2) = dlk(2) * ((COS(dl_lat0) * SIN(dlphi(2))) - (SIN(dl_lat0) * COS(dlphi(2)) * COS(dllam(2) - dl_long0))) |
---|
88 | ! |
---|
89 | dlx(3) = dlk(3) * COS(dlphi(3)) * SIN(dllam(3) - dl_long0) |
---|
90 | dly(3) = dlk(3) * ((COS(dl_lat0) * SIN(dlphi(3))) - (SIN(dl_lat0) * COS(dlphi(3)) * COS(dllam(3) - dl_long0))) |
---|
91 | ! |
---|
92 | dlx(4) = dlk(4) * COS(dlphi(4)) * SIN(dllam(4) - dl_long0) |
---|
93 | dly(4) = dlk(4) * ((COS(dl_lat0) * SIN(dlphi(4))) - (SIN(dl_lat0) * COS(dlphi(4)) * COS(dllam(4) - dl_long0))) |
---|
94 | ! |
---|
95 | dmixlam(1) = smixgrd%glam(kji,kjj) |
---|
96 | dmixlam(2) = smixgrd%glam(kji+1*nn_rhox*klon,kjj+1*nn_rhoy*klat) |
---|
97 | dmixlam(3) = smixgrd%glam(kji+2*nn_rhox*klon,kjj+2*nn_rhoy*klat) |
---|
98 | dmixlam(4) = smixgrd%glam(kji+3*nn_rhox*klon,kjj+3*nn_rhoy*klat) |
---|
99 | ! |
---|
100 | dmixphi(1) = smixgrd%gphi(kji,kjj) |
---|
101 | dmixphi(2) = smixgrd%gphi(kji+1*nn_rhox*klon,kjj+1*nn_rhoy*klat) |
---|
102 | dmixphi(3) = smixgrd%gphi(kji+2*nn_rhox*klon,kjj+2*nn_rhoy*klat) |
---|
103 | dmixphi(4) = smixgrd%gphi(kji+3*nn_rhox*klon,kjj+3*nn_rhoy*klat) |
---|
104 | ! |
---|
105 | smixgrd%glam(kji,kjj) = dlx(1) |
---|
106 | smixgrd%glam(kji+1*nn_rhox*klon,kjj+1*nn_rhoy*klat) = dlx(2) |
---|
107 | smixgrd%glam(kji+2*nn_rhox*klon,kjj+2*nn_rhoy*klat) = dlx(3) |
---|
108 | smixgrd%glam(kji+3*nn_rhox*klon,kjj+3*nn_rhoy*klat) = dlx(4) |
---|
109 | ! |
---|
110 | smixgrd%gphi(kji,kjj) = dly(1) |
---|
111 | smixgrd%gphi(kji+1*nn_rhox*klon,kjj+1*nn_rhoy*klat) = dly(2) |
---|
112 | smixgrd%gphi(kji+2*nn_rhox*klon,kjj+2*nn_rhoy*klat) = dly(3) |
---|
113 | smixgrd%gphi(kji+3*nn_rhox*klon,kjj+3*nn_rhoy*klat) = dly(4) |
---|
114 | ! |
---|
115 | END SUBROUTINE stereo_projection |
---|
116 | ! |
---|
117 | ! |
---|
118 | ! |
---|
119 | !******************************************************** |
---|
120 | ! SUBROUTINE stereo_projection_inv * |
---|
121 | ! * |
---|
122 | ! * |
---|
123 | ! CALLED from here * |
---|
124 | !******************************************************** |
---|
125 | SUBROUTINE stereo_projection_inv(kji,kjj,kjk,knorth_pole,kway) |
---|
126 | ! |
---|
127 | INTEGER,INTENT(IN) :: kji,kjj,kjk |
---|
128 | LOGICAL,INTENT(IN) :: knorth_pole |
---|
129 | INTEGER,INTENT(IN) :: kway |
---|
130 | INTEGER :: klon, klat |
---|
131 | REAL*8,DIMENSION(5) :: dlx, dly |
---|
132 | REAL*8,DIMENSION(5) :: dllam, dlphi |
---|
133 | REAL*8,DIMENSION(5) :: dlro, dlc |
---|
134 | REAL*8 :: dl_long0, dl_lat0 |
---|
135 | ! |
---|
136 | dl_lat0 = dlat0 * PI/180. |
---|
137 | dl_long0 = dlong0 * PI/180. |
---|
138 | ! |
---|
139 | IF(kway.EQ.1)THEN |
---|
140 | klon = 1 |
---|
141 | klat = 0 |
---|
142 | ELSEIF(kway.EQ.2) THEN |
---|
143 | klon = 0 |
---|
144 | klat = 1 |
---|
145 | ENDIF |
---|
146 | ! |
---|
147 | ! |
---|
148 | dlx(1) = smixgrd%glam(kji,kjj) |
---|
149 | dlx(2) = smixgrd%glam(kji+1*nn_rhox*klon,kjj+1*nn_rhoy*klat) |
---|
150 | dlx(3) = smixgrd%glam(kji+2*nn_rhox*klon,kjj+2*nn_rhoy*klat) |
---|
151 | dlx(4) = smixgrd%glam(kji+3*nn_rhox*klon,kjj+3*nn_rhoy*klat) |
---|
152 | dlx(5) = smixgrd%glam(kji+(nn_rhox+kjk)*klon,kjj+(nn_rhoy+kjk)*klat) |
---|
153 | ! |
---|
154 | dly(1) = smixgrd%gphi(kji,kjj) |
---|
155 | dly(2) = smixgrd%gphi(kji+1*nn_rhox*klon,kjj+1*nn_rhoy*klat) |
---|
156 | dly(3) = smixgrd%gphi(kji+2*nn_rhox*klon,kjj+2*nn_rhoy*klat) |
---|
157 | dly(4) = smixgrd%gphi(kji+3*nn_rhox*klon,kjj+3*nn_rhoy*klat) |
---|
158 | dly(5) = smixgrd%gphi(kji+(nn_rhox+kjk)*klon,kjj+(nn_rhoy+kjk)*klat) |
---|
159 | ! |
---|
160 | dlro(1) = SQRT(dlx(1)*dlx(1) + dly(1)*dly(1)) |
---|
161 | dlro(2) = SQRT(dlx(2)*dlx(2) + dly(2)*dly(2)) |
---|
162 | dlro(3) = SQRT(dlx(3)*dlx(3) + dly(3)*dly(3)) |
---|
163 | dlro(4) = SQRT(dlx(4)*dlx(4) + dly(4)*dly(4)) |
---|
164 | dlro(5) = SQRT(dlx(5)*dlx(5) + dly(5)*dly(5)) |
---|
165 | ! |
---|
166 | dlc(1) = 2*ATAN(dlro(1)/(2*dray)) |
---|
167 | dlc(2) = 2*ATAN(dlro(2)/(2*dray)) |
---|
168 | dlc(3) = 2*ATAN(dlro(3)/(2*dray)) |
---|
169 | dlc(4) = 2*ATAN(dlro(4)/(2*dray)) |
---|
170 | dlc(5) = 2*ATAN(dlro(5)/(2*dray)) |
---|
171 | ! |
---|
172 | dlphi(1) = ASIN(COS(dlc(1))*SIN(dl_lat0) + (dly(1)*SIN(dlc(1))*COS(dl_lat0) / dlro(1))) |
---|
173 | dlphi(2) = ASIN(COS(dlc(2))*SIN(dl_lat0) + (dly(2)*SIN(dlc(2))*COS(dl_lat0) / dlro(2))) |
---|
174 | dlphi(3) = ASIN(COS(dlc(3))*SIN(dl_lat0) + (dly(3)*SIN(dlc(3))*COS(dl_lat0) / dlro(3))) |
---|
175 | dlphi(4) = ASIN(COS(dlc(4))*SIN(dl_lat0) + (dly(4)*SIN(dlc(4))*COS(dl_lat0) / dlro(4))) |
---|
176 | dlphi(5) = ASIN(COS(dlc(5))*SIN(dl_lat0) + (dly(5)*SIN(dlc(5))*COS(dl_lat0) / dlro(5))) |
---|
177 | ! |
---|
178 | dlphi(:) = dlphi(:) * 180./PI |
---|
179 | ! |
---|
180 | dllam(1) = dl_long0 + ATAN(dlx(1)*SIN(dlc(1) / ((dlro(1)*COS(dl_lat0)*COS(dlc(1)))-(dly(1)*SIN(dl_lat0)*SIN(dlc(1)))))) |
---|
181 | dllam(2) = dl_long0 + ATAN(dlx(2)*SIN(dlc(2) / ((dlro(2)*COS(dl_lat0)*COS(dlc(2)))-(dly(2)*SIN(dl_lat0)*SIN(dlc(2)))))) |
---|
182 | dllam(3) = dl_long0 + ATAN(dlx(3)*SIN(dlc(3) / ((dlro(3)*COS(dl_lat0)*COS(dlc(3)))-(dly(3)*SIN(dl_lat0)*SIN(dlc(3)))))) |
---|
183 | dllam(4) = dl_long0 + ATAN(dlx(4)*SIN(dlc(4) / ((dlro(4)*COS(dl_lat0)*COS(dlc(4)))-(dly(4)*SIN(dl_lat0)*SIN(dlc(4)))))) |
---|
184 | dllam(5) = dl_long0 + ATAN(dlx(5)*SIN(dlc(5) / ((dlro(5)*COS(dl_lat0)*COS(dlc(5)))-(dly(5)*SIN(dl_lat0)*SIN(dlc(5)))))) |
---|
185 | ! |
---|
186 | dllam(1) = dllam(1) * 180./PI |
---|
187 | dllam(2) = dllam(2) * 180./PI |
---|
188 | dllam(3) = dllam(3) * 180./PI |
---|
189 | dllam(4) = dllam(4) * 180./PI |
---|
190 | ! |
---|
191 | dllam(5) = dllam(5) * 180./PI |
---|
192 | ! |
---|
193 | IF(kway.EQ.1)THEN |
---|
194 | IF(idisc.EQ.1) THEN |
---|
195 | dllam(5) = dllam(5) + 180. |
---|
196 | ELSEIF(idisc.EQ.2) THEN |
---|
197 | dllam(5) = dllam(5) - 180. |
---|
198 | ELSEIF(idisc.EQ.3) THEN |
---|
199 | dllam(5) = dllam(5) - 180. |
---|
200 | ELSE |
---|
201 | dllam(5) = dllam(5) |
---|
202 | ENDIF |
---|
203 | ! |
---|
204 | ELSEIF(kway.EQ.2)THEN |
---|
205 | IF(smixgrd%glam(kji,kjj+1*nn_rhoy).GT.0.AND.smixgrd%glam(kji,kjj+2*nn_rhoy).GT.0)THEN |
---|
206 | IF(dllam(5).LT.0.)THEN |
---|
207 | dllam(5) = dllam(5) + 180 |
---|
208 | ENDIF |
---|
209 | ELSEIF(smixgrd%glam(kji,kjj+1*nn_rhoy).LT.0.AND.smixgrd%glam(kji,kjj+2*nn_rhoy).LT.0)THEN |
---|
210 | IF(dllam(5).GT.0.)THEN |
---|
211 | dllam(5) = dllam(5) - 180 |
---|
212 | ENDIF |
---|
213 | ELSEIF(smixgrd%glam(kji,kjj+1*nn_rhoy)*smixgrd%glam(kji,kjj+2*nn_rhoy).LT.0)THEN |
---|
214 | IF(dllam(5).LT.0.)THEN |
---|
215 | dllam(5) = dllam(5) + 180 |
---|
216 | ELSEIF(dllam(5).GT.0.)THEN |
---|
217 | dllam(5) = dllam(5) - 180 |
---|
218 | ENDIF |
---|
219 | ENDIF |
---|
220 | ENDIF |
---|
221 | ! |
---|
222 | smixgrd%glam(kji,kjj) = dmixlam(1) |
---|
223 | smixgrd%glam(kji+1*nn_rhox*klon,kjj+1*nn_rhoy*klat) = dmixlam(2) |
---|
224 | smixgrd%glam(kji+2*nn_rhox*klon,kjj+2*nn_rhoy*klat) = dmixlam(3) |
---|
225 | smixgrd%glam(kji+3*nn_rhox*klon,kjj+3*nn_rhoy*klat) = dmixlam(4) |
---|
226 | ! |
---|
227 | smixgrd%glam(kji+(nn_rhox+kjk)*klon,kjj+(nn_rhoy+kjk)*klat) = dllam(5) |
---|
228 | ! |
---|
229 | smixgrd%gphi(kji,kjj) = dmixphi(1) |
---|
230 | smixgrd%gphi(kji+1*nn_rhox*klon,kjj+1*nn_rhoy*klat) = dmixphi(2) |
---|
231 | smixgrd%gphi(kji+2*nn_rhox*klon,kjj+2*nn_rhoy*klat) = dmixphi(3) |
---|
232 | smixgrd%gphi(kji+3*nn_rhox*klon,kjj+3*nn_rhoy*klat) = dmixphi(4) |
---|
233 | ! |
---|
234 | smixgrd%gphi(kji+(nn_rhox+kjk)*klon,kjj+(nn_rhoy+kjk)*klat) = dlphi(5) |
---|
235 | ! |
---|
236 | END SUBROUTINE stereo_projection_inv |
---|
237 | ! |
---|
238 | ! |
---|
239 | ! |
---|
240 | END MODULE projection |
---|