[2136] | 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 |
---|