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
146 real(
dp),
private :: mapprojection_ec_lat
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
234 mapprojection_m_lat = undef
235 mapprojection_ec_lat = undef
240 read(
io_fid_conf,nml=param_mapprojection,iostat=ierr)
243 log_info(
"MAPPROJECTION_setup",*)
'Not found namelist. Default used.'
244 elseif( ierr > 0 )
then
245 log_error(
"MAPPROJECTION_setup",*)
'Not appropriate names in namelist PARAM_MAPPROJECTION. Check!'
248 log_nml(param_mapprojection)
250 if (
prc_twod .and. mapprojection_type .ne.
'NONE' )
then
251 log_error(
"MAPPROJECTION_setup",*)
'MAPPROJECTION_type must be "NONE" for 2D experiment'
256 log_info(
"MAPPROJECTION_setup",*)
'Map projection information'
257 log_info_cont(
'(1x,A,F15.3)')
'Basepoint(x) [m] : ', mapprojection_basepoint_x
258 log_info_cont(
'(1x,A,F15.3)')
'Basepoint(y) [m] : ', mapprojection_basepoint_y
259 log_info_cont(*)
'Map projection type : ', trim(mapprojection_type)
272 select case(mapprojection_type)
274 log_info_cont(*)
'=> NO map projection'
282 log_info_cont(*)
'=> Lambert Conformal projection'
292 log_info_cont(*)
'=> Polar Stereographic projection'
303 log_info_cont(*)
'=> Mercator projection'
313 log_info_cont(*)
'=> Equidistant Cylindrical projection'
324 log_error(
"MAPPROJECTION_setup",*)
'Unsupported MAPPROJECTION_type. STOP'
329 mapprojection_mappingparam )
343 select case( info%mapping_name )
347 param%basepoint_lon = info%longitude_of_projection_origin * d2r
348 case(
"lambert_conformal_conic" )
351 param%basepoint_lon = info%longitude_of_central_meridian * d2r
352 case(
"polar_stereographic" )
355 param%basepoint_lon = info%straight_vertical_longitude_from_pole * d2r
359 param%basepoint_lon = info%longitude_of_projection_origin * d2r
360 case(
"equirectangular" )
363 param%basepoint_lon = info%longitude_of_central_meridian * d2r
365 log_error(
"MAPPROJECTION_set_param",*)
'Unsupported mapping type. STOP'
369 param%basepoint_x = info%false_easting
370 param%basepoint_y = info%false_northing
372 param%basepoint_lat = info%latitude_of_projection_origin * d2r
374 param%rotation = info%rotation * d2r
375 param%rot_fact_sin = sin(param%rotation)
376 param%rot_fact_cos = cos(param%rotation)
387 real(RP),
intent(in) :: x, y
388 real(RP),
intent(out) :: lon, lat
392 mapprojection_mappingparam, &
398 subroutine mapprojection_xy2lonlat_2d_initialized( &
399 IA, IS, IE, JA, JS, JE, &
403 integer,
intent(in) :: IA, IS, IE
404 integer,
intent(in) :: JA, JS, JE
406 real(RP),
intent(in) :: x(IA,JA)
407 real(RP),
intent(in) :: y(IA,JA)
408 real(RP),
intent(out) :: lon(IA,JA)
409 real(RP),
intent(out) :: lat(IA,JA)
422 end subroutine mapprojection_xy2lonlat_2d_initialized
424 subroutine mapprojection_xy2lonlat_0d_param( &
432 real(RP),
intent(in) :: x, y
433 character(len=*),
intent(in) :: mapping_name
436 real(RP),
intent(out) :: lon, lat
439 select case( mapping_name )
444 case(
"lambert_conformal_conic" )
448 case(
"polar_stereographic" )
456 case(
"equirectangular" )
461 log_error(
"MAPPROJECTION_xy2lonlat_0D_param",*)
'Unsupported mapping type. STOP: ', trim(mapping_name)
466 end subroutine mapprojection_xy2lonlat_0d_param
469 IA, IS, IE, JA, JS, JE, &
477 integer,
intent(in) :: IA, IS, IE
478 integer,
intent(in) :: JA, JS, JE
480 real(RP),
intent(in) :: x(IA,JA)
481 real(RP),
intent(in) :: y(IA,JA)
482 character(len=*),
intent(in) :: mapping_name
485 real(RP),
intent(out) :: lon(IA,JA)
486 real(RP),
intent(out) :: lat(IA,JA)
491 select case( mapping_name )
501 case(
"lambert_conformal_conic" )
510 case(
"polar_stereographic" )
528 case(
"equirectangular" )
538 log_error(
"MAPPROJECTION_xy2lonlat_2D_param",*)
'Unsupported mapping type. STOP: ', trim(mapping_name)
551 real(RP),
intent(in) :: lon, lat
552 real(RP),
intent(out) :: x, y
555 call lonlat2xy( lon, lat, &
556 mapprojection_mappingparam, &
562 subroutine mapprojection_lonlat2xy_2d_initialized( &
563 IA, IS, IE, JA, JS, JE, &
567 integer,
intent(in) :: IA, IS, IE
568 integer,
intent(in) :: JA, JS, JE
570 real(RP),
intent(in) :: lon(IA,JA)
571 real(RP),
intent(in) :: lat(IA,JA)
572 real(RP),
intent(out) :: x(IA,JA)
573 real(RP),
intent(out) :: y(IA,JA)
586 end subroutine mapprojection_lonlat2xy_2d_initialized
588 subroutine mapprojection_lonlat2xy_0d_param( &
596 real(RP),
intent(in) :: lon, lat
597 character(len=*),
intent(in) :: mapping_name
600 real(RP),
intent(out) :: x, y
603 select case( mapping_name )
608 case(
"lambert_conformal_conic" )
612 case(
"polar_stereographic" )
620 case(
"equirectangular" )
625 log_error(
"MAPPROJECTION_lonlat2xy_0D_param",*)
'Unsupported mapping type. STOP: ', trim(mapping_name)
630 end subroutine mapprojection_lonlat2xy_0d_param
633 IA, IS, IE, JA, JS, JE, &
641 integer,
intent(in) :: IA, IS, IE
642 integer,
intent(in) :: JA, JS, JE
644 real(RP),
intent(in) :: lon(IA,JA)
645 real(RP),
intent(in) :: lat(IA,JA)
646 character(len=*),
intent(in) :: mapping_name
649 real(RP),
intent(out) :: x(IA,JA)
650 real(RP),
intent(out) :: y(IA,JA)
655 select case( mapping_name )
665 case(
"lambert_conformal_conic" )
674 case(
"polar_stereographic" )
692 case(
"equirectangular" )
702 log_error(
"MAPPROJECTION_lonlat2xy_2D_param",*)
'Unsupported mapping type. STOP: ', trim(mapping_name)
712 IA, IS, IE, JA, JS, JE, &
716 integer,
intent(in) :: IA, IS, IE
717 integer,
intent(in) :: JA, JS, JE
719 real(RP),
intent(in) :: lat(IA,JA)
720 real(RP),
intent(out) :: m1 (IA,JA)
721 real(RP),
intent(out) :: m2 (IA,JA)
729 call mapfactor( lat(i,j), &
730 mapprojection_mappingparam, &
738 subroutine mapprojection_mapfactor_param( &
739 IA, IS, IE, JA, JS, JE, &
747 integer,
intent(in) :: IA, IS, IE
748 integer,
intent(in) :: JA, JS, JE
750 real(RP),
intent(in) :: lat(IA,JA)
751 character(len=*),
intent(in) :: mapping_name
754 real(RP),
intent(out) :: m1 (IA,JA)
755 real(RP),
intent(out) :: m2 (IA,JA)
760 select case( mapping_name )
770 case(
"lambert_conformal_conic" )
779 case(
"polar_stereographic" )
797 case(
"equirectangular" )
807 log_error(
"MAPPROJECTION_mapfactor_param",*)
'Unsupported mapping type. STOP: ', trim(mapping_name)
812 end subroutine mapprojection_mapfactor_param
818 IA, IS, IE, JA, JS, JE, &
822 integer,
intent(in) :: IA, IS, IE
823 integer,
intent(in) :: JA, JS, JE
825 real(RP),
intent(in) :: lon(IA,JA)
826 real(RP),
intent(in) :: lat(IA,JA)
827 real(RP),
intent(out) :: rotc_cos(IA,JA)
828 real(RP),
intent(out) :: rotc_sin(IA,JA)
836 call rotcoef( lon(i,j), lat(i,j), &
837 mapprojection_mappingparam, &
838 rotc_cos(i,j), rotc_sin(i,j) )
845 subroutine mapprojection_rotcoef_param( &
846 IA, IS, IE, JA, JS, JE, &
854 integer,
intent(in) :: IA, IS, IE
855 integer,
intent(in) :: JA, JS, JE
857 real(RP),
intent(in) :: lon(IA,JA)
858 real(RP),
intent(in) :: lat(IA,JA)
859 character(len=*),
intent(in) :: mapping_name
862 real(RP),
intent(out) :: rotc_cos(IA,JA)
863 real(RP),
intent(out) :: rotc_sin(IA,JA)
868 select case( mapping_name )
875 rotc_cos(i,j), rotc_sin(i,j) )
878 case(
"lambert_conformal_conic" )
884 rotc_cos(i,j), rotc_sin(i,j) )
887 case(
"polar_stereographic" )
893 rotc_cos(i,j), rotc_sin(i,j) )
902 rotc_cos(i,j), rotc_sin(i,j) )
905 case(
"equirectangular" )
911 rotc_cos(i,j), rotc_sin(i,j) )
915 log_error(
"MAPPROJECTION_rotcoef_param",*)
'Unsupported mapping type. STOP: ', trim(mapping_name)
920 end subroutine mapprojection_rotcoef_param
940 real(
rp),
intent(in) :: x, y
942 real(
rp),
intent(out) :: lon, lat
945 real(
dp) :: gno1, gno2
947 real(
dp) :: gmmc, gmms
948 real(
dp) :: latc, lats
951 if ( x==undef .or. y==undef )
then
958 call xy_rotate( x, y, &
962 gno1 = ( x - param%basepoint_x ) / radius
963 gno2 = ( y - param%basepoint_y ) / radius
965 rho = sqrt( gno1 * gno1 + gno2 * gno2 )
968 if ( rho == 0.0_dp )
then
969 lon = param%basepoint_lon
970 lat = param%basepoint_lat
974 latc = cos(param%basepoint_lat)
975 lats = sin(param%basepoint_lat)
976 lon = param%basepoint_lon &
977 + atan( gno1*gmms / ( rho * latc * gmmc &
978 - gno2 * lats * gmms ) )
979 lat = asin( lats * gmmc &
980 + gno2*latc * gmms / rho )
993 real(
rp),
intent(in) :: lon, lat
995 real(
rp),
intent(out) :: x, y
998 real(
dp) :: lat_d, lat0_d, dlon
999 real(
dp) :: lat_d_c, lat_d_s, lat0_d_c, lat0_d_s, dlon_c, dlon_s
1000 real(
dp) :: gno1, gno2
1005 if ( lon==undef .or. lat==undef )
then
1012 lat_d = real(lat,kind=
dp)
1013 lat0_d = param%basepoint_lat
1015 dlon = lon - param%basepoint_lon
1017 lat_d_c = cos(lat_d)
1018 lat_d_s = sin(lat_d)
1019 lat0_d_c = cos(lat0_d)
1020 lat0_d_s = sin(lat0_d)
1024 cos_gmm = lat0_d_s * lat_d_s &
1025 + lat0_d_c * lat_d_c * dlon_c
1027 gno1 = ( lat_d_c * dlon_s ) / cos_gmm
1028 gno2 = ( lat0_d_c * lat_d_s &
1029 - lat0_d_s * lat_d_c * dlon_c ) / cos_gmm
1031 xx = gno1 * radius + param%basepoint_x
1032 yy = gno2 * radius + param%basepoint_y
1034 call xy_unrotate( xx, yy, &
1050 real(
rp),
intent(in) :: lat
1052 real(
rp),
intent(out) :: m1, m2
1066 rotc_cos, rotc_sin )
1068 real(
rp),
intent(in) :: lon, lat
1070 real(
rp),
intent(out) :: rotc_cos, rotc_sin
1073 rotc_cos = param%rot_fact_cos
1074 rotc_sin = param%rot_fact_sin
1090 real(
dp),
parameter :: eps = 1d-16
1091 real(
dp) :: lc_lat1, lc_lat2
1092 real(
dp) :: basepoint_lat
1093 real(
dp) :: basepoint_x, basepoint_y
1095 real(
dp) :: lat1rot, lat2rot
1096 real(
dp) :: dlon, latrot, dist
1099 lc_lat1 = info%standard_parallel(1)
1100 lc_lat2 = info%standard_parallel(2)
1102 basepoint_lat = info%latitude_of_projection_origin
1103 basepoint_x = info%false_easting
1104 basepoint_y = info%false_northing
1106 if ( lc_lat1 > lc_lat2 )
then
1107 log_error(
"MAPPROJECTION_get_param_LambertConformal",*)
'Please set LC_lat1 <= LC_lat2 in degree. STOP'
1112 param%hemisphere = sign(1.0_dp,lc_lat1+lc_lat2)
1114 lat1rot = 0.5_dp*pi - param%hemisphere * lc_lat1 * d2r
1115 lat2rot = 0.5_dp*pi - param%hemisphere * lc_lat2 * d2r
1116 if ( lc_lat2 - lc_lat1 > eps )
then
1118 param%c = ( log( sin(lat1rot) ) - log( sin(lat2rot) ) ) &
1119 / ( log( tan(0.5_dp*lat1rot) ) - log( tan(0.5_dp*lat2rot) ) )
1121 param%c = cos(lat1rot)
1124 param%fact = sin(lat1rot) / param%c / tan(0.5_dp*lat1rot)**param%c
1129 latrot = 0.5_dp*pi - param%hemisphere * basepoint_lat * d2r
1131 dist = param%fact * radius * tan(0.5_dp*latrot)**param%c
1133 param%x = basepoint_x - dist * sin(param%c*dlon)
1134 param%y = basepoint_y + param%hemisphere * dist * cos(param%c*dlon)
1137 log_info(
"MAPPROJECTION_get_param_LambertConformal",*)
'Input parameters'
1138 log_info_cont(*)
'LC_lat1 = ', lc_lat1
1139 log_info_cont(*)
'LC_lat2 = ', lc_lat2
1140 log_info_cont(*)
'hemisphere = ', param%hemisphere
1141 log_info_cont(*)
'LC_c = ', param%c
1142 log_info_cont(*)
'LC_fact = ', param%fact
1143 log_info_cont(*)
'pole_x = ', param%x
1144 log_info_cont(*)
'pole_y = ', param%y
1158 real(
rp),
intent(in) :: x, y
1160 real(
rp),
intent(out) :: lon, lat
1162 real(
dp) :: xx, yy, dist
1165 if ( x==undef .or. y==undef )
then
1172 call xy_rotate( x, y, &
1175 xx = ( xx - param%x ) / radius / param%fact
1176 yy = - param%hemisphere * ( yy - param%y ) / radius / param%fact
1178 dist = sqrt( xx*xx + yy*yy )
1180 lon = param%basepoint_lon + atan2(xx,yy) / param%c
1183 lat = param%hemisphere * ( 0.5_dp*pi - 2.0_dp*atan( dist**(1.0_dp/param%c) ) )
1197 real(
rp),
intent(in) :: lon, lat
1199 real(
rp),
intent(out) :: x, y
1202 real(
dp) :: dlon, latrot, dist
1205 if ( lon==undef .or. lat==undef )
then
1212 dlon = lon - param%basepoint_lon
1213 if ( dlon > pi ) dlon = dlon - pi*2.0_dp
1214 if ( dlon < -pi ) dlon = dlon + pi*2.0_dp
1216 latrot = 0.5_dp*pi - param%hemisphere * lat
1218 dist = param%fact * radius * tan(0.5_dp*latrot)**param%c
1220 xx = param%x + dist * sin(param%c*dlon)
1221 yy = param%y - param%hemisphere * dist * cos(param%c*dlon)
1223 call xy_unrotate( xx, yy, &
1237 real(
rp),
intent(in) :: lat
1239 real(
rp),
intent(out) :: m1, m2
1244 latrot = 0.5_dp*pi - param%hemisphere * lat
1246 m1 = param%fact / sin(latrot) * param%c * tan(0.5_dp*latrot)**param%c
1256 rotc_cos, rotc_sin )
1258 real(
rp),
intent(in) :: lon, lat
1260 real(
rp),
intent(out) :: rotc_cos, rotc_sin
1266 dlon = lon - param%basepoint_lon
1267 if ( dlon > pi ) dlon = dlon - pi*2.0_dp
1268 if ( dlon < -pi ) dlon = dlon + pi*2.0_dp
1270 alpha = - param%c * dlon * param%hemisphere &
1273 rotc_cos = cos( alpha )
1274 rotc_sin = sin( alpha )
1289 real(
dp) :: basepoint_lat
1290 real(
dp) :: basepoint_x, basepoint_y
1293 real(
dp) :: dlon, latrot, dist
1296 ps_lat = info%standard_parallel(1)
1298 basepoint_lat = info%latitude_of_projection_origin
1299 basepoint_x = info%false_easting
1300 basepoint_y = info%false_northing
1303 param%hemisphere = sign(1.0_dp,ps_lat)
1305 lat0 = param%hemisphere * ps_lat * d2r
1308 param%fact = 1.0_dp + sin(lat0)
1313 latrot = 0.5_dp*pi - param%hemisphere * basepoint_lat * d2r
1315 dist = param%fact * radius * tan(0.5_dp*latrot)
1317 param%x = basepoint_x - dist * sin(dlon)
1318 param%y = basepoint_y + param%hemisphere * dist * cos(dlon)
1321 log_info(
"MAPPROJECTION_get_param_PolarStereographic",*)
'PS_lat1 = ', ps_lat
1322 log_info(
"MAPPROJECTION_get_param_PolarStereographic",*)
'hemisphere = ', param%hemisphere
1323 log_info(
"MAPPROJECTION_get_param_PolarStereographic",*)
'PS_fact = ', param%fact
1324 log_info(
"MAPPROJECTION_get_param_PolarStereographic",*)
'pole_x = ', param%x
1325 log_info(
"MAPPROJECTION_get_param_PolarStereographic",*)
'pole_y = ', param%y
1339 real(
rp),
intent(in) :: x, y
1341 real(
rp),
intent(out) :: lon, lat
1343 real(
dp) :: xx, yy, dist
1346 if ( x==undef .or. y==undef )
then
1353 call xy_rotate( x, y, &
1357 xx = ( xx - param%x ) / radius / param%fact
1358 yy = - param%hemisphere * ( yy - param%y ) / radius / param%fact
1360 dist = sqrt( xx*xx + yy*yy )
1362 lon = param%basepoint_lon + atan2(xx,yy)
1363 lat = param%hemisphere * ( 0.5_dp*pi - 2.0_dp*atan(dist) )
1377 real(
rp),
intent(in) :: lon, lat
1379 real(
rp),
intent(out) :: x, y
1381 real(
dp) :: dlon, latrot, dist
1385 if ( lon==undef .or. lat==undef )
then
1392 dlon = lon - param%basepoint_lon
1394 latrot = 0.5_dp*pi - param%hemisphere * lat
1396 dist = param%fact * radius * tan(0.5_dp*latrot)
1398 xx = param%x + dist * sin(dlon)
1399 yy = param%y - param%hemisphere * dist * cos(dlon)
1401 call xy_unrotate( xx, yy, &
1415 real(
rp),
intent(in) :: lat
1417 real(
rp),
intent(out) :: m1, m2
1420 m1 = param%fact / ( 1.0_dp + sin(param%hemisphere*lat) )
1430 rotc_cos, rotc_sin )
1432 real(
rp),
intent(in) :: lon, lat
1434 real(
rp),
intent(out) :: rotc_cos, rotc_sin
1440 dlon = lon - param%basepoint_lon
1441 if ( dlon > pi ) dlon = dlon - pi*2.0_dp
1442 if ( dlon < -pi ) dlon = dlon + pi*2.0_dp
1444 alpha = - dlon * param%hemisphere &
1447 rotc_cos = cos( alpha )
1448 rotc_sin = sin( alpha )
1465 real(
dp) :: basepoint_lat
1466 real(
dp) :: basepoint_x, basepoint_y
1469 real(
dp) :: latrot, dist
1472 m_lat = info%standard_parallel(1)
1474 basepoint_lat = info%latitude_of_projection_origin
1475 basepoint_x = info%false_easting
1476 basepoint_y = info%false_northing
1481 param%fact = cos(lat0)
1483 if ( param%fact == 0.0_dp )
then
1484 log_error(
"MAPPROJECTION_get_param_Mercator",*)
'M_lat cannot be set to pole point! value=', m_lat
1489 latrot = 0.5_dp*pi - basepoint_lat * d2r
1491 dist = 1.0_dp / tan(0.5_dp*latrot)
1493 param%x = basepoint_x
1494 param%y = basepoint_y - radius * param%fact * log(dist)
1497 log_info(
"MAPPROJECTION_get_param_Mercator",*)
'M_lat = ', m_lat
1498 log_info(
"MAPPROJECTION_get_param_Mercator",*)
'M_fact = ', param%fact
1499 log_info(
"MAPPROJECTION_get_param_Mercator",*)
'eq_x = ', param%x
1500 log_info(
"MAPPROJECTION_get_param_Mercator",*)
'eq_y = ', param%y
1514 real(
rp),
intent(in) :: x, y
1516 real(
rp),
intent(out) :: lon, lat
1521 if ( x==undef .or. y==undef )
then
1528 call xy_rotate( x, y, &
1532 xx = ( xx - param%x ) / radius / param%fact
1533 yy = ( yy - param%y ) / radius / param%fact
1535 lon = xx + param%basepoint_lon
1537 lat = 0.5_dp*pi - 2.0_dp*atan( 1.0_dp/exp(yy) )
1551 real(
rp),
intent(in) :: lon, lat
1553 real(
rp),
intent(out) :: x, y
1555 real(
dp) :: dlon, latrot, dist
1558 if ( lon==undef .or. lat==undef )
then
1565 dlon = lon - param%basepoint_lon
1567 latrot = 0.5_dp*pi - lat
1569 dist = 1.0_dp / tan(0.5_dp*latrot)
1571 x = param%x + radius * param%fact * dlon
1572 y = param%y + radius * param%fact * log(dist)
1584 real(
rp),
intent(in) :: lat
1586 real(
rp),
intent(out) :: m1, m2
1589 m1 = param%fact / cos(real(lat,kind=
dp))
1599 rotc_cos, rotc_sin )
1601 real(
rp),
intent(in) :: lon, lat
1603 real(
rp),
intent(out) :: rotc_cos, rotc_sin
1606 rotc_cos = param%rot_fact_cos
1607 rotc_sin = param%rot_fact_sin
1624 real(
dp) :: basepoint_lat
1625 real(
dp) :: basepoint_x, basepoint_y
1630 ec_lat = info%standard_parallel(1)
1632 basepoint_lat = info%latitude_of_projection_origin
1633 basepoint_x = info%false_easting
1634 basepoint_y = info%false_northing
1639 param%fact = cos(lat0)
1641 if ( param%fact == 0.0_dp )
then
1642 log_error(
"MAPPROJECTION_get_param_EquidistantCylindrical",*)
'EC_lat cannot be set to pole point! value=', ec_lat
1646 param%x = basepoint_x
1647 param%y = basepoint_y - radius * basepoint_lat * d2r
1650 log_info(
"MAPPROJECTION_get_param_EquidistantCylindrical",*)
'EC_lat = ', ec_lat
1651 log_info(
"MAPPROJECTION_get_param_EquidistantCylindrical",*)
'EC_fact = ', param%fact
1652 log_info(
"MAPPROJECTION_get_param_EquidistantCylindrical",*)
'eq_x = ', param%x
1653 log_info(
"MAPPROJECTION_get_param_EquidistantCylindrical",*)
'eq_y = ', param%y
1669 real(
rp),
intent(in) :: x, y
1671 real(
rp),
intent(out) :: lon, lat
1676 if ( x==undef .or. y==undef )
then
1683 call xy_rotate( x, y, &
1687 xx = ( xx - param%x ) / radius / param%fact
1688 yy = ( yy - param%y ) / radius
1690 lon = xx + param%basepoint_lon
1693 if ( abs(lat) > 0.5_dp*pi )
then
1694 log_error(
"MAPPROJECTION_EquidistantCylindrical_xy2lonlat",*)
'Invalid latitude range! value=', lat
1710 real(
rp),
intent(in) :: lon, lat
1712 real(
rp),
intent(out) :: x, y
1718 if ( lon==undef .or. lat==undef )
then
1725 dlon = lon - param%basepoint_lon
1727 xx = param%x + radius * param%fact * dlon
1728 yy = param%y + radius * lat
1730 call xy_unrotate( xx, yy, &
1744 real(
rp),
intent(in) :: lat
1746 real(
rp),
intent(out) :: m1, m2
1748 real(
dp) :: mm1, mm2
1751 mm1 = param%fact / cos(real(lat,kind=
dp))
1754 m1 = sqrt( (param%rot_fact_cos * mm1)**2 + (param%rot_fact_sin * mm2)**2 )
1755 m2 = sqrt( (param%rot_fact_cos * mm2)**2 + (param%rot_fact_sin * mm1)**2 )
1764 rotc_cos, rotc_sin )
1766 real(
rp),
intent(in) :: lon, lat
1768 real(
rp),
intent(out) :: rotc_cos, rotc_sin
1771 rotc_cos = param%rot_fact_cos
1772 rotc_sin = param%rot_fact_sin
1779 subroutine xy_rotate( &
1784 real(
rp),
intent(in) :: x, y
1786 real(
dp),
intent(out) :: xx, yy
1790 xd = x - param%basepoint_x
1791 yd = y - param%basepoint_y
1793 xx = param%basepoint_x + xd * param%rot_fact_cos &
1794 - yd * param%rot_fact_sin
1795 yy = param%basepoint_y + yd * param%rot_fact_cos &
1796 + xd * param%rot_fact_sin
1799 end subroutine xy_rotate
1801 subroutine xy_unrotate( &
1806 real(
dp),
intent(in) :: xx, yy
1808 real(
rp),
intent(out) :: x, y
1810 real(
dp) :: xxd, yyd
1812 xxd = xx - param%basepoint_x
1813 yyd = yy - param%basepoint_y
1815 x = param%basepoint_x + ( xxd * param%rot_fact_cos &
1816 + yyd * param%rot_fact_sin )
1817 y = param%basepoint_y + ( yyd * param%rot_fact_cos &
1818 - xxd * param%rot_fact_sin )
1821 end subroutine xy_unrotate