SCALE-RM
scale_mapprojection.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
10 !-------------------------------------------------------------------------------
11 #include "scalelib.h"
13  !-----------------------------------------------------------------------------
14  !
15  !++ used modules
16  !
17  use scale_precision
18  use scale_io
19  use scale_prof
20  !-----------------------------------------------------------------------------
21  implicit none
22  private
23  !-----------------------------------------------------------------------------
24  !
25  !++ Public procedure
26  !
27  public :: mapprojection_setup
28 
29  public :: mapprojection_get_param
35 
36  public :: mapprojection_xy2lonlat
42 
43  public :: mapprojection_lonlat2xy
49 
50  public :: mapprojection_mapfactor
56 
57  public :: mapprojection_rotcoef
63 
64  interface mapprojection_xy2lonlat
66  module procedure mapprojection_xy2lonlat_2d_initialized
67  module procedure mapprojection_xy2lonlat_0d_param
68  module procedure mapprojection_xy2lonlat_2d_param
69  end interface
70  interface mapprojection_lonlat2xy
72  module procedure mapprojection_lonlat2xy_2d_initialized
73  module procedure mapprojection_lonlat2xy_0d_param
74  module procedure mapprojection_lonlat2xy_2d_param
75  end interface
76 
77  interface mapprojection_mapfactor
79  module procedure mapprojection_mapfactor_param
80  end interface mapprojection_mapfactor
81 
82  interface mapprojection_rotcoef
83  module procedure mapprojection_rotcoef_initialized
84  module procedure mapprojection_rotcoef_param
85  end interface mapprojection_rotcoef
86  !-----------------------------------------------------------------------------
87  !
88  !++ Public parameters & variables
89  !
90  real(rp), public :: mapprojection_basepoint_lon = 135.221_rp ! position of base point in real world [deg]
91  real(rp), public :: mapprojection_basepoint_lat = 34.653_rp ! position of base point in real world [deg]
92 
93  type, public :: mappinginfo
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)
102  real(dp) :: rotation
103  end type mappinginfo
104  type(mappinginfo), public, save :: mapprojection_mappinginfo
105 
106  type, public :: mappingparam
107  real(dp) :: basepoint_x
108  real(dp) :: basepoint_y
109  real(dp) :: basepoint_lon
110  real(dp) :: basepoint_lat
111  real(dp) :: rotation
112  real(dp) :: rot_fact_sin
113  real(dp) :: rot_fact_cos
114  real(dp) :: hemisphere ! hemisphere flag: 1=north, -1=south
115  real(dp) :: x ! position of the pole or equator in the model [m]
116  real(dp) :: y !
117  real(dp) :: fact ! pre-calc factor
118  real(dp) :: c ! conformal factor
119  end type mappingparam
120 
121  !-----------------------------------------------------------------------------
122  !
123  !++ Private procedure
124  !
125  !-----------------------------------------------------------------------------
126  !
127  !++ Private parameters & variables
128  !
129  character(len=H_SHORT), private :: mapprojection_type = 'NONE'
130  ! 'NONE'
131  ! 'LC'
132  ! 'PS'
133  ! 'MER'
134  ! 'EC'
135 
136  real(dp), private :: mapprojection_basepoint_x ! position of base point in the model [m]
137  real(dp), private :: mapprojection_basepoint_y ! position of base point in the model [m]
138  real(dp), private :: mapprojection_basepoint_lonrad
139  real(dp), private :: mapprojection_basepoint_latrad
140  real(dp), private :: mapprojection_rotation = 0.0_dp ! rotation factor
141 
142  real(dp), private :: mapprojection_lc_lat1 = 30.0_dp ! standard latitude1 for L.C. projection [deg]
143  real(dp), private :: mapprojection_lc_lat2 = 60.0_dp ! standard latitude2 for L.C. projection [deg]
144  real(dp), private :: mapprojection_ps_lat ! standard latitude1 for P.S. projection [deg]
145  real(dp), private :: mapprojection_m_lat ! standard latitude1 for Mer. projection [deg]
146  real(dp), private :: mapprojection_ec_lat ! standard latitude1 for E.C. projection [deg]
147 
148  type(mappingparam), private :: mapprojection_mappingparam
149 
150  real(dp), private :: pi
151  real(dp), private :: d2r
152  real(dp), private :: radius
153 
154  interface
155  subroutine xy2lonlat_s( x, y, param, lon, lat )
156  use scale_precision
157  import mappingparam
158  real(rp), intent(in) :: x, y
159  type(mappingparam), intent(in) :: param
160  real(rp), intent(out) :: lon, lat
161  end subroutine xy2lonlat_s
162  subroutine lonlat2xy_s( lon, lat, param, x, y )
163  use scale_precision
164  import mappingparam
165  real(RP), intent(in) :: lon, lat
166  type(mappingparam), intent(in) :: param
167  real(RP), intent(out) :: x, y
168  end subroutine lonlat2xy_s
169  subroutine mapfactor_s( lat, param, m1, m2 )
170  use scale_precision
171  import mappingparam
172  real(RP), intent(in) :: lat
173  type(mappingparam), intent(in) :: param
174  real(RP), intent(out) :: m1, m2
175  end subroutine mapfactor_s
176  subroutine rotcoef_s( lon, lat, param, rotc_cos, rotc_sin )
177  use scale_precision
178  import mappingparam
179  real(RP), intent(in) :: lon, lat
180  type(mappingparam), intent(in) :: param
181  real(RP), intent(out) :: rotc_cos, rotc_sin
182  end subroutine rotcoef_s
183  end interface
184  procedure(xy2lonlat_s), pointer :: xy2lonlat => null()
185  procedure(lonlat2xy_s), pointer :: lonlat2xy => null()
186  procedure(mapfactor_s), pointer :: mapfactor => null()
187  procedure(rotcoef_s), pointer :: rotcoef => null()
188 
189  !-----------------------------------------------------------------------------
190 contains
191  !-----------------------------------------------------------------------------
193  subroutine mapprojection_setup( DOMAIN_CENTER_X, DOMAIN_CENTER_Y )
194  use scale_prc, only: &
195  prc_abort
196  use scale_prc_cartesc, only: &
197  prc_twod
198  use scale_const, only: &
199  undef => const_undef, &
200  pi_rp => const_pi, &
201  d2r_rp => const_d2r, &
202  radius_rp => const_radius
203  implicit none
204 
205  real(rp), intent(in) :: domain_center_x
206  real(rp), intent(in) :: domain_center_y
207 
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, &
219  mapprojection_ec_lat
220 
221  integer :: ierr
222  !---------------------------------------------------------------------------
223 
224  log_newline
225  log_info("MAPPROJECTION_setup",*) 'Setup'
226 
227  pi = real(pi_rp, kind=dp)
228  d2r = real(d2r_rp, kind=dp)
229  radius = real(radius_rp, kind=dp)
230 
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
236 
237 
238  !--- read namelist
239  rewind(io_fid_conf)
240  read(io_fid_conf,nml=param_mapprojection,iostat=ierr)
241 
242  if( ierr < 0 ) then !--- missing
243  log_info("MAPPROJECTION_setup",*) 'Not found namelist. Default used.'
244  elseif( ierr > 0 ) then !--- fatal error
245  log_error("MAPPROJECTION_setup",*) 'Not appropriate names in namelist PARAM_MAPPROJECTION. Check!'
246  call prc_abort
247  endif
248  log_nml(param_mapprojection)
249 
250  if ( prc_twod .and. mapprojection_type .ne. 'NONE' ) then
251  log_error("MAPPROJECTION_setup",*) 'MAPPROJECTION_type must be "NONE" for 2D experiment'
252  call prc_abort
253  end if
254 
255  log_newline
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)
260 
261  mapprojection_mappinginfo%longitude_of_central_meridian = undef
262  mapprojection_mappinginfo%straight_vertical_longitude_from_pole = undef
263  mapprojection_mappinginfo%standard_parallel(:) = undef
264 
265  mapprojection_mappinginfo%longitude_of_projection_origin = mapprojection_basepoint_lon
266  mapprojection_mappinginfo%latitude_of_projection_origin = mapprojection_basepoint_lat
267  mapprojection_mappinginfo%false_easting = mapprojection_basepoint_x
268  mapprojection_mappinginfo%false_northing = mapprojection_basepoint_y
269  mapprojection_mappinginfo%rotation = mapprojection_rotation
270 
271 
272  select case(mapprojection_type)
273  case('NONE')
274  log_info_cont(*) '=> NO map projection'
275  mapprojection_mappinginfo%mapping_name = ""
276 
278  lonlat2xy => mapprojection_lonlat2xy_none
279  mapfactor => mapprojection_mapfactor_none
280  rotcoef => mapprojection_rotcoef_none
281  case('LC')
282  log_info_cont(*) '=> Lambert Conformal projection'
283  mapprojection_mappinginfo%mapping_name = "lambert_conformal_conic"
284  mapprojection_mappinginfo%standard_parallel(:) = (/ mapprojection_lc_lat1, mapprojection_lc_lat2 /)
285  mapprojection_mappinginfo%longitude_of_central_meridian = mapprojection_basepoint_lon
286 
291  case('PS')
292  log_info_cont(*) '=> Polar Stereographic projection'
293  if( mapprojection_ps_lat == undef ) mapprojection_ps_lat = mapprojection_basepoint_lat
294  mapprojection_mappinginfo%mapping_name = "polar_stereographic"
295  mapprojection_mappinginfo%straight_vertical_longitude_from_pole = mapprojection_basepoint_lon
296  mapprojection_mappinginfo%standard_parallel(1) = mapprojection_ps_lat
297 
302  case('MER')
303  log_info_cont(*) '=> Mercator projection'
304  if( mapprojection_m_lat == undef ) mapprojection_m_lat = mapprojection_basepoint_lat
305  mapprojection_mappinginfo%mapping_name = "mercator"
306  mapprojection_mappinginfo%standard_parallel(1) = mapprojection_m_lat
307 
312  case('EC')
313  log_info_cont(*) '=> Equidistant Cylindrical projection'
314  if( mapprojection_ec_lat == undef ) mapprojection_ec_lat = mapprojection_basepoint_lat
315  mapprojection_mappinginfo%mapping_name = "equirectangular"
316  mapprojection_mappinginfo%standard_parallel(1) = mapprojection_ec_lat
317  mapprojection_mappinginfo%longitude_of_central_meridian = mapprojection_basepoint_lon
318 
323  case default
324  log_error("MAPPROJECTION_setup",*) 'Unsupported MAPPROJECTION_type. STOP'
325  call prc_abort
326  endselect
327 
329  mapprojection_mappingparam ) ! (out)
330 
331  return
332  end subroutine mapprojection_setup
333 
334  subroutine mapprojection_get_param( &
335  info, &
336  param )
337  use scale_prc, only: &
338  prc_abort
339  implicit none
340  type(mappinginfo), intent(in) :: info
341  type(mappingparam), intent(out) :: param
342 
343  select case( info%mapping_name )
344  case( "" )
346 
347  param%basepoint_lon = info%longitude_of_projection_origin * d2r
348  case( "lambert_conformal_conic" )
349  call mapprojection_get_param_lambertconformal( info, & ! (in)
350  param ) ! (out)
351  param%basepoint_lon = info%longitude_of_central_meridian * d2r
352  case( "polar_stereographic" )
353  call mapprojection_get_param_polarstereographic( info, & ! (in)
354  param ) ! (out)
355  param%basepoint_lon = info%straight_vertical_longitude_from_pole * d2r
356  case( "mercator" )
357  call mapprojection_get_param_mercator( info, & ! (in)
358  param ) ! (out)
359  param%basepoint_lon = info%longitude_of_projection_origin * d2r
360  case( "equirectangular" )
362  param ) ! (out)
363  param%basepoint_lon = info%longitude_of_central_meridian * d2r
364  case default
365  log_error("MAPPROJECTION_set_param",*) 'Unsupported mapping type. STOP'
366  call prc_abort
367  endselect
368 
369  param%basepoint_x = info%false_easting
370  param%basepoint_y = info%false_northing
371 
372  param%basepoint_lat = info%latitude_of_projection_origin * d2r
373 
374  param%rotation = info%rotation * d2r
375  param%rot_fact_sin = sin(param%rotation)
376  param%rot_fact_cos = cos(param%rotation)
377 
378  return
379  end subroutine mapprojection_get_param
380 
381  !-----------------------------------------------------------------------------
384  x, y, &
385  lon, lat )
386  implicit none
387  real(RP), intent(in) :: x, y
388  real(RP), intent(out) :: lon, lat ! [rad]
389  !---------------------------------------------------------------------------
390 
391  call xy2lonlat( x, y, & ! (in)
392  mapprojection_mappingparam, & ! (in)
393  lon, lat ) ! (out)
394 
395  return
397 
398  subroutine mapprojection_xy2lonlat_2d_initialized( &
399  IA, IS, IE, JA, JS, JE, &
400  x, y, &
401  lon, lat )
402  implicit none
403  integer, intent(in) :: IA, IS, IE
404  integer, intent(in) :: JA, JS, JE
405 
406  real(RP), intent(in) :: x(IA,JA)
407  real(RP), intent(in) :: y(IA,JA)
408  real(RP), intent(out) :: lon(IA,JA) ! [rad]
409  real(RP), intent(out) :: lat(IA,JA) ! [rad]
410 
411  integer :: i, j
412  !---------------------------------------------------------------------------
413 
414  !$omp parallel do
415  do j = js, je
416  do i = is, ie
417  call mapprojection_xy2lonlat_0d_initialized( x(i,j), y(i,j), lon(i,j), lat(i,j) )
418  end do
419  end do
420 
421  return
422  end subroutine mapprojection_xy2lonlat_2d_initialized
423 
424  subroutine mapprojection_xy2lonlat_0d_param( &
425  x, y, &
426  mapping_name, &
427  mapping_param, &
428  lon, lat )
429  use scale_prc, only: &
430  prc_abort
431  implicit none
432  real(RP), intent(in) :: x, y
433  character(len=*), intent(in) :: mapping_name
434  type(mappingparam), intent(in) :: mapping_param
435 
436  real(RP), intent(out) :: lon, lat ! [rad]
437  !---------------------------------------------------------------------------
438 
439  select case( mapping_name )
440  case( "" )
441  call mapprojection_xy2lonlat_none( x, y, & ! (in)
442  mapping_param, & ! (in)
443  lon, lat ) ! (out)
444  case( "lambert_conformal_conic" )
445  call mapprojection_xy2lonlat_lambertconformal( x, y, & ! (in)
446  mapping_param, & ! (in)
447  lon, lat ) ! (out)
448  case( "polar_stereographic" )
449  call mapprojection_xy2lonlat_polarstereographic( x, y, & ! (in)
450  mapping_param, & ! (in)
451  lon, lat ) ! (out)
452  case( "mercator" )
453  call mapprojection_xy2lonlat_mercator( x, y, & ! (in)
454  mapping_param, & ! (in)
455  lon, lat ) ! (out)
456  case( "equirectangular" )
458  mapping_param, & ! (in)
459  lon, lat ) ! (out)
460  case default
461  log_error("MAPPROJECTION_xy2lonlat_0D_param",*) 'Unsupported mapping type. STOP: ', trim(mapping_name)
462  call prc_abort
463  endselect
464 
465  return
466  end subroutine mapprojection_xy2lonlat_0d_param
467 
469  IA, IS, IE, JA, JS, JE, &
470  x, y, &
471  mapping_name, &
472  mapping_param, &
473  lon, lat )
474  use scale_prc, only: &
475  prc_abort
476  implicit none
477  integer, intent(in) :: IA, IS, IE
478  integer, intent(in) :: JA, JS, JE
479 
480  real(RP), intent(in) :: x(IA,JA)
481  real(RP), intent(in) :: y(IA,JA)
482  character(len=*), intent(in) :: mapping_name
483  type(mappingparam), intent(in) :: mapping_param
484 
485  real(RP), intent(out) :: lon(IA,JA) ! [rad]
486  real(RP), intent(out) :: lat(IA,JA) ! [rad]
487 
488  integer :: i, j
489  !---------------------------------------------------------------------------
490 
491  select case( mapping_name )
492  case( "" )
493  !$omp parallel do
494  do j = js, je
495  do i = is, ie
496  call mapprojection_xy2lonlat_none( x(i,j), y(i,j), & ! (in)
497  mapping_param, & ! (in)
498  lon(i,j), lat(i,j) ) ! (out)
499  end do
500  end do
501  case( "lambert_conformal_conic" )
502  !$omp parallel do
503  do j = js, je
504  do i = is, ie
505  call mapprojection_xy2lonlat_lambertconformal( x(i,j), y(i,j), & ! (in)
506  mapping_param, & ! (in)
507  lon(i,j), lat(i,j) ) ! (out)
508  end do
509  end do
510  case( "polar_stereographic" )
511  !$omp parallel do
512  do j = js, je
513  do i = is, ie
514  call mapprojection_xy2lonlat_polarstereographic( x(i,j), y(i,j), & ! (in)
515  mapping_param, & ! (in)
516  lon(i,j), lat(i,j) ) ! (out)
517  end do
518  end do
519  case( "mercator" )
520  !$omp parallel do
521  do j = js, je
522  do i = is, ie
523  call mapprojection_xy2lonlat_mercator( x(i,j), y(i,j), & ! (in)
524  mapping_param, & ! (in)
525  lon(i,j), lat(i,j) ) ! (out)
526  end do
527  end do
528  case( "equirectangular" )
529  !$omp parallel do
530  do j = js, je
531  do i = is, ie
532  call mapprojection_xy2lonlat_equidistantcylindrical( x(i,j), y(i,j), & ! (in)
533  mapping_param, & ! (in)
534  lon(i,j), lat(i,j) ) ! (out)
535  end do
536  end do
537  case default
538  log_error("MAPPROJECTION_xy2lonlat_2D_param",*) 'Unsupported mapping type. STOP: ', trim(mapping_name)
539  call prc_abort
540  endselect
541 
542  return
543  end subroutine mapprojection_xy2lonlat_2d_param
544 
545  !-----------------------------------------------------------------------------
548  lon, lat, &
549  x, y )
550  implicit none
551  real(RP), intent(in) :: lon, lat ! [rad]
552  real(RP), intent(out) :: x, y
553  !---------------------------------------------------------------------------
554 
555  call lonlat2xy( lon, lat, & ! (in)
556  mapprojection_mappingparam, & ! (in)
557  x, y ) ! (out)
558 
559  return
561 
562  subroutine mapprojection_lonlat2xy_2d_initialized( &
563  IA, IS, IE, JA, JS, JE, &
564  lon, lat, &
565  x, y )
566  implicit none
567  integer, intent(in) :: IA, IS, IE
568  integer, intent(in) :: JA, JS, JE
569 
570  real(RP), intent(in) :: lon(IA,JA) ! [rad]
571  real(RP), intent(in) :: lat(IA,JA) ! [rad]
572  real(RP), intent(out) :: x(IA,JA)
573  real(RP), intent(out) :: y(IA,JA)
574 
575  integer :: i, j
576  !---------------------------------------------------------------------------
577 
578  !$omp parallel do
579  do j = js, je
580  do i = is, ie
581  call mapprojection_lonlat2xy_0d_initialized( lon(i,j), lat(i,j), x(i,j), y(i,j) )
582  end do
583  end do
584 
585  return
586  end subroutine mapprojection_lonlat2xy_2d_initialized
587 
588  subroutine mapprojection_lonlat2xy_0d_param( &
589  lon, lat, &
590  mapping_name, &
591  mapping_param, &
592  x, y )
593  use scale_prc, only: &
594  prc_abort
595  implicit none
596  real(RP), intent(in) :: lon, lat ! [rad]
597  character(len=*), intent(in) :: mapping_name
598  type(mappingparam), intent(in) :: mapping_param
599 
600  real(RP), intent(out) :: x, y
601  !---------------------------------------------------------------------------
602 
603  select case( mapping_name )
604  case( "" )
605  call mapprojection_lonlat2xy_none( lon, lat, & ! (in)
606  mapping_param, & ! (in)
607  x, y ) ! (out)
608  case( "lambert_conformal_conic" )
609  call mapprojection_lonlat2xy_lambertconformal( lon, lat, & ! (in)
610  mapping_param, & ! (in)
611  x, y ) ! (out)
612  case( "polar_stereographic" )
613  call mapprojection_lonlat2xy_polarstereographic( lon, lat, & ! (in)
614  mapping_param, & ! (in)
615  x, y ) ! (out)
616  case( "mercator" )
617  call mapprojection_lonlat2xy_mercator( lon, lat, & ! (in)
618  mapping_param, & ! (in)
619  x, y ) ! (out)
620  case( "equirectangular" )
621  call mapprojection_lonlat2xy_equidistantcylindrical( lon, lat, & ! (in)
622  mapping_param, & ! (in)
623  x, y ) ! (out)
624  case default
625  log_error("MAPPROJECTION_lonlat2xy_0D_param",*) 'Unsupported mapping type. STOP: ', trim(mapping_name)
626  call prc_abort
627  endselect
628 
629  return
630  end subroutine mapprojection_lonlat2xy_0d_param
631 
633  IA, IS, IE, JA, JS, JE, &
634  lon, lat, &
635  mapping_name, &
636  mapping_param, &
637  x, y )
638  use scale_prc, only: &
639  prc_abort
640  implicit none
641  integer, intent(in) :: IA, IS, IE
642  integer, intent(in) :: JA, JS, JE
643 
644  real(RP), intent(in) :: lon(IA,JA) ! [rad]
645  real(RP), intent(in) :: lat(IA,JA) ! [rad]
646  character(len=*), intent(in) :: mapping_name
647  type(mappingparam), intent(in) :: mapping_param
648 
649  real(RP), intent(out) :: x(IA,JA)
650  real(RP), intent(out) :: y(IA,JA)
651 
652  integer :: i, j
653  !---------------------------------------------------------------------------
654 
655  select case( mapping_name )
656  case( "" )
657  !$omp parallel do
658  do j = js, je
659  do i = is, ie
660  call mapprojection_lonlat2xy_none( lon(i,j), lat(i,j), & ! (in)
661  mapping_param, & ! (in)
662  x(i,j), y(i,j) ) ! (out)
663  end do
664  end do
665  case( "lambert_conformal_conic" )
666  !$omp parallel do
667  do j = js, je
668  do i = is, ie
669  call mapprojection_lonlat2xy_lambertconformal( lon(i,j), lat(i,j), & ! (in)
670  mapping_param, & ! (in)
671  x(i,j), y(i,j) ) ! (out)
672  end do
673  end do
674  case( "polar_stereographic" )
675  !$omp parallel do
676  do j = js, je
677  do i = is, ie
678  call mapprojection_lonlat2xy_polarstereographic( lon(i,j), lat(i,j), & ! (in)
679  mapping_param, & ! (in)
680  x(i,j), y(i,j) ) ! (out)
681  end do
682  end do
683  case( "mercator" )
684  !$omp parallel do
685  do j = js, je
686  do i = is, ie
687  call mapprojection_lonlat2xy_mercator( lon(i,j), lat(i,j), & ! (in)
688  mapping_param, & ! (in)
689  x(i,j), y(i,j) ) ! (out)
690  end do
691  end do
692  case( "equirectangular" )
693  !$omp parallel do
694  do j = js, je
695  do i = is, ie
696  call mapprojection_lonlat2xy_equidistantcylindrical( lon(i,j), lat(i,j), & ! (in)
697  mapping_param, & ! (in)
698  x(i,j), y(i,j) ) ! (out)
699  end do
700  end do
701  case default
702  log_error("MAPPROJECTION_lonlat2xy_2D_param",*) 'Unsupported mapping type. STOP: ', trim(mapping_name)
703  call prc_abort
704  endselect
705 
706  return
707  end subroutine mapprojection_lonlat2xy_2d_param
708 
709  !-----------------------------------------------------------------------------
712  IA, IS, IE, JA, JS, JE, &
713  lat, &
714  m1, m2 )
715  implicit none
716  integer, intent(in) :: IA, IS, IE
717  integer, intent(in) :: JA, JS, JE
718 
719  real(RP), intent(in) :: lat(IA,JA) ! [rad]
720  real(RP), intent(out) :: m1 (IA,JA)
721  real(RP), intent(out) :: m2 (IA,JA)
722 
723  integer :: i, j
724  !---------------------------------------------------------------------------
725 
726  !$omp parallel do
727  do j = js, je
728  do i = is, ie
729  call mapfactor( lat(i,j), & ! (in)
730  mapprojection_mappingparam, & ! (in)
731  m1(i,j), m2(i,j) ) ! (out)
732  end do
733  end do
734 
735  return
737 
738  subroutine mapprojection_mapfactor_param( &
739  IA, IS, IE, JA, JS, JE, &
740  lat, &
741  mapping_name, &
742  mapping_param, &
743  m1, m2 )
744  use scale_prc, only: &
745  prc_abort
746  implicit none
747  integer, intent(in) :: IA, IS, IE
748  integer, intent(in) :: JA, JS, JE
749 
750  real(RP), intent(in) :: lat(IA,JA) ! [rad]
751  character(len=*), intent(in) :: mapping_name
752  type(mappingparam), intent(in) :: mapping_param
753 
754  real(RP), intent(out) :: m1 (IA,JA)
755  real(RP), intent(out) :: m2 (IA,JA)
756 
757  integer :: i, j
758  !---------------------------------------------------------------------------
759 
760  select case( mapping_name )
761  case( "" )
762  !$omp parallel do
763  do j = js, je
764  do i = is, ie
765  call mapprojection_mapfactor_none( lat(i,j), & ! (in)
766  mapping_param, & ! (in)
767  m1(i,j), m2(i,j) ) ! (out)
768  end do
769  end do
770  case( "lambert_conformal_conic" )
771  !$omp parallel do
772  do j = js, je
773  do i = is, ie
774  call mapprojection_mapfactor_lambertconformal( lat(i,j), & ! (in)
775  mapping_param, & ! (in)
776  m1(i,j), m2(i,j) ) ! (out)
777  end do
778  end do
779  case( "polar_stereographic" )
780  !$omp parallel do
781  do j = js, je
782  do i = is, ie
783  call mapprojection_mapfactor_polarstereographic( lat(i,j), & ! (in)
784  mapping_param, & ! (in)
785  m1(i,j), m2(i,j) ) ! (out)
786  end do
787  end do
788  case( "mercator" )
789  !$omp parallel do
790  do j = js, je
791  do i = is, ie
792  call mapprojection_mapfactor_mercator( lat(i,j), & ! (in)
793  mapping_param, & ! (in)
794  m1(i,j), m2(i,j) ) ! (out)
795  end do
796  end do
797  case( "equirectangular" )
798  !$omp parallel do
799  do j = js, je
800  do i = is, ie
801  call mapprojection_mapfactor_equidistantcylindrical( lat(i,j), & ! (in)
802  mapping_param, & ! (in)
803  m1(i,j), m2(i,j) ) ! (out)
804  end do
805  end do
806  case default
807  log_error("MAPPROJECTION_mapfactor_param",*) 'Unsupported mapping type. STOP: ', trim(mapping_name)
808  call prc_abort
809  endselect
810 
811  return
812  end subroutine mapprojection_mapfactor_param
813 
814  !-----------------------------------------------------------------------------
818  IA, IS, IE, JA, JS, JE, &
819  lon, lat, &
820  rotc_cos, rotc_sin )
821  implicit none
822  integer, intent(in) :: IA, IS, IE
823  integer, intent(in) :: JA, JS, JE
824 
825  real(RP), intent(in) :: lon(IA,JA) ! [rad]
826  real(RP), intent(in) :: lat(IA,JA) ! [rad]
827  real(RP), intent(out) :: rotc_cos(IA,JA)
828  real(RP), intent(out) :: rotc_sin(IA,JA)
829 
830  integer :: i, j
831  !---------------------------------------------------------------------------
832 
833  !$omp parallel do
834  do j = js, je
835  do i = is, ie
836  call rotcoef( lon(i,j), lat(i,j), & ! (in)
837  mapprojection_mappingparam, & ! (in)
838  rotc_cos(i,j), rotc_sin(i,j) ) ! (out)
839  end do
840  end do
841 
842  return
843  end subroutine mapprojection_rotcoef_initialized
844 
845  subroutine mapprojection_rotcoef_param( &
846  IA, IS, IE, JA, JS, JE, &
847  lon, lat, &
848  mapping_name, &
849  mapping_param, &
850  rotc_cos, rotc_sin )
851  use scale_prc, only: &
852  prc_abort
853  implicit none
854  integer, intent(in) :: IA, IS, IE
855  integer, intent(in) :: JA, JS, JE
856 
857  real(RP), intent(in) :: lon(IA,JA) ! [rad]
858  real(RP), intent(in) :: lat(IA,JA) ! [rad]
859  character(len=*), intent(in) :: mapping_name
860  type(mappingparam), intent(in) :: mapping_param
861 
862  real(RP), intent(out) :: rotc_cos(IA,JA)
863  real(RP), intent(out) :: rotc_sin(IA,JA)
864 
865  integer :: i, j
866  !---------------------------------------------------------------------------
867 
868  select case( mapping_name )
869  case( "" )
870  !$omp parallel do
871  do j = js, je
872  do i = is, ie
873  call mapprojection_rotcoef_none( lon(i,j), lat(i,j), & ! (in)
874  mapping_param, & ! (in)
875  rotc_cos(i,j), rotc_sin(i,j) ) ! (out)
876  end do
877  end do
878  case( "lambert_conformal_conic" )
879  !$omp parallel do
880  do j = js, je
881  do i = is, ie
882  call mapprojection_rotcoef_lambertconformal( lon(i,j), lat(i,j), & ! (in)
883  mapping_param, & ! (in)
884  rotc_cos(i,j), rotc_sin(i,j) ) ! (out)
885  end do
886  end do
887  case( "polar_stereographic" )
888  !$omp parallel do
889  do j = js, je
890  do i = is, ie
891  call mapprojection_rotcoef_polarstereographic( lon(i,j), lat(i,j), & ! (in)
892  mapping_param, & ! (in)
893  rotc_cos(i,j), rotc_sin(i,j) ) ! (out)
894  end do
895  end do
896  case( "mercator" )
897  !$omp parallel do
898  do j = js, je
899  do i = is, ie
900  call mapprojection_rotcoef_mercator( lon(i,j), lat(i,j), & ! (in)
901  mapping_param, & ! (in)
902  rotc_cos(i,j), rotc_sin(i,j) ) ! (out)
903  end do
904  end do
905  case( "equirectangular" )
906  !$omp parallel do
907  do j = js, je
908  do i = is, ie
909  call mapprojection_rotcoef_equidistantcylindrical( lon(i,j), lat(i,j), & ! (in)
910  mapping_param, & ! (in)
911  rotc_cos(i,j), rotc_sin(i,j) ) ! (out)
912  end do
913  end do
914  case default
915  log_error("MAPPROJECTION_rotcoef_param",*) 'Unsupported mapping type. STOP: ', trim(mapping_name)
916  call prc_abort
917  endselect
918 
919  return
920  end subroutine mapprojection_rotcoef_param
921 
922  !-----------------------------------------------------------------------------
925  implicit none
926  !---------------------------------------------------------------------------
927 
928  return
929  end subroutine mapprojection_get_param_none
930 
931  !-----------------------------------------------------------------------------
933  subroutine mapprojection_xy2lonlat_none( &
934  x, y, &
935  param, &
936  lon, lat )
937  use scale_const, only: &
938  undef => const_undef
939  implicit none
940  real(rp), intent(in) :: x, y
941  type(mappingparam), intent(in) :: param
942  real(rp), intent(out) :: lon, lat ! [rad]
943 
944  real(dp) :: xx, yy
945  real(dp) :: gno1, gno2
946  real(dp) :: rho, gmm
947  real(dp) :: gmmc, gmms
948  real(dp) :: latc, lats
949  !---------------------------------------------------------------------------
950 
951  if ( x==undef .or. y==undef ) then
952  lon = undef
953  lat = undef
954 
955  return
956  end if
957 
958  call xy_rotate( x, y, & ! (in)
959  param, & ! (in)
960  xx, yy ) ! (out)
961 
962  gno1 = ( x - param%basepoint_x ) / radius
963  gno2 = ( y - param%basepoint_y ) / radius
964 
965  rho = sqrt( gno1 * gno1 + gno2 * gno2 )
966  gmm = atan( rho )
967 
968  if ( rho == 0.0_dp ) then
969  lon = param%basepoint_lon
970  lat = param%basepoint_lat
971  else
972  gmmc = cos(gmm)
973  gmms = sin(gmm)
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 )
981  endif
982 
983  return
984  end subroutine mapprojection_xy2lonlat_none
985 
986  subroutine mapprojection_lonlat2xy_none( &
987  lon, lat, &
988  param, &
989  x, y )
990  use scale_const, only: &
991  undef => const_undef
992  implicit none
993  real(rp), intent(in) :: lon, lat ! [rad]
994  type(mappingparam), intent(in) :: param
995  real(rp), intent(out) :: x, y
996 
997  real(dp) :: xx, yy
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
1001  real(dp) :: cos_gmm
1002  !---------------------------------------------------------------------------
1003  ! http://mathworld.wolfram.com/GnomonicProjection.html
1004 
1005  if ( lon==undef .or. lat==undef ) then
1006  x = undef
1007  y = undef
1008 
1009  return
1010  end if
1011 
1012  lat_d = real(lat,kind=dp)
1013  lat0_d = param%basepoint_lat
1014 
1015  dlon = lon - param%basepoint_lon
1016 
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)
1021  dlon_c = cos(dlon)
1022  dlon_s = sin(dlon)
1023 
1024  cos_gmm = lat0_d_s * lat_d_s &
1025  + lat0_d_c * lat_d_c * dlon_c
1026 
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
1030 
1031  xx = gno1 * radius + param%basepoint_x
1032  yy = gno2 * radius + param%basepoint_y
1033 
1034  call xy_unrotate( xx, yy, & ! (in)
1035  param, & ! (in)
1036  x, y ) ! (out)
1037 
1038  return
1039  end subroutine mapprojection_lonlat2xy_none
1040 
1041  !-----------------------------------------------------------------------------
1043  subroutine mapprojection_mapfactor_none( &
1044  lat, &
1045  param, &
1046  m1, m2 )
1047  use scale_const, only: &
1048  undef => const_undef
1049  implicit none
1050  real(rp), intent(in) :: lat ! [rad]
1051  type(mappingparam), intent(in) :: param
1052  real(rp), intent(out) :: m1, m2
1053  !---------------------------------------------------------------------------
1054 
1055  m1 = 1.0_rp
1056  m2 = m1
1057 
1058  return
1059  end subroutine mapprojection_mapfactor_none
1060 
1061  !-----------------------------------------------------------------------------
1063  subroutine mapprojection_rotcoef_none( &
1064  lon, lat, &
1065  param, &
1066  rotc_cos, rotc_sin )
1067  implicit none
1068  real(rp), intent(in) :: lon, lat
1069  type(mappingparam), intent(in) :: param
1070  real(rp), intent(out) :: rotc_cos, rotc_sin
1071  !---------------------------------------------------------------------------
1072 
1073  rotc_cos = param%rot_fact_cos
1074  rotc_sin = param%rot_fact_sin
1075 
1076  return
1077  end subroutine mapprojection_rotcoef_none
1078 
1079  !-----------------------------------------------------------------------------
1082  info, &
1083  param )
1084  use scale_prc, only: &
1085  prc_abort
1086  implicit none
1087  type(mappinginfo), intent(in) :: info
1088  type(mappingparam), intent(out) :: param
1089 
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
1094 
1095  real(dp) :: lat1rot, lat2rot
1096  real(dp) :: dlon, latrot, dist
1097  !---------------------------------------------------------------------------
1098 
1099  lc_lat1 = info%standard_parallel(1)
1100  lc_lat2 = info%standard_parallel(2)
1101 
1102  basepoint_lat = info%latitude_of_projection_origin
1103  basepoint_x = info%false_easting
1104  basepoint_y = info%false_northing
1105 
1106  if ( lc_lat1 > lc_lat2 ) then
1107  log_error("MAPPROJECTION_get_param_LambertConformal",*) 'Please set LC_lat1 <= LC_lat2 in degree. STOP'
1108  call prc_abort
1109  endif
1110 
1111  ! check hemisphere: 1=north, -1=south
1112  param%hemisphere = sign(1.0_dp,lc_lat1+lc_lat2)
1113 
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
1117  ! calc conformal factor c
1118  param%c = ( log( sin(lat1rot) ) - log( sin(lat2rot) ) ) &
1119  / ( log( tan(0.5_dp*lat1rot) ) - log( tan(0.5_dp*lat2rot) ) )
1120  else
1121  param%c = cos(lat1rot)
1122  end if
1123  ! pre-calc factor
1124  param%fact = sin(lat1rot) / param%c / tan(0.5_dp*lat1rot)**param%c
1125 
1126  ! calc (x,y) at pole point
1127  dlon = 0.0_dp
1128 
1129  latrot = 0.5_dp*pi - param%hemisphere * basepoint_lat * d2r
1130 
1131  dist = param%fact * radius * tan(0.5_dp*latrot)**param%c
1132 
1133  param%x = basepoint_x - dist * sin(param%c*dlon)
1134  param%y = basepoint_y + param%hemisphere * dist * cos(param%c*dlon)
1135 
1136  log_newline
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
1145 
1146  return
1148 
1149  !-----------------------------------------------------------------------------
1152  x, y, &
1153  param, &
1154  lon, lat )
1155  use scale_const, only: &
1156  undef => const_undef
1157  implicit none
1158  real(rp), intent(in) :: x, y
1159  type(mappingparam), intent(in) :: param
1160  real(rp), intent(out) :: lon, lat ! [rad]
1161 
1162  real(dp) :: xx, yy, dist
1163  !---------------------------------------------------------------------------
1164 
1165  if ( x==undef .or. y==undef ) then
1166  lon = undef
1167  lat = undef
1168 
1169  return
1170  end if
1171 
1172  call xy_rotate( x, y, & ! (in)
1173  param, & ! (in)
1174  xx, yy ) ! (out)
1175  xx = ( xx - param%x ) / radius / param%fact
1176  yy = - param%hemisphere * ( yy - param%y ) / radius / param%fact
1177 
1178  dist = sqrt( xx*xx + yy*yy )
1179 
1180  lon = param%basepoint_lon + atan2(xx,yy) / param%c
1181 
1182  ! check hemisphere: 1=north, -1=south
1183  lat = param%hemisphere * ( 0.5_dp*pi - 2.0_dp*atan( dist**(1.0_dp/param%c) ) )
1184 
1185  return
1187 
1188  !-----------------------------------------------------------------------------
1191  lon, lat, &
1192  param, &
1193  x, y )
1194  use scale_const, only: &
1195  undef => const_undef
1196  implicit none
1197  real(rp), intent(in) :: lon, lat ! [rad]
1198  type(mappingparam), intent(in) :: param
1199  real(rp), intent(out) :: x, y
1200 
1201  real(dp) :: xx, yy
1202  real(dp) :: dlon, latrot, dist
1203  !---------------------------------------------------------------------------
1204 
1205  if ( lon==undef .or. lat==undef ) then
1206  x = undef
1207  y = undef
1208 
1209  return
1210  end if
1211 
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
1215 
1216  latrot = 0.5_dp*pi - param%hemisphere * lat
1217 
1218  dist = param%fact * radius * tan(0.5_dp*latrot)**param%c
1219 
1220  xx = param%x + dist * sin(param%c*dlon)
1221  yy = param%y - param%hemisphere * dist * cos(param%c*dlon)
1222 
1223  call xy_unrotate( xx, yy, & ! (in)
1224  param, & ! (in)
1225  x, y ) ! (out)
1226 
1227  return
1229 
1230  !-----------------------------------------------------------------------------
1233  lat, &
1234  param, &
1235  m1, m2 )
1236  implicit none
1237  real(rp), intent(in) :: lat ! [rad]
1238  type(mappingparam), intent(in) :: param
1239  real(rp), intent(out) :: m1, m2
1240 
1241  real(dp) :: latrot
1242  !---------------------------------------------------------------------------
1243 
1244  latrot = 0.5_dp*pi - param%hemisphere * lat
1245 
1246  m1 = param%fact / sin(latrot) * param%c * tan(0.5_dp*latrot)**param%c
1247  m2 = m1
1248 
1249  return
1251 
1252  !-----------------------------------------------------------------------------
1254  lon, lat, &
1255  param, &
1256  rotc_cos, rotc_sin )
1257  implicit none
1258  real(rp), intent(in) :: lon, lat
1259  type(mappingparam), intent(in) :: param
1260  real(rp), intent(out) :: rotc_cos, rotc_sin
1261 
1262  real(dp) :: dlon
1263  real(dp) :: alpha
1264  !---------------------------------------------------------------------------
1265 
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
1269 
1270  alpha = - param%c * dlon * param%hemisphere &
1271  + param%rotation
1272 
1273  rotc_cos = cos( alpha )
1274  rotc_sin = sin( alpha )
1275 
1276  return
1278 
1279  !-----------------------------------------------------------------------------
1282  info, &
1283  param )
1284  implicit none
1285  type(mappinginfo), intent(in) :: info
1286  type(mappingparam), intent(out) :: param
1287 
1288  real(dp) :: ps_lat
1289  real(dp) :: basepoint_lat
1290  real(dp) :: basepoint_x, basepoint_y
1291 
1292  real(dp) :: lat0
1293  real(dp) :: dlon, latrot, dist
1294  !---------------------------------------------------------------------------
1295 
1296  ps_lat = info%standard_parallel(1)
1297 
1298  basepoint_lat = info%latitude_of_projection_origin
1299  basepoint_x = info%false_easting
1300  basepoint_y = info%false_northing
1301 
1302  ! check hemisphere: 1=north, -1=south
1303  param%hemisphere = sign(1.0_dp,ps_lat)
1304 
1305  lat0 = param%hemisphere * ps_lat * d2r
1306 
1307  ! pre-calc factor
1308  param%fact = 1.0_dp + sin(lat0)
1309 
1310  ! calc (x,y) at pole point
1311  dlon = 0.0_dp
1312 
1313  latrot = 0.5_dp*pi - param%hemisphere * basepoint_lat * d2r
1314 
1315  dist = param%fact * radius * tan(0.5_dp*latrot)
1316 
1317  param%x = basepoint_x - dist * sin(dlon)
1318  param%y = basepoint_y + param%hemisphere * dist * cos(dlon)
1319 
1320  log_newline
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
1326 
1327  return
1329 
1330  !-----------------------------------------------------------------------------
1333  x, y, &
1334  param, &
1335  lon, lat )
1336  use scale_const, only: &
1337  undef => const_undef
1338  implicit none
1339  real(rp), intent(in) :: x, y
1340  type(mappingparam), intent(in) :: param
1341  real(rp), intent(out) :: lon, lat ! [rad]
1342 
1343  real(dp) :: xx, yy, dist
1344  !---------------------------------------------------------------------------
1345 
1346  if ( x==undef .or. y==undef ) then
1347  lon = undef
1348  lat = undef
1349 
1350  return
1351  end if
1352 
1353  call xy_rotate( x, y, & ! (in)
1354  param, & ! (in)
1355  xx, yy ) ! (out)
1356 
1357  xx = ( xx - param%x ) / radius / param%fact
1358  yy = - param%hemisphere * ( yy - param%y ) / radius / param%fact
1359 
1360  dist = sqrt( xx*xx + yy*yy )
1361 
1362  lon = param%basepoint_lon + atan2(xx,yy)
1363  lat = param%hemisphere * ( 0.5_dp*pi - 2.0_dp*atan(dist) )
1364 
1365  return
1367 
1368  !-----------------------------------------------------------------------------
1371  lon, lat, &
1372  param, &
1373  x, y )
1374  use scale_const, only: &
1375  undef => const_undef
1376  implicit none
1377  real(rp), intent(in) :: lon, lat ! [rad]
1378  type(mappingparam), intent(in) :: param
1379  real(rp), intent(out) :: x, y
1380 
1381  real(dp) :: dlon, latrot, dist
1382  real(dp) :: xx, yy
1383  !---------------------------------------------------------------------------
1384 
1385  if ( lon==undef .or. lat==undef ) then
1386  x = undef
1387  y = undef
1388 
1389  return
1390  end if
1391 
1392  dlon = lon - param%basepoint_lon
1393 
1394  latrot = 0.5_dp*pi - param%hemisphere * lat
1395 
1396  dist = param%fact * radius * tan(0.5_dp*latrot)
1397 
1398  xx = param%x + dist * sin(dlon)
1399  yy = param%y - param%hemisphere * dist * cos(dlon)
1400 
1401  call xy_unrotate( xx, yy, & ! (in)
1402  param, & ! (in)
1403  x, y ) ! (out)
1404 
1405  return
1407 
1408  !-----------------------------------------------------------------------------
1411  lat, &
1412  param, &
1413  m1, m2 )
1414  implicit none
1415  real(rp), intent(in) :: lat ! [rad]
1416  type(mappingparam), intent(in) :: param
1417  real(rp), intent(out) :: m1, m2
1418  !---------------------------------------------------------------------------
1419 
1420  m1 = param%fact / ( 1.0_dp + sin(param%hemisphere*lat) )
1421  m2 = m1
1422 
1423  return
1425 
1426  !-----------------------------------------------------------------------------
1428  lon, lat, &
1429  param, &
1430  rotc_cos, rotc_sin )
1431  implicit none
1432  real(rp), intent(in) :: lon, lat
1433  type(mappingparam), intent(in) :: param
1434  real(rp), intent(out) :: rotc_cos, rotc_sin
1435 
1436  real(dp) :: dlon
1437  real(dp) :: alpha
1438  !---------------------------------------------------------------------------
1439 
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
1443 
1444  alpha = - dlon * param%hemisphere &
1445  + param%rotation
1446 
1447  rotc_cos = cos( alpha )
1448  rotc_sin = sin( alpha )
1449 
1450  return
1452 
1453  !-----------------------------------------------------------------------------
1455  subroutine mapprojection_get_param_mercator( &
1456  info, &
1457  param )
1458  use scale_prc, only: &
1459  prc_abort
1460  implicit none
1461  type(mappinginfo), intent(in) :: info
1462  type(mappingparam), intent(out) :: param
1463 
1464  real(dp) :: m_lat
1465  real(dp) :: basepoint_lat
1466  real(dp) :: basepoint_x, basepoint_y
1467 
1468  real(dp) :: lat0
1469  real(dp) :: latrot, dist
1470  !---------------------------------------------------------------------------
1471 
1472  m_lat = info%standard_parallel(1)
1473 
1474  basepoint_lat = info%latitude_of_projection_origin
1475  basepoint_x = info%false_easting
1476  basepoint_y = info%false_northing
1477 
1478  lat0 = m_lat * d2r
1479 
1480  ! pre-calc factor
1481  param%fact = cos(lat0)
1482 
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
1485  call prc_abort
1486  endif
1487 
1488  ! calc (x,y) at (lon,lat) = (base,0)
1489  latrot = 0.5_dp*pi - basepoint_lat * d2r
1490 
1491  dist = 1.0_dp / tan(0.5_dp*latrot)
1492 
1493  param%x = basepoint_x
1494  param%y = basepoint_y - radius * param%fact * log(dist)
1495 
1496  log_newline
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
1501 
1502  return
1503  end subroutine mapprojection_get_param_mercator
1504 
1505  !-----------------------------------------------------------------------------
1507  subroutine mapprojection_xy2lonlat_mercator( &
1508  x, y, &
1509  param, &
1510  lon, lat )
1511  use scale_const, only: &
1512  undef => const_undef
1513  implicit none
1514  real(rp), intent(in) :: x, y
1515  type(mappingparam), intent(in) :: param
1516  real(rp), intent(out) :: lon, lat ! [rad]
1517 
1518  real(dp) :: xx, yy
1519  !---------------------------------------------------------------------------
1520 
1521  if ( x==undef .or. y==undef ) then
1522  lon = undef
1523  lat = undef
1524 
1525  return
1526  end if
1527 
1528  call xy_rotate( x, y, & ! (in)
1529  param, & ! (in)
1530  xx, yy ) ! (out)
1531 
1532  xx = ( xx - param%x ) / radius / param%fact
1533  yy = ( yy - param%y ) / radius / param%fact
1534 
1535  lon = xx + param%basepoint_lon
1536 
1537  lat = 0.5_dp*pi - 2.0_dp*atan( 1.0_dp/exp(yy) )
1538 
1539  return
1540  end subroutine mapprojection_xy2lonlat_mercator
1541 
1542  !-----------------------------------------------------------------------------
1544  subroutine mapprojection_lonlat2xy_mercator( &
1545  lon, lat, &
1546  param, &
1547  x, y )
1548  use scale_const, only: &
1549  undef => const_undef
1550  implicit none
1551  real(rp), intent(in) :: lon, lat ! [rad]
1552  type(mappingparam), intent(in) :: param
1553  real(rp), intent(out) :: x, y
1554 
1555  real(dp) :: dlon, latrot, dist
1556  !---------------------------------------------------------------------------
1557 
1558  if ( lon==undef .or. lat==undef ) then
1559  x = undef
1560  y = undef
1561 
1562  return
1563  end if
1564 
1565  dlon = lon - param%basepoint_lon
1566 
1567  latrot = 0.5_dp*pi - lat
1568 
1569  dist = 1.0_dp / tan(0.5_dp*latrot)
1570 
1571  x = param%x + radius * param%fact * dlon
1572  y = param%y + radius * param%fact * log(dist)
1573 
1574  return
1575  end subroutine mapprojection_lonlat2xy_mercator
1576 
1577  !-----------------------------------------------------------------------------
1579  subroutine mapprojection_mapfactor_mercator( &
1580  lat, &
1581  param, &
1582  m1, m2 )
1583  implicit none
1584  real(rp), intent(in) :: lat ! [rad]
1585  type(mappingparam), intent(in) :: param
1586  real(rp), intent(out) :: m1, m2
1587  !---------------------------------------------------------------------------
1588 
1589  m1 = param%fact / cos(real(lat,kind=dp))
1590  m2 = m1
1591 
1592  return
1593  end subroutine mapprojection_mapfactor_mercator
1594 
1595  !-----------------------------------------------------------------------------
1596  subroutine mapprojection_rotcoef_mercator( &
1597  lon, lat, &
1598  param, &
1599  rotc_cos, rotc_sin )
1600  implicit none
1601  real(rp), intent(in) :: lon, lat
1602  type(mappingparam), intent(in) :: param
1603  real(rp), intent(out) :: rotc_cos, rotc_sin
1604  !---------------------------------------------------------------------------
1605 
1606  rotc_cos = param%rot_fact_cos
1607  rotc_sin = param%rot_fact_sin
1608 
1609  return
1610  end subroutine mapprojection_rotcoef_mercator
1611 
1612  !-----------------------------------------------------------------------------
1615  info, &
1616  param )
1617  use scale_prc, only: &
1618  prc_abort
1619  implicit none
1620  type(mappinginfo), intent(in) :: info
1621  type(mappingparam), intent(out) :: param
1622 
1623  real(dp) :: ec_lat
1624  real(dp) :: basepoint_lat
1625  real(dp) :: basepoint_x, basepoint_y
1626 
1627  real(dp) :: lat0
1628  !---------------------------------------------------------------------------
1629 
1630  ec_lat = info%standard_parallel(1)
1631 
1632  basepoint_lat = info%latitude_of_projection_origin
1633  basepoint_x = info%false_easting
1634  basepoint_y = info%false_northing
1635 
1636  lat0 = ec_lat * d2r
1637 
1638  ! pre-calc factor
1639  param%fact = cos(lat0)
1640 
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
1643  call prc_abort
1644  endif
1645 
1646  param%x = basepoint_x
1647  param%y = basepoint_y - radius * basepoint_lat * d2r
1648 
1649  log_newline
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
1654 
1655  return
1657 
1658  !-----------------------------------------------------------------------------
1661  x, y, &
1662  param, &
1663  lon, lat )
1664  use scale_prc, only: &
1665  prc_abort
1666  use scale_const, only: &
1667  undef => const_undef
1668  implicit none
1669  real(rp), intent(in) :: x, y
1670  type(mappingparam), intent(in) :: param
1671  real(rp), intent(out) :: lon, lat ! [rad]
1672 
1673  real(dp) :: xx, yy
1674  !---------------------------------------------------------------------------
1675 
1676  if ( x==undef .or. y==undef ) then
1677  lon = undef
1678  lat = undef
1679 
1680  return
1681  end if
1682 
1683  call xy_rotate( x, y, & ! (in)
1684  param, & ! (in)
1685  xx, yy ) ! (out)
1686 
1687  xx = ( xx - param%x ) / radius / param%fact
1688  yy = ( yy - param%y ) / radius
1689 
1690  lon = xx + param%basepoint_lon
1691  lat = yy
1692 
1693  if ( abs(lat) > 0.5_dp*pi ) then
1694  log_error("MAPPROJECTION_EquidistantCylindrical_xy2lonlat",*) 'Invalid latitude range! value=', lat
1695  call prc_abort
1696  endif
1697 
1698  return
1700 
1701  !-----------------------------------------------------------------------------
1704  lon, lat, &
1705  param, &
1706  x, y )
1707  use scale_const, only: &
1708  undef => const_undef
1709  implicit none
1710  real(rp), intent(in) :: lon, lat ! [rad]
1711  type(mappingparam), intent(in) :: param
1712  real(rp), intent(out) :: x, y
1713 
1714  real(dp) :: dlon
1715  real(dp) :: xx, yy
1716  !---------------------------------------------------------------------------
1717 
1718  if ( lon==undef .or. lat==undef ) then
1719  x = undef
1720  y = undef
1721 
1722  return
1723  end if
1724 
1725  dlon = lon - param%basepoint_lon
1726 
1727  xx = param%x + radius * param%fact * dlon
1728  yy = param%y + radius * lat
1729 
1730  call xy_unrotate( xx, yy, & ! (in)
1731  param, & ! (in)
1732  x, y ) ! (out)
1733 
1734  return
1736 
1737  !-----------------------------------------------------------------------------
1740  lat, &
1741  param, &
1742  m1, m2 )
1743  implicit none
1744  real(rp), intent(in) :: lat ! [rad]
1745  type(mappingparam), intent(in) :: param
1746  real(rp), intent(out) :: m1, m2
1747 
1748  real(dp) :: mm1, mm2
1749  !---------------------------------------------------------------------------
1750 
1751  mm1 = param%fact / cos(real(lat,kind=dp))
1752  mm2 = 1.0_dp
1753 
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 )
1756 
1757  return
1759 
1760  !-----------------------------------------------------------------------------
1762  lon, lat, &
1763  param, &
1764  rotc_cos, rotc_sin )
1765  implicit none
1766  real(rp), intent(in) :: lon, lat
1767  type(mappingparam), intent(in) :: param
1768  real(rp), intent(out) :: rotc_cos, rotc_sin
1769  !---------------------------------------------------------------------------
1770 
1771  rotc_cos = param%rot_fact_cos
1772  rotc_sin = param%rot_fact_sin
1773 
1774  return
1776 
1777 ! private
1778 
1779  subroutine xy_rotate( &
1780  x, y, &
1781  param, &
1782  xx, yy )
1783  implicit none
1784  real(rp), intent(in) :: x, y
1785  type(mappingparam), intent(in) :: param
1786  real(dp), intent(out) :: xx, yy
1787 
1788  real(dp) :: xd, yd
1789 
1790  xd = x - param%basepoint_x
1791  yd = y - param%basepoint_y
1792 
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
1797 
1798  return
1799  end subroutine xy_rotate
1800 
1801  subroutine xy_unrotate( &
1802  xx, yy, &
1803  param, &
1804  x, y )
1805  implicit none
1806  real(dp), intent(in) :: xx, yy
1807  type(mappingparam), intent(in) :: param
1808  real(rp), intent(out) :: x, y
1809 
1810  real(dp) :: xxd, yyd
1811 
1812  xxd = xx - param%basepoint_x
1813  yyd = yy - param%basepoint_y
1814 
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 )
1819 
1820  return
1821  end subroutine xy_unrotate
1822 
1823 end module scale_mapprojection
scale_mapprojection::lonlat2xy_s
Definition: scale_mapprojection.F90:162
scale_mapprojection::mapprojection_xy2lonlat_0d_initialized
subroutine mapprojection_xy2lonlat_0d_initialized(x, y, lon, lat)
(x,y) -> (lon,lat)
Definition: scale_mapprojection.F90:386
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:350
scale_mapprojection::mapprojection_mapfactor_mercator
subroutine, public mapprojection_mapfactor_mercator(lat, param, m1, m2)
Mercator projection: (lon,lat) -> (m1=m2)
Definition: scale_mapprojection.F90:1583
scale_mapprojection::mapprojection_get_param_lambertconformal
subroutine, public mapprojection_get_param_lambertconformal(info, param)
Lambert Conformal projection.
Definition: scale_mapprojection.F90:1084
scale_mapprojection::mapfactor_s
Definition: scale_mapprojection.F90:169
scale_mapprojection::mapprojection_xy2lonlat_none
subroutine, public mapprojection_xy2lonlat_none(x, y, param, lon, lat)
No projection, lon,lat are determined by gnomonic projection: (x,y) -> (lon,lat)
Definition: scale_mapprojection.F90:937
scale_precision
module PRECISION
Definition: scale_precision.F90:14
scale_mapprojection::mapprojection_get_param_none
subroutine, public mapprojection_get_param_none
No projection.
Definition: scale_mapprojection.F90:925
scale_mapprojection::mapprojection_lonlat2xy_equidistantcylindrical
subroutine, public mapprojection_lonlat2xy_equidistantcylindrical(lon, lat, param, x, y)
Equidistant Cylindrical projection: (lon,lat) -> (x,y)
Definition: scale_mapprojection.F90:1707
scale_mapprojection::mapprojection_xy2lonlat_2d_param
subroutine mapprojection_xy2lonlat_2d_param(IA, IS, IE, JA, JS, JE, x, y, mapping_name, mapping_param, lon, lat)
Definition: scale_mapprojection.F90:474
scale_mapprojection::mapprojection_get_param_mercator
subroutine, public mapprojection_get_param_mercator(info, param)
Mercator projection.
Definition: scale_mapprojection.F90:1458
scale_mapprojection::mapprojection_lonlat2xy_none
subroutine, public mapprojection_lonlat2xy_none(lon, lat, param, x, y)
Definition: scale_mapprojection.F90:990
scale_mapprojection::mapprojection_xy2lonlat_lambertconformal
subroutine, public mapprojection_xy2lonlat_lambertconformal(x, y, param, lon, lat)
Lambert Conformal projection: (x,y) -> (lon,lat)
Definition: scale_mapprojection.F90:1155
scale_mapprojection::mapprojection_rotcoef_equidistantcylindrical
subroutine, public mapprojection_rotcoef_equidistantcylindrical(lon, lat, param, rotc_cos, rotc_sin)
Definition: scale_mapprojection.F90:1765
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_mapprojection::mapprojection_lonlat2xy_mercator
subroutine, public mapprojection_lonlat2xy_mercator(lon, lat, param, x, y)
Mercator projection: (lon,lat) -> (x,y)
Definition: scale_mapprojection.F90:1548
scale_mapprojection::mapprojection_basepoint_lat
real(rp), public mapprojection_basepoint_lat
Definition: scale_mapprojection.F90:91
scale_precision::rp
integer, parameter, public rp
Definition: scale_precision.F90:41
scale_mapprojection
module Map projection
Definition: scale_mapprojection.F90:12
scale_mapprojection::rotcoef_s
Definition: scale_mapprojection.F90:176
scale_io
module STDIO
Definition: scale_io.F90:10
scale_mapprojection::mapprojection_mapfactor_lambertconformal
subroutine, public mapprojection_mapfactor_lambertconformal(lat, param, m1, m2)
Lambert Conformal projection: (lon,lat) -> (m1=m2)
Definition: scale_mapprojection.F90:1236
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_mapprojection::mapprojection_xy2lonlat_polarstereographic
subroutine, public mapprojection_xy2lonlat_polarstereographic(x, y, param, lon, lat)
Polar Stereographic projection: (x,y) -> (lon,lat)
Definition: scale_mapprojection.F90:1336
scale_mapprojection::mapprojection_mapfactor_none
subroutine, public mapprojection_mapfactor_none(lat, param, m1, m2)
No projection: m1=m2=1.
Definition: scale_mapprojection.F90:1047
scale_mapprojection::xy2lonlat
procedure(xy2lonlat_s), pointer xy2lonlat
Definition: scale_mapprojection.F90:184
scale_prc_cartesc
module process / cartesC
Definition: scale_prc_cartesC.F90:11
scale_mapprojection::mapprojection_setup
subroutine, public mapprojection_setup(DOMAIN_CENTER_X, DOMAIN_CENTER_Y)
Setup.
Definition: scale_mapprojection.F90:194
scale_mapprojection::mapprojection_rotcoef_initialized
subroutine mapprojection_rotcoef_initialized(IA, IS, IE, JA, JS, JE, lon, lat, rotc_cos, rotc_sin)
u(lat,lon) = cos u(x,y) - sin v(x,y) v(lat,lon) = sin u(x,y) + cos v(x,y)
Definition: scale_mapprojection.F90:821
scale_prof
module profiler
Definition: scale_prof.F90:11
scale_mapprojection::mapprojection_rotcoef_lambertconformal
subroutine, public mapprojection_rotcoef_lambertconformal(lon, lat, param, rotc_cos, rotc_sin)
Definition: scale_mapprojection.F90:1257
scale_mapprojection::mapprojection_mapfactor_polarstereographic
subroutine, public mapprojection_mapfactor_polarstereographic(lat, param, m1, m2)
Polar Stereographic projection: (lon,lat) -> (m1=m2)
Definition: scale_mapprojection.F90:1414
scale_precision::dp
integer, parameter, public dp
Definition: scale_precision.F90:32
scale_mapprojection::mapprojection_basepoint_lon
real(rp), public mapprojection_basepoint_lon
Definition: scale_mapprojection.F90:90
scale_const::const_pi
real(rp), parameter, public const_pi
pi
Definition: scale_const.F90:32
scale_mapprojection::mapprojection_mapfactor_equidistantcylindrical
subroutine, public mapprojection_mapfactor_equidistantcylindrical(lat, param, m1, m2)
Equidistant Cylindrical projection: (lon,lat) -> (m1,m2)
Definition: scale_mapprojection.F90:1743
scale_mapprojection::mapprojection_xy2lonlat_mercator
subroutine, public mapprojection_xy2lonlat_mercator(x, y, param, lon, lat)
Mercator projection: (x,y) -> (lon,lat)
Definition: scale_mapprojection.F90:1511
scale_mapprojection::mapprojection_mappinginfo
type(mappinginfo), save, public mapprojection_mappinginfo
Definition: scale_mapprojection.F90:104
scale_mapprojection::mapprojection_mapfactor_initialized
subroutine mapprojection_mapfactor_initialized(IA, IS, IE, JA, JS, JE, lat, m1, m2)
(x,y) -> (lon,lat)
Definition: scale_mapprojection.F90:715
scale_mapprojection::mapprojection_get_param
subroutine, public mapprojection_get_param(info, param)
Definition: scale_mapprojection.F90:337
scale_const::const_radius
real(rp), public const_radius
radius of the planet [m]
Definition: scale_const.F90:47
scale_mapprojection::mapprojection_rotcoef_polarstereographic
subroutine, public mapprojection_rotcoef_polarstereographic(lon, lat, param, rotc_cos, rotc_sin)
Definition: scale_mapprojection.F90:1431
scale_mapprojection::mapprojection_lonlat2xy_polarstereographic
subroutine, public mapprojection_lonlat2xy_polarstereographic(lon, lat, param, x, y)
Polar Stereographic projection: (lon,lat) -> (x,y)
Definition: scale_mapprojection.F90:1374
scale_const::const_d2r
real(rp), public const_d2r
degree to radian
Definition: scale_const.F90:33
scale_mapprojection::mapprojection_get_param_polarstereographic
subroutine, public mapprojection_get_param_polarstereographic(info, param)
Polar Stereographic projection.
Definition: scale_mapprojection.F90:1284
scale_mapprojection::mapprojection_lonlat2xy_lambertconformal
subroutine, public mapprojection_lonlat2xy_lambertconformal(lon, lat, param, x, y)
Lambert Conformal projection: (lon,lat) -> (x,y)
Definition: scale_mapprojection.F90:1194
scale_mapprojection::mapprojection_xy2lonlat_equidistantcylindrical
subroutine, public mapprojection_xy2lonlat_equidistantcylindrical(x, y, param, lon, lat)
Equidistant Cylindrical projection: (x,y) -> (lon,lat)
Definition: scale_mapprojection.F90:1664
scale_mapprojection::mapprojection_lonlat2xy_2d_param
subroutine mapprojection_lonlat2xy_2d_param(IA, IS, IE, JA, JS, JE, lon, lat, mapping_name, mapping_param, x, y)
Definition: scale_mapprojection.F90:638
scale_mapprojection::mapprojection_get_param_equidistantcylindrical
subroutine, public mapprojection_get_param_equidistantcylindrical(info, param)
Equidistant Cylindrical projection.
Definition: scale_mapprojection.F90:1617
scale_mapprojection::mapprojection_rotcoef_none
subroutine, public mapprojection_rotcoef_none(lon, lat, param, rotc_cos, rotc_sin)
No projection:
Definition: scale_mapprojection.F90:1067
scale_const::const_undef
real(rp), public const_undef
Definition: scale_const.F90:43
scale_io::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:57
scale_mapprojection::mappingparam
Definition: scale_mapprojection.F90:106
scale_prc_cartesc::prc_twod
logical, public prc_twod
2D experiment
Definition: scale_prc_cartesC.F90:56
scale_mapprojection::mapprojection_rotcoef_mercator
subroutine, public mapprojection_rotcoef_mercator(lon, lat, param, rotc_cos, rotc_sin)
Definition: scale_mapprojection.F90:1600
scale_mapprojection::mappinginfo
Definition: scale_mapprojection.F90:93
scale_mapprojection::mapprojection_lonlat2xy_0d_initialized
subroutine mapprojection_lonlat2xy_0d_initialized(lon, lat, x, y)
(lon,lat) -> (x,y)
Definition: scale_mapprojection.F90:550