source: branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/TOOLS/GRIDGEN/src/projection.f90 @ 5445

Last change on this file since 5445 was 5445, checked in by davestorkey, 5 years ago

Clear SVN keywords from 2015/dev_r5021_UKMO1_CICE_coupling branch.

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.