36 public :: mapprojection_xy2lonlat
43 public :: mapprojection_lonlat2xy
50 public :: mapprojection_mapfactor
57 public :: mapprojection_rotcoef
64 interface mapprojection_xy2lonlat
66 module procedure mapprojection_xy2lonlat_2d_initialized
67 module procedure mapprojection_xy2lonlat_0d_param
70 interface mapprojection_lonlat2xy
72 module procedure mapprojection_lonlat2xy_2d_initialized
73 module procedure mapprojection_lonlat2xy_0d_param
77 interface mapprojection_mapfactor
79 module procedure mapprojection_mapfactor_param
80 end interface mapprojection_mapfactor
82 interface mapprojection_rotcoef
84 module procedure mapprojection_rotcoef_param
85 end interface mapprojection_rotcoef
94 character(len=H_SHORT) :: mapping_name
95 real(
dp) :: false_easting
96 real(
dp) :: false_northing
97 real(
dp) :: longitude_of_central_meridian
98 real(
dp) :: longitude_of_projection_origin
99 real(
dp) :: latitude_of_projection_origin
100 real(
dp) :: straight_vertical_longitude_from_pole
101 real(
dp) :: standard_parallel(2)
107 real(
dp) :: basepoint_x
108 real(
dp) :: basepoint_y
109 real(
dp) :: basepoint_lon
110 real(
dp) :: basepoint_lat
112 real(
dp) :: rot_fact_sin
113 real(
dp) :: rot_fact_cos
114 real(
dp) :: hemisphere
129 character(len=H_SHORT),
private :: mapprojection_type =
'NONE'
136 real(
dp),
private :: mapprojection_basepoint_x
137 real(
dp),
private :: mapprojection_basepoint_y
138 real(
dp),
private :: mapprojection_basepoint_lonrad
139 real(
dp),
private :: mapprojection_basepoint_latrad
140 real(
dp),
private :: mapprojection_rotation = 0.0_dp
142 real(
dp),
private :: mapprojection_lc_lat1 = 30.0_dp
143 real(
dp),
private :: mapprojection_lc_lat2 = 60.0_dp
144 real(
dp),
private :: mapprojection_ps_lat
145 real(
dp),
private :: mapprojection_m_lat = 0.0_dp
146 real(
dp),
private :: mapprojection_ec_lat = 0.0_dp
148 type(
mappingparam),
private :: mapprojection_mappingparam
150 real(
dp),
private :: pi
151 real(
dp),
private :: d2r
152 real(
dp),
private :: radius
155 subroutine xy2lonlat_s( x, y, param, lon, lat )
158 real(
rp),
intent(in) :: x, y
160 real(
rp),
intent(out) :: lon, lat
161 end subroutine xy2lonlat_s
165 real(RP),
intent(in) :: lon, lat
166 type(mappingparam),
intent(in) :: param
167 real(RP),
intent(out) :: x, y
172 real(RP),
intent(in) :: lat
173 type(mappingparam),
intent(in) :: param
174 real(RP),
intent(out) :: m1, m2
176 subroutine rotcoef_s( lon, lat, param, rotc_cos, rotc_sin )
179 real(RP),
intent(in) :: lon, lat
180 type(mappingparam),
intent(in) :: param
181 real(RP),
intent(out) :: rotc_cos, rotc_sin
185 procedure(
lonlat2xy_s),
pointer :: lonlat2xy => null()
186 procedure(
mapfactor_s),
pointer :: mapfactor => null()
187 procedure(
rotcoef_s),
pointer :: rotcoef => null()
205 real(
rp),
intent(in) :: domain_center_x
206 real(
rp),
intent(in) :: domain_center_y
208 namelist / param_mapprojection / &
211 mapprojection_basepoint_x, &
212 mapprojection_basepoint_y, &
213 mapprojection_type, &
214 mapprojection_rotation, &
215 mapprojection_lc_lat1, &
216 mapprojection_lc_lat2, &
217 mapprojection_ps_lat, &
218 mapprojection_m_lat, &
225 log_info(
"MAPPROJECTION_setup",*)
'Setup'
227 pi = real(pi_rp, kind=
dp)
228 d2r = real(d2r_rp, kind=
dp)
229 radius = real(radius_rp, kind=
dp)
231 mapprojection_basepoint_x = domain_center_x
232 mapprojection_basepoint_y = domain_center_y
233 mapprojection_ps_lat = undef
238 read(
io_fid_conf,nml=param_mapprojection,iostat=ierr)
241 log_info(
"MAPPROJECTION_setup",*)
'Not found namelist. Default used.'
242 elseif( ierr > 0 )
then
243 log_error(
"MAPPROJECTION_setup",*)
'Not appropriate names in namelist PARAM_MAPPROJECTION. Check!'
246 log_nml(param_mapprojection)
248 if (
prc_twod .and. mapprojection_type .ne.
'NONE' )
then
249 log_error(
"MAPPROJECTION_setup",*)
'MAPPROJECTION_type must be "NONE" for 2D experiment'
254 log_info(
"MAPPROJECTION_setup",*)
'Map projection information'
255 log_info_cont(
'(1x,A,F15.3)')
'Basepoint(x) [m] : ', mapprojection_basepoint_x
256 log_info_cont(
'(1x,A,F15.3)')
'Basepoint(y) [m] : ', mapprojection_basepoint_y
257 log_info_cont(*)
'Map projection type : ', trim(mapprojection_type)
270 select case(mapprojection_type)
272 log_info_cont(*)
'=> NO map projection'
280 log_info_cont(*)
'=> Lambert Conformal projection'
290 log_info_cont(*)
'=> Polar Stereographic projection'
301 log_info_cont(*)
'=> Mercator projection'
310 log_info_cont(*)
'=> Equidistant Cylindrical projection'
320 log_error(
"MAPPROJECTION_setup",*)
'Unsupported MAPPROJECTION_type. STOP'
325 mapprojection_mappingparam )
339 select case( info%mapping_name )
343 param%basepoint_lon = info%longitude_of_projection_origin * d2r
344 case(
"lambert_conformal_conic" )
347 param%basepoint_lon = info%longitude_of_central_meridian * d2r
348 case(
"polar_stereographic" )
351 param%basepoint_lon = info%straight_vertical_longitude_from_pole * d2r
355 param%basepoint_lon = info%longitude_of_projection_origin * d2r
356 case(
"equirectangular" )
359 param%basepoint_lon = info%longitude_of_central_meridian * d2r
361 log_error(
"MAPPROJECTION_set_param",*)
'Unsupported mapping type. STOP'
365 param%basepoint_x = info%false_easting
366 param%basepoint_y = info%false_northing
368 param%basepoint_lat = info%latitude_of_projection_origin * d2r
370 param%rotation = info%rotation * d2r
371 param%rot_fact_sin = sin(param%rotation)
372 param%rot_fact_cos = cos(param%rotation)
383 real(RP),
intent(in) :: x, y
384 real(RP),
intent(out) :: lon, lat
388 mapprojection_mappingparam, &
394 subroutine mapprojection_xy2lonlat_2d_initialized( &
395 IA, IS, IE, JA, JS, JE, &
399 integer,
intent(in) :: IA, IS, IE
400 integer,
intent(in) :: JA, JS, JE
402 real(RP),
intent(in) :: x(IA,JA)
403 real(RP),
intent(in) :: y(IA,JA)
404 real(RP),
intent(out) :: lon(IA,JA)
405 real(RP),
intent(out) :: lat(IA,JA)
418 end subroutine mapprojection_xy2lonlat_2d_initialized
420 subroutine mapprojection_xy2lonlat_0d_param( &
428 real(RP),
intent(in) :: x, y
429 character(len=*),
intent(in) :: mapping_name
432 real(RP),
intent(out) :: lon, lat
435 select case( mapping_name )
440 case(
"lambert_conformal_conic" )
444 case(
"polar_stereographic" )
452 case(
"equirectangular" )
457 log_error(
"MAPPROJECTION_xy2lonlat_0D_param",*)
'Unsupported mapping type. STOP: ', trim(mapping_name)
462 end subroutine mapprojection_xy2lonlat_0d_param
465 IA, IS, IE, JA, JS, JE, &
473 integer,
intent(in) :: IA, IS, IE
474 integer,
intent(in) :: JA, JS, JE
476 real(RP),
intent(in) :: x(IA,JA)
477 real(RP),
intent(in) :: y(IA,JA)
478 character(len=*),
intent(in) :: mapping_name
481 real(RP),
intent(out) :: lon(IA,JA)
482 real(RP),
intent(out) :: lat(IA,JA)
487 select case( mapping_name )
497 case(
"lambert_conformal_conic" )
506 case(
"polar_stereographic" )
524 case(
"equirectangular" )
534 log_error(
"MAPPROJECTION_xy2lonlat_2D_param",*)
'Unsupported mapping type. STOP: ', trim(mapping_name)
547 real(RP),
intent(in) :: lon, lat
548 real(RP),
intent(out) :: x, y
551 call lonlat2xy( lon, lat, &
552 mapprojection_mappingparam, &
558 subroutine mapprojection_lonlat2xy_2d_initialized( &
559 IA, IS, IE, JA, JS, JE, &
563 integer,
intent(in) :: IA, IS, IE
564 integer,
intent(in) :: JA, JS, JE
566 real(RP),
intent(in) :: lon(IA,JA)
567 real(RP),
intent(in) :: lat(IA,JA)
568 real(RP),
intent(out) :: x(IA,JA)
569 real(RP),
intent(out) :: y(IA,JA)
582 end subroutine mapprojection_lonlat2xy_2d_initialized
584 subroutine mapprojection_lonlat2xy_0d_param( &
592 real(RP),
intent(in) :: lon, lat
593 character(len=*),
intent(in) :: mapping_name
596 real(RP),
intent(out) :: x, y
599 select case( mapping_name )
604 case(
"lambert_conformal_conic" )
608 case(
"polar_stereographic" )
616 case(
"equirectangular" )
621 log_error(
"MAPPROJECTION_lonlat2xy_0D_param",*)
'Unsupported mapping type. STOP: ', trim(mapping_name)
626 end subroutine mapprojection_lonlat2xy_0d_param
629 IA, IS, IE, JA, JS, JE, &
637 integer,
intent(in) :: IA, IS, IE
638 integer,
intent(in) :: JA, JS, JE
640 real(RP),
intent(in) :: lon(IA,JA)
641 real(RP),
intent(in) :: lat(IA,JA)
642 character(len=*),
intent(in) :: mapping_name
645 real(RP),
intent(out) :: x(IA,JA)
646 real(RP),
intent(out) :: y(IA,JA)
651 select case( mapping_name )
661 case(
"lambert_conformal_conic" )
670 case(
"polar_stereographic" )
688 case(
"equirectangular" )
698 log_error(
"MAPPROJECTION_lonlat2xy_2D_param",*)
'Unsupported mapping type. STOP: ', trim(mapping_name)
708 IA, IS, IE, JA, JS, JE, &
712 integer,
intent(in) :: IA, IS, IE
713 integer,
intent(in) :: JA, JS, JE
715 real(RP),
intent(in) :: lat(IA,JA)
716 real(RP),
intent(out) :: m1 (IA,JA)
717 real(RP),
intent(out) :: m2 (IA,JA)
725 call mapfactor( lat(i,j), &
726 mapprojection_mappingparam, &
734 subroutine mapprojection_mapfactor_param( &
735 IA, IS, IE, JA, JS, JE, &
743 integer,
intent(in) :: IA, IS, IE
744 integer,
intent(in) :: JA, JS, JE
746 real(RP),
intent(in) :: lat(IA,JA)
747 character(len=*),
intent(in) :: mapping_name
750 real(RP),
intent(out) :: m1 (IA,JA)
751 real(RP),
intent(out) :: m2 (IA,JA)
756 select case( mapping_name )
766 case(
"lambert_conformal_conic" )
775 case(
"polar_stereographic" )
793 case(
"equirectangular" )
803 log_error(
"MAPPROJECTION_mapfactor_param",*)
'Unsupported mapping type. STOP: ', trim(mapping_name)
808 end subroutine mapprojection_mapfactor_param
814 IA, IS, IE, JA, JS, JE, &
818 integer,
intent(in) :: IA, IS, IE
819 integer,
intent(in) :: JA, JS, JE
821 real(RP),
intent(in) :: lon(IA,JA)
822 real(RP),
intent(in) :: lat(IA,JA)
823 real(RP),
intent(out) :: rotc_cos(IA,JA)
824 real(RP),
intent(out) :: rotc_sin(IA,JA)
832 call rotcoef( lon(i,j), lat(i,j), &
833 mapprojection_mappingparam, &
834 rotc_cos(i,j), rotc_sin(i,j) )
841 subroutine mapprojection_rotcoef_param( &
842 IA, IS, IE, JA, JS, JE, &
850 integer,
intent(in) :: IA, IS, IE
851 integer,
intent(in) :: JA, JS, JE
853 real(RP),
intent(in) :: lon(IA,JA)
854 real(RP),
intent(in) :: lat(IA,JA)
855 character(len=*),
intent(in) :: mapping_name
858 real(RP),
intent(out) :: rotc_cos(IA,JA)
859 real(RP),
intent(out) :: rotc_sin(IA,JA)
864 select case( mapping_name )
871 rotc_cos(i,j), rotc_sin(i,j) )
874 case(
"lambert_conformal_conic" )
880 rotc_cos(i,j), rotc_sin(i,j) )
883 case(
"polar_stereographic" )
889 rotc_cos(i,j), rotc_sin(i,j) )
898 rotc_cos(i,j), rotc_sin(i,j) )
901 case(
"equirectangular" )
907 rotc_cos(i,j), rotc_sin(i,j) )
911 log_error(
"MAPPROJECTION_rotcoef_param",*)
'Unsupported mapping type. STOP: ', trim(mapping_name)
916 end subroutine mapprojection_rotcoef_param
936 real(
rp),
intent(in) :: x, y
938 real(
rp),
intent(out) :: lon, lat
941 real(
dp) :: gno1, gno2
943 real(
dp) :: gmmc, gmms
944 real(
dp) :: latc, lats
947 if ( x==undef .or. y==undef )
then
954 call xy_rotate( x, y, &
958 gno1 = ( x - param%basepoint_x ) / radius
959 gno2 = ( y - param%basepoint_y ) / radius
961 rho = sqrt( gno1 * gno1 + gno2 * gno2 )
964 if ( rho == 0.0_dp )
then
965 lon = param%basepoint_lon
966 lat = param%basepoint_lat
970 latc = cos(param%basepoint_lat)
971 lats = sin(param%basepoint_lat)
972 lon = param%basepoint_lon &
973 + atan( gno1*gmms / ( rho * latc * gmmc &
974 - gno2 * lats * gmms ) )
975 lat = asin( lats * gmmc &
976 + gno2*latc * gmms / rho )
989 real(
rp),
intent(in) :: lon, lat
991 real(
rp),
intent(out) :: x, y
994 real(
dp) :: lat_d, lat0_d, dlon
995 real(
dp) :: lat_d_c, lat_d_s, lat0_d_c, lat0_d_s, dlon_c, dlon_s
996 real(
dp) :: gno1, gno2
1001 if ( lon==undef .or. lat==undef )
then
1008 lat_d = real(lat,kind=
dp)
1009 lat0_d = param%basepoint_lat
1011 dlon = lon - param%basepoint_lon
1013 lat_d_c = cos(lat_d)
1014 lat_d_s = sin(lat_d)
1015 lat0_d_c = cos(lat0_d)
1016 lat0_d_s = sin(lat0_d)
1020 cos_gmm = lat0_d_s * lat_d_s &
1021 + lat0_d_c * lat_d_c * dlon_c
1023 gno1 = ( lat_d_c * dlon_s ) / cos_gmm
1024 gno2 = ( lat0_d_c * lat_d_s &
1025 - lat0_d_s * lat_d_c * dlon_c ) / cos_gmm
1027 xx = gno1 * radius + param%basepoint_x
1028 yy = gno2 * radius + param%basepoint_y
1030 call xy_unrotate( xx, yy, &
1046 real(
rp),
intent(in) :: lat
1048 real(
rp),
intent(out) :: m1, m2
1062 rotc_cos, rotc_sin )
1064 real(
rp),
intent(in) :: lon, lat
1066 real(
rp),
intent(out) :: rotc_cos, rotc_sin
1069 rotc_cos = param%rot_fact_cos
1070 rotc_sin = param%rot_fact_sin
1086 real(
dp) :: lc_lat1, lc_lat2
1087 real(
dp) :: basepoint_lat
1088 real(
dp) :: basepoint_x, basepoint_y
1090 real(
dp) :: lat1rot, lat2rot
1091 real(
dp) :: dlon, latrot, dist
1094 lc_lat1 = info%standard_parallel(1)
1095 lc_lat2 = info%standard_parallel(2)
1097 basepoint_lat = info%latitude_of_projection_origin
1098 basepoint_x = info%false_easting
1099 basepoint_y = info%false_northing
1101 if ( lc_lat1 >= lc_lat2 )
then
1102 log_error(
"MAPPROJECTION_get_param_LambertConformal",*)
'Please set LC_lat1 < LC_lat2 in degree. STOP'
1107 param%hemisphere = sign(1.0_dp,lc_lat1+lc_lat2)
1109 lat1rot = 0.5_dp*pi - param%hemisphere * lc_lat1 * d2r
1110 lat2rot = 0.5_dp*pi - param%hemisphere * lc_lat2 * d2r
1113 param%c = ( log( sin(lat1rot) ) - log( sin(lat2rot) ) ) &
1114 / ( log( tan(0.5_dp*lat1rot) ) - log( tan(0.5_dp*lat2rot) ) )
1117 param%fact = sin(lat1rot) / param%c / tan(0.5_dp*lat1rot)**param%c
1122 latrot = 0.5_dp*pi - param%hemisphere * basepoint_lat * d2r
1124 dist = param%fact * radius * tan(0.5_dp*latrot)**param%c
1126 param%x = basepoint_x - dist * sin(param%c*dlon)
1127 param%y = basepoint_y + param%hemisphere * dist * cos(param%c*dlon)
1130 log_info(
"MAPPROJECTION_get_param_LambertConformal",*)
'Input parameters'
1131 log_info_cont(*)
'LC_lat1 = ', lc_lat1
1132 log_info_cont(*)
'LC_lat2 = ', lc_lat2
1133 log_info_cont(*)
'hemisphere = ', param%hemisphere
1134 log_info_cont(*)
'LC_c = ', param%c
1135 log_info_cont(*)
'LC_fact = ', param%fact
1136 log_info_cont(*)
'pole_x = ', param%x
1137 log_info_cont(*)
'pole_y = ', param%y
1151 real(
rp),
intent(in) :: x, y
1153 real(
rp),
intent(out) :: lon, lat
1155 real(
dp) :: xx, yy, dist
1158 if ( x==undef .or. y==undef )
then
1165 call xy_rotate( x, y, &
1168 xx = ( xx - param%x ) / radius / param%fact
1169 yy = - param%hemisphere * ( yy - param%y ) / radius / param%fact
1171 dist = sqrt( xx*xx + yy*yy )
1173 lon = param%basepoint_lon + atan2(xx,yy) / param%c
1176 lat = param%hemisphere * ( 0.5_dp*pi - 2.0_dp*atan( dist**(1.0_dp/param%c) ) )
1190 real(
rp),
intent(in) :: lon, lat
1192 real(
rp),
intent(out) :: x, y
1195 real(
dp) :: dlon, latrot, dist
1198 if ( lon==undef .or. lat==undef )
then
1205 dlon = lon - param%basepoint_lon
1206 if ( dlon > pi ) dlon = dlon - pi*2.0_dp
1207 if ( dlon < -pi ) dlon = dlon + pi*2.0_dp
1209 latrot = 0.5_dp*pi - param%hemisphere * lat
1211 dist = param%fact * radius * tan(0.5_dp*latrot)**param%c
1213 xx = param%x + dist * sin(param%c*dlon)
1214 yy = param%y - param%hemisphere * dist * cos(param%c*dlon)
1216 call xy_unrotate( xx, yy, &
1230 real(
rp),
intent(in) :: lat
1232 real(
rp),
intent(out) :: m1, m2
1237 latrot = 0.5_dp*pi - param%hemisphere * lat
1239 m1 = param%fact / sin(latrot) * param%c * tan(0.5_dp*latrot)**param%c
1249 rotc_cos, rotc_sin )
1251 real(
rp),
intent(in) :: lon, lat
1253 real(
rp),
intent(out) :: rotc_cos, rotc_sin
1259 dlon = lon - param%basepoint_lon
1260 if ( dlon > pi ) dlon = dlon - pi*2.0_dp
1261 if ( dlon < -pi ) dlon = dlon + pi*2.0_dp
1263 alpha = - param%c * dlon * param%hemisphere &
1266 rotc_cos = cos( alpha )
1267 rotc_sin = sin( alpha )
1282 real(
dp) :: basepoint_lat
1283 real(
dp) :: basepoint_x, basepoint_y
1286 real(
dp) :: dlon, latrot, dist
1289 ps_lat = info%standard_parallel(1)
1291 basepoint_lat = info%latitude_of_projection_origin
1292 basepoint_x = info%false_easting
1293 basepoint_y = info%false_northing
1296 param%hemisphere = sign(1.0_dp,ps_lat)
1298 lat0 = param%hemisphere * ps_lat * d2r
1301 param%fact = 1.0_dp + sin(lat0)
1306 latrot = 0.5_dp*pi - param%hemisphere * basepoint_lat * d2r
1308 dist = param%fact * radius * tan(0.5_dp*latrot)
1310 param%x = basepoint_x - dist * sin(dlon)
1311 param%y = basepoint_y + param%hemisphere * dist * cos(dlon)
1314 log_info(
"MAPPROJECTION_get_param_PolarStereographic",*)
'PS_lat1 = ', ps_lat
1315 log_info(
"MAPPROJECTION_get_param_PolarStereographic",*)
'hemisphere = ', param%hemisphere
1316 log_info(
"MAPPROJECTION_get_param_PolarStereographic",*)
'PS_fact = ', param%fact
1317 log_info(
"MAPPROJECTION_get_param_PolarStereographic",*)
'pole_x = ', param%x
1318 log_info(
"MAPPROJECTION_get_param_PolarStereographic",*)
'pole_y = ', param%y
1332 real(
rp),
intent(in) :: x, y
1334 real(
rp),
intent(out) :: lon, lat
1336 real(
dp) :: xx, yy, dist
1339 if ( x==undef .or. y==undef )
then
1346 call xy_rotate( x, y, &
1350 xx = ( xx - param%x ) / radius / param%fact
1351 yy = - param%hemisphere * ( yy - param%y ) / radius / param%fact
1353 dist = sqrt( xx*xx + yy*yy )
1355 lon = param%basepoint_lon + atan2(xx,yy)
1356 lat = param%hemisphere * ( 0.5_dp*pi - 2.0_dp*atan(dist) )
1370 real(
rp),
intent(in) :: lon, lat
1372 real(
rp),
intent(out) :: x, y
1374 real(
dp) :: dlon, latrot, dist
1378 if ( lon==undef .or. lat==undef )
then
1385 dlon = lon - param%basepoint_lon
1387 latrot = 0.5_dp*pi - param%hemisphere * lat
1389 dist = param%fact * radius * tan(0.5_dp*latrot)
1391 xx = param%x + dist * sin(dlon)
1392 yy = param%y - param%hemisphere * dist * cos(dlon)
1394 call xy_unrotate( xx, yy, &
1408 real(
rp),
intent(in) :: lat
1410 real(
rp),
intent(out) :: m1, m2
1413 m1 = param%fact / ( 1.0_dp + sin(param%hemisphere*lat) )
1423 rotc_cos, rotc_sin )
1425 real(
rp),
intent(in) :: lon, lat
1427 real(
rp),
intent(out) :: rotc_cos, rotc_sin
1433 dlon = lon - param%basepoint_lon
1434 if ( dlon > pi ) dlon = dlon - pi*2.0_dp
1435 if ( dlon < -pi ) dlon = dlon + pi*2.0_dp
1437 alpha = - dlon * param%hemisphere &
1440 rotc_cos = cos( alpha )
1441 rotc_sin = sin( alpha )
1458 real(
dp) :: basepoint_lat
1459 real(
dp) :: basepoint_x, basepoint_y
1462 real(
dp) :: latrot, dist
1465 m_lat = info%standard_parallel(1)
1467 basepoint_lat = info%latitude_of_projection_origin
1468 basepoint_x = info%false_easting
1469 basepoint_y = info%false_northing
1474 param%fact = cos(lat0)
1476 if ( param%fact == 0.0_dp )
then
1477 log_error(
"MAPPROJECTION_get_param_Mercator",*)
'M_lat cannot be set to pole point! value=', m_lat
1482 latrot = 0.5_dp*pi - basepoint_lat * d2r
1484 dist = 1.0_dp / tan(0.5_dp*latrot)
1486 param%x = basepoint_x
1487 param%y = basepoint_y - radius * param%fact * log(dist)
1490 log_info(
"MAPPROJECTION_get_param_Mercator",*)
'M_lat = ', m_lat
1491 log_info(
"MAPPROJECTION_get_param_Mercator",*)
'M_fact = ', param%fact
1492 log_info(
"MAPPROJECTION_get_param_Mercator",*)
'eq_x = ', param%x
1493 log_info(
"MAPPROJECTION_get_param_Mercator",*)
'eq_y = ', param%y
1507 real(
rp),
intent(in) :: x, y
1509 real(
rp),
intent(out) :: lon, lat
1514 if ( x==undef .or. y==undef )
then
1521 call xy_rotate( x, y, &
1525 xx = ( xx - param%x ) / radius / param%fact
1526 yy = ( yy - param%y ) / radius / param%fact
1528 lon = xx + param%basepoint_lon
1530 lat = 0.5_dp*pi - 2.0_dp*atan( 1.0_dp/exp(yy) )
1544 real(
rp),
intent(in) :: lon, lat
1546 real(
rp),
intent(out) :: x, y
1548 real(
dp) :: dlon, latrot, dist
1551 if ( lon==undef .or. lat==undef )
then
1558 dlon = lon - param%basepoint_lon
1560 latrot = 0.5_dp*pi - lat
1562 dist = 1.0_dp / tan(0.5_dp*latrot)
1564 x = param%x + radius * param%fact * dlon
1565 y = param%y + radius * param%fact * log(dist)
1577 real(
rp),
intent(in) :: lat
1579 real(
rp),
intent(out) :: m1, m2
1582 m1 = param%fact / cos(real(lat,kind=
dp))
1592 rotc_cos, rotc_sin )
1594 real(
rp),
intent(in) :: lon, lat
1596 real(
rp),
intent(out) :: rotc_cos, rotc_sin
1599 rotc_cos = param%rot_fact_cos
1600 rotc_sin = param%rot_fact_sin
1617 real(
dp) :: basepoint_lat
1618 real(
dp) :: basepoint_x, basepoint_y
1623 ec_lat = info%standard_parallel(1)
1625 basepoint_lat = info%latitude_of_projection_origin
1626 basepoint_x = info%false_easting
1627 basepoint_y = info%false_northing
1632 param%fact = cos(lat0)
1634 if ( param%fact == 0.0_dp )
then
1635 log_error(
"MAPPROJECTION_get_param_EquidistantCylindrical",*)
'EC_lat cannot be set to pole point! value=', ec_lat
1639 param%x = basepoint_x
1640 param%y = basepoint_y - radius * basepoint_lat * d2r
1643 log_info(
"MAPPROJECTION_get_param_EquidistantCylindrical",*)
'EC_lat = ', ec_lat
1644 log_info(
"MAPPROJECTION_get_param_EquidistantCylindrical",*)
'EC_fact = ', param%fact
1645 log_info(
"MAPPROJECTION_get_param_EquidistantCylindrical",*)
'eq_x = ', param%x
1646 log_info(
"MAPPROJECTION_get_param_EquidistantCylindrical",*)
'eq_y = ', param%y
1662 real(
rp),
intent(in) :: x, y
1664 real(
rp),
intent(out) :: lon, lat
1669 if ( x==undef .or. y==undef )
then
1676 call xy_rotate( x, y, &
1680 xx = ( xx - param%x ) / radius / param%fact
1681 yy = ( yy - param%y ) / radius
1683 lon = xx + param%basepoint_lon
1686 if ( abs(lat) > 0.5_dp*pi )
then
1687 log_error(
"MAPPROJECTION_EquidistantCylindrical_xy2lonlat",*)
'Invalid latitude range! value=', lat
1703 real(
rp),
intent(in) :: lon, lat
1705 real(
rp),
intent(out) :: x, y
1711 if ( lon==undef .or. lat==undef )
then
1718 dlon = lon - param%basepoint_lon
1720 xx = param%x + radius * param%fact * dlon
1721 yy = param%y + radius * lat
1723 call xy_unrotate( xx, yy, &
1737 real(
rp),
intent(in) :: lat
1739 real(
rp),
intent(out) :: m1, m2
1741 real(
dp) :: mm1, mm2
1744 mm1 = param%fact / cos(real(lat,kind=
dp))
1747 m1 = sqrt( (param%rot_fact_cos * mm1)**2 + (param%rot_fact_sin * mm2)**2 )
1748 m2 = sqrt( (param%rot_fact_cos * mm2)**2 + (param%rot_fact_sin * mm1)**2 )
1757 rotc_cos, rotc_sin )
1759 real(
rp),
intent(in) :: lon, lat
1761 real(
rp),
intent(out) :: rotc_cos, rotc_sin
1764 rotc_cos = param%rot_fact_cos
1765 rotc_sin = param%rot_fact_sin
1772 subroutine xy_rotate( &
1777 real(
rp),
intent(in) :: x, y
1779 real(
dp),
intent(out) :: xx, yy
1783 xd = x - param%basepoint_x
1784 yd = y - param%basepoint_y
1786 xx = param%basepoint_x + xd * param%rot_fact_cos &
1787 - yd * param%rot_fact_sin
1788 yy = param%basepoint_y + yd * param%rot_fact_cos &
1789 + xd * param%rot_fact_sin
1792 end subroutine xy_rotate
1794 subroutine xy_unrotate( &
1799 real(
dp),
intent(in) :: xx, yy
1801 real(
rp),
intent(out) :: x, y
1803 real(
dp) :: xxd, yyd
1805 xxd = xx - param%basepoint_x
1806 yyd = yy - param%basepoint_y
1808 x = param%basepoint_x + ( xxd * param%rot_fact_cos &
1809 + yyd * param%rot_fact_sin )
1810 y = param%basepoint_y + ( yyd * param%rot_fact_cos &
1811 - xxd * param%rot_fact_sin )
1814 end subroutine xy_unrotate