source: CONFIG_DEVT/IPSLCM6.5_work_ENSEMBLES/modeles/NEMO/tools/GRIDGEN/src/projection.f90 @ 5501

Last change on this file since 5501 was 5501, checked in by aclsce, 4 years ago

First import of IPSLCM6.5_work_ENSEMBLES working configuration

File size: 9.3 KB
Line 
1MODULE 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  !
240END MODULE projection
Note: See TracBrowser for help on using the repository browser.