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
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 = 0.0_dp ! standard latitude1 for Mer. projection [deg]
146  real(dp), private :: mapprojection_ec_lat = 0.0_dp ! 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 
235 
236  !--- read namelist
237  rewind(io_fid_conf)
238  read(io_fid_conf,nml=param_mapprojection,iostat=ierr)
239 
240  if( ierr < 0 ) then !--- missing
241  log_info("MAPPROJECTION_setup",*) 'Not found namelist. Default used.'
242  elseif( ierr > 0 ) then !--- fatal error
243  log_error("MAPPROJECTION_setup",*) 'Not appropriate names in namelist PARAM_MAPPROJECTION. Check!'
244  call prc_abort
245  endif
246  log_nml(param_mapprojection)
247 
248  if ( prc_twod .and. mapprojection_type .ne. 'NONE' ) then
249  log_error("MAPPROJECTION_setup",*) 'MAPPROJECTION_type must be "NONE" for 2D experiment'
250  call prc_abort
251  end if
252 
253  log_newline
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)
258 
259  mapprojection_mappinginfo%longitude_of_central_meridian = undef
260  mapprojection_mappinginfo%straight_vertical_longitude_from_pole = undef
261  mapprojection_mappinginfo%standard_parallel(:) = undef
262 
263  mapprojection_mappinginfo%longitude_of_projection_origin = mapprojection_basepoint_lon
264  mapprojection_mappinginfo%latitude_of_projection_origin = mapprojection_basepoint_lat
265  mapprojection_mappinginfo%false_easting = mapprojection_basepoint_x
266  mapprojection_mappinginfo%false_northing = mapprojection_basepoint_y
267  mapprojection_mappinginfo%rotation = mapprojection_rotation
268 
269 
270  select case(mapprojection_type)
271  case('NONE')
272  log_info_cont(*) '=> NO map projection'
273  mapprojection_mappinginfo%mapping_name = ""
274 
276  lonlat2xy => mapprojection_lonlat2xy_none
277  mapfactor => mapprojection_mapfactor_none
278  rotcoef => mapprojection_rotcoef_none
279  case('LC')
280  log_info_cont(*) '=> Lambert Conformal projection'
281  mapprojection_mappinginfo%mapping_name = "lambert_conformal_conic"
282  mapprojection_mappinginfo%standard_parallel(:) = (/ mapprojection_lc_lat1, mapprojection_lc_lat2 /)
283  mapprojection_mappinginfo%longitude_of_central_meridian = mapprojection_basepoint_lon
284 
289  case('PS')
290  log_info_cont(*) '=> Polar Stereographic projection'
291  if( mapprojection_ps_lat == undef ) mapprojection_ps_lat = mapprojection_basepoint_lat
292  mapprojection_mappinginfo%mapping_name = "polar_stereographic"
293  mapprojection_mappinginfo%straight_vertical_longitude_from_pole = mapprojection_basepoint_lon
294  mapprojection_mappinginfo%standard_parallel(1) = mapprojection_ps_lat
295 
300  case('MER')
301  log_info_cont(*) '=> Mercator projection'
302  mapprojection_mappinginfo%mapping_name = "mercator"
303  mapprojection_mappinginfo%standard_parallel(1) = mapprojection_m_lat
304 
309  case('EC')
310  log_info_cont(*) '=> Equidistant Cylindrical projection'
311  mapprojection_mappinginfo%mapping_name = "equirectangular"
312  mapprojection_mappinginfo%standard_parallel(1) = mapprojection_ec_lat
313  mapprojection_mappinginfo%longitude_of_central_meridian = mapprojection_basepoint_lon
314 
319  case default
320  log_error("MAPPROJECTION_setup",*) 'Unsupported MAPPROJECTION_type. STOP'
321  call prc_abort
322  endselect
323 
325  mapprojection_mappingparam ) ! (out)
326 
327  return
328  end subroutine mapprojection_setup
329 
330  subroutine mapprojection_get_param( &
331  info, &
332  param )
333  use scale_prc, only: &
334  prc_abort
335  implicit none
336  type(mappinginfo), intent(in) :: info
337  type(mappingparam), intent(out) :: param
338 
339  select case( info%mapping_name )
340  case( "" )
342 
343  param%basepoint_lon = info%longitude_of_projection_origin * d2r
344  case( "lambert_conformal_conic" )
345  call mapprojection_get_param_lambertconformal( info, & ! (in)
346  param ) ! (out)
347  param%basepoint_lon = info%longitude_of_central_meridian * d2r
348  case( "polar_stereographic" )
349  call mapprojection_get_param_polarstereographic( info, & ! (in)
350  param ) ! (out)
351  param%basepoint_lon = info%straight_vertical_longitude_from_pole * d2r
352  case( "mercator" )
353  call mapprojection_get_param_mercator( info, & ! (in)
354  param ) ! (out)
355  param%basepoint_lon = info%longitude_of_projection_origin * d2r
356  case( "equirectangular" )
358  param ) ! (out)
359  param%basepoint_lon = info%longitude_of_central_meridian * d2r
360  case default
361  log_error("MAPPROJECTION_set_param",*) 'Unsupported mapping type. STOP'
362  call prc_abort
363  endselect
364 
365  param%basepoint_x = info%false_easting
366  param%basepoint_y = info%false_northing
367 
368  param%basepoint_lat = info%latitude_of_projection_origin * d2r
369 
370  param%rotation = info%rotation * d2r
371  param%rot_fact_sin = sin(param%rotation)
372  param%rot_fact_cos = cos(param%rotation)
373 
374  return
375  end subroutine mapprojection_get_param
376 
377  !-----------------------------------------------------------------------------
380  x, y, &
381  lon, lat )
382  implicit none
383  real(RP), intent(in) :: x, y
384  real(RP), intent(out) :: lon, lat ! [rad]
385  !---------------------------------------------------------------------------
386 
387  call xy2lonlat( x, y, & ! (in)
388  mapprojection_mappingparam, & ! (in)
389  lon, lat ) ! (out)
390 
391  return
393 
394  subroutine mapprojection_xy2lonlat_2d_initialized( &
395  IA, IS, IE, JA, JS, JE, &
396  x, y, &
397  lon, lat )
398  implicit none
399  integer, intent(in) :: IA, IS, IE
400  integer, intent(in) :: JA, JS, JE
401 
402  real(RP), intent(in) :: x(IA,JA)
403  real(RP), intent(in) :: y(IA,JA)
404  real(RP), intent(out) :: lon(IA,JA) ! [rad]
405  real(RP), intent(out) :: lat(IA,JA) ! [rad]
406 
407  integer :: i, j
408  !---------------------------------------------------------------------------
409 
410  !$omp parallel do
411  do j = js, je
412  do i = is, ie
413  call mapprojection_xy2lonlat_0d_initialized( x(i,j), y(i,j), lon(i,j), lat(i,j) )
414  end do
415  end do
416 
417  return
418  end subroutine mapprojection_xy2lonlat_2d_initialized
419 
420  subroutine mapprojection_xy2lonlat_0d_param( &
421  x, y, &
422  mapping_name, &
423  mapping_param, &
424  lon, lat )
425  use scale_prc, only: &
426  prc_abort
427  implicit none
428  real(RP), intent(in) :: x, y
429  character(len=*), intent(in) :: mapping_name
430  type(mappingparam), intent(in) :: mapping_param
431 
432  real(RP), intent(out) :: lon, lat ! [rad]
433  !---------------------------------------------------------------------------
434 
435  select case( mapping_name )
436  case( "" )
437  call mapprojection_xy2lonlat_none( x, y, & ! (in)
438  mapping_param, & ! (in)
439  lon, lat ) ! (out)
440  case( "lambert_conformal_conic" )
441  call mapprojection_xy2lonlat_lambertconformal( x, y, & ! (in)
442  mapping_param, & ! (in)
443  lon, lat ) ! (out)
444  case( "polar_stereographic" )
445  call mapprojection_xy2lonlat_polarstereographic( x, y, & ! (in)
446  mapping_param, & ! (in)
447  lon, lat ) ! (out)
448  case( "mercator" )
449  call mapprojection_xy2lonlat_mercator( x, y, & ! (in)
450  mapping_param, & ! (in)
451  lon, lat ) ! (out)
452  case( "equirectangular" )
454  mapping_param, & ! (in)
455  lon, lat ) ! (out)
456  case default
457  log_error("MAPPROJECTION_xy2lonlat_0D_param",*) 'Unsupported mapping type. STOP: ', trim(mapping_name)
458  call prc_abort
459  endselect
460 
461  return
462  end subroutine mapprojection_xy2lonlat_0d_param
463 
465  IA, IS, IE, JA, JS, JE, &
466  x, y, &
467  mapping_name, &
468  mapping_param, &
469  lon, lat )
470  use scale_prc, only: &
471  prc_abort
472  implicit none
473  integer, intent(in) :: IA, IS, IE
474  integer, intent(in) :: JA, JS, JE
475 
476  real(RP), intent(in) :: x(IA,JA)
477  real(RP), intent(in) :: y(IA,JA)
478  character(len=*), intent(in) :: mapping_name
479  type(mappingparam), intent(in) :: mapping_param
480 
481  real(RP), intent(out) :: lon(IA,JA) ! [rad]
482  real(RP), intent(out) :: lat(IA,JA) ! [rad]
483 
484  integer :: i, j
485  !---------------------------------------------------------------------------
486 
487  select case( mapping_name )
488  case( "" )
489  !$omp parallel do
490  do j = js, je
491  do i = is, ie
492  call mapprojection_xy2lonlat_none( x(i,j), y(i,j), & ! (in)
493  mapping_param, & ! (in)
494  lon(i,j), lat(i,j) ) ! (out)
495  end do
496  end do
497  case( "lambert_conformal_conic" )
498  !$omp parallel do
499  do j = js, je
500  do i = is, ie
501  call mapprojection_xy2lonlat_lambertconformal( x(i,j), y(i,j), & ! (in)
502  mapping_param, & ! (in)
503  lon(i,j), lat(i,j) ) ! (out)
504  end do
505  end do
506  case( "polar_stereographic" )
507  !$omp parallel do
508  do j = js, je
509  do i = is, ie
510  call mapprojection_xy2lonlat_polarstereographic( x(i,j), y(i,j), & ! (in)
511  mapping_param, & ! (in)
512  lon(i,j), lat(i,j) ) ! (out)
513  end do
514  end do
515  case( "mercator" )
516  !$omp parallel do
517  do j = js, je
518  do i = is, ie
519  call mapprojection_xy2lonlat_mercator( x(i,j), y(i,j), & ! (in)
520  mapping_param, & ! (in)
521  lon(i,j), lat(i,j) ) ! (out)
522  end do
523  end do
524  case( "equirectangular" )
525  !$omp parallel do
526  do j = js, je
527  do i = is, ie
528  call mapprojection_xy2lonlat_equidistantcylindrical( x(i,j), y(i,j), & ! (in)
529  mapping_param, & ! (in)
530  lon(i,j), lat(i,j) ) ! (out)
531  end do
532  end do
533  case default
534  log_error("MAPPROJECTION_xy2lonlat_2D_param",*) 'Unsupported mapping type. STOP: ', trim(mapping_name)
535  call prc_abort
536  endselect
537 
538  return
539  end subroutine mapprojection_xy2lonlat_2d_param
540 
541  !-----------------------------------------------------------------------------
544  lon, lat, &
545  x, y )
546  implicit none
547  real(RP), intent(in) :: lon, lat ! [rad]
548  real(RP), intent(out) :: x, y
549  !---------------------------------------------------------------------------
550 
551  call lonlat2xy( lon, lat, & ! (in)
552  mapprojection_mappingparam, & ! (in)
553  x, y ) ! (out)
554 
555  return
557 
558  subroutine mapprojection_lonlat2xy_2d_initialized( &
559  IA, IS, IE, JA, JS, JE, &
560  lon, lat, &
561  x, y )
562  implicit none
563  integer, intent(in) :: IA, IS, IE
564  integer, intent(in) :: JA, JS, JE
565 
566  real(RP), intent(in) :: lon(IA,JA) ! [rad]
567  real(RP), intent(in) :: lat(IA,JA) ! [rad]
568  real(RP), intent(out) :: x(IA,JA)
569  real(RP), intent(out) :: y(IA,JA)
570 
571  integer :: i, j
572  !---------------------------------------------------------------------------
573 
574  !$omp parallel do
575  do j = js, je
576  do i = is, ie
577  call mapprojection_lonlat2xy_0d_initialized( lon(i,j), lat(i,j), x(i,j), y(i,j) )
578  end do
579  end do
580 
581  return
582  end subroutine mapprojection_lonlat2xy_2d_initialized
583 
584  subroutine mapprojection_lonlat2xy_0d_param( &
585  lon, lat, &
586  mapping_name, &
587  mapping_param, &
588  x, y )
589  use scale_prc, only: &
590  prc_abort
591  implicit none
592  real(RP), intent(in) :: lon, lat ! [rad]
593  character(len=*), intent(in) :: mapping_name
594  type(mappingparam), intent(in) :: mapping_param
595 
596  real(RP), intent(out) :: x, y
597  !---------------------------------------------------------------------------
598 
599  select case( mapping_name )
600  case( "" )
601  call mapprojection_lonlat2xy_none( lon, lat, & ! (in)
602  mapping_param, & ! (in)
603  x, y ) ! (out)
604  case( "lambert_conformal_conic" )
605  call mapprojection_lonlat2xy_lambertconformal( lon, lat, & ! (in)
606  mapping_param, & ! (in)
607  x, y ) ! (out)
608  case( "polar_stereographic" )
609  call mapprojection_lonlat2xy_polarstereographic( lon, lat, & ! (in)
610  mapping_param, & ! (in)
611  x, y ) ! (out)
612  case( "mercator" )
613  call mapprojection_lonlat2xy_mercator( lon, lat, & ! (in)
614  mapping_param, & ! (in)
615  x, y ) ! (out)
616  case( "equirectangular" )
617  call mapprojection_lonlat2xy_equidistantcylindrical( lon, lat, & ! (in)
618  mapping_param, & ! (in)
619  x, y ) ! (out)
620  case default
621  log_error("MAPPROJECTION_lonlat2xy_0D_param",*) 'Unsupported mapping type. STOP: ', trim(mapping_name)
622  call prc_abort
623  endselect
624 
625  return
626  end subroutine mapprojection_lonlat2xy_0d_param
627 
629  IA, IS, IE, JA, JS, JE, &
630  lon, lat, &
631  mapping_name, &
632  mapping_param, &
633  x, y )
634  use scale_prc, only: &
635  prc_abort
636  implicit none
637  integer, intent(in) :: IA, IS, IE
638  integer, intent(in) :: JA, JS, JE
639 
640  real(RP), intent(in) :: lon(IA,JA) ! [rad]
641  real(RP), intent(in) :: lat(IA,JA) ! [rad]
642  character(len=*), intent(in) :: mapping_name
643  type(mappingparam), intent(in) :: mapping_param
644 
645  real(RP), intent(out) :: x(IA,JA)
646  real(RP), intent(out) :: y(IA,JA)
647 
648  integer :: i, j
649  !---------------------------------------------------------------------------
650 
651  select case( mapping_name )
652  case( "" )
653  !$omp parallel do
654  do j = js, je
655  do i = is, ie
656  call mapprojection_lonlat2xy_none( lon(i,j), lat(i,j), & ! (in)
657  mapping_param, & ! (in)
658  x(i,j), y(i,j) ) ! (out)
659  end do
660  end do
661  case( "lambert_conformal_conic" )
662  !$omp parallel do
663  do j = js, je
664  do i = is, ie
665  call mapprojection_lonlat2xy_lambertconformal( lon(i,j), lat(i,j), & ! (in)
666  mapping_param, & ! (in)
667  x(i,j), y(i,j) ) ! (out)
668  end do
669  end do
670  case( "polar_stereographic" )
671  !$omp parallel do
672  do j = js, je
673  do i = is, ie
674  call mapprojection_lonlat2xy_polarstereographic( lon(i,j), lat(i,j), & ! (in)
675  mapping_param, & ! (in)
676  x(i,j), y(i,j) ) ! (out)
677  end do
678  end do
679  case( "mercator" )
680  !$omp parallel do
681  do j = js, je
682  do i = is, ie
683  call mapprojection_lonlat2xy_mercator( lon(i,j), lat(i,j), & ! (in)
684  mapping_param, & ! (in)
685  x(i,j), y(i,j) ) ! (out)
686  end do
687  end do
688  case( "equirectangular" )
689  !$omp parallel do
690  do j = js, je
691  do i = is, ie
692  call mapprojection_lonlat2xy_equidistantcylindrical( lon(i,j), lat(i,j), & ! (in)
693  mapping_param, & ! (in)
694  x(i,j), y(i,j) ) ! (out)
695  end do
696  end do
697  case default
698  log_error("MAPPROJECTION_lonlat2xy_2D_param",*) 'Unsupported mapping type. STOP: ', trim(mapping_name)
699  call prc_abort
700  endselect
701 
702  return
703  end subroutine mapprojection_lonlat2xy_2d_param
704 
705  !-----------------------------------------------------------------------------
708  IA, IS, IE, JA, JS, JE, &
709  lat, &
710  m1, m2 )
711  implicit none
712  integer, intent(in) :: IA, IS, IE
713  integer, intent(in) :: JA, JS, JE
714 
715  real(RP), intent(in) :: lat(IA,JA) ! [rad]
716  real(RP), intent(out) :: m1 (IA,JA)
717  real(RP), intent(out) :: m2 (IA,JA)
718 
719  integer :: i, j
720  !---------------------------------------------------------------------------
721 
722  !$omp parallel do
723  do j = js, je
724  do i = is, ie
725  call mapfactor( lat(i,j), & ! (in)
726  mapprojection_mappingparam, & ! (in)
727  m1(i,j), m2(i,j) ) ! (out)
728  end do
729  end do
730 
731  return
733 
734  subroutine mapprojection_mapfactor_param( &
735  IA, IS, IE, JA, JS, JE, &
736  lat, &
737  mapping_name, &
738  mapping_param, &
739  m1, m2 )
740  use scale_prc, only: &
741  prc_abort
742  implicit none
743  integer, intent(in) :: IA, IS, IE
744  integer, intent(in) :: JA, JS, JE
745 
746  real(RP), intent(in) :: lat(IA,JA) ! [rad]
747  character(len=*), intent(in) :: mapping_name
748  type(mappingparam), intent(in) :: mapping_param
749 
750  real(RP), intent(out) :: m1 (IA,JA)
751  real(RP), intent(out) :: m2 (IA,JA)
752 
753  integer :: i, j
754  !---------------------------------------------------------------------------
755 
756  select case( mapping_name )
757  case( "" )
758  !$omp parallel do
759  do j = js, je
760  do i = is, ie
761  call mapprojection_mapfactor_none( lat(i,j), & ! (in)
762  mapping_param, & ! (in)
763  m1(i,j), m2(i,j) ) ! (out)
764  end do
765  end do
766  case( "lambert_conformal_conic" )
767  !$omp parallel do
768  do j = js, je
769  do i = is, ie
770  call mapprojection_mapfactor_lambertconformal( lat(i,j), & ! (in)
771  mapping_param, & ! (in)
772  m1(i,j), m2(i,j) ) ! (out)
773  end do
774  end do
775  case( "polar_stereographic" )
776  !$omp parallel do
777  do j = js, je
778  do i = is, ie
779  call mapprojection_mapfactor_polarstereographic( lat(i,j), & ! (in)
780  mapping_param, & ! (in)
781  m1(i,j), m2(i,j) ) ! (out)
782  end do
783  end do
784  case( "mercator" )
785  !$omp parallel do
786  do j = js, je
787  do i = is, ie
788  call mapprojection_mapfactor_mercator( lat(i,j), & ! (in)
789  mapping_param, & ! (in)
790  m1(i,j), m2(i,j) ) ! (out)
791  end do
792  end do
793  case( "equirectangular" )
794  !$omp parallel do
795  do j = js, je
796  do i = is, ie
797  call mapprojection_mapfactor_equidistantcylindrical( lat(i,j), & ! (in)
798  mapping_param, & ! (in)
799  m1(i,j), m2(i,j) ) ! (out)
800  end do
801  end do
802  case default
803  log_error("MAPPROJECTION_mapfactor_param",*) 'Unsupported mapping type. STOP: ', trim(mapping_name)
804  call prc_abort
805  endselect
806 
807  return
808  end subroutine mapprojection_mapfactor_param
809 
810  !-----------------------------------------------------------------------------
814  IA, IS, IE, JA, JS, JE, &
815  lon, lat, &
816  rotc_cos, rotc_sin )
817  implicit none
818  integer, intent(in) :: IA, IS, IE
819  integer, intent(in) :: JA, JS, JE
820 
821  real(RP), intent(in) :: lon(IA,JA) ! [rad]
822  real(RP), intent(in) :: lat(IA,JA) ! [rad]
823  real(RP), intent(out) :: rotc_cos(IA,JA)
824  real(RP), intent(out) :: rotc_sin(IA,JA)
825 
826  integer :: i, j
827  !---------------------------------------------------------------------------
828 
829  !$omp parallel do
830  do j = js, je
831  do i = is, ie
832  call rotcoef( lon(i,j), lat(i,j), & ! (in)
833  mapprojection_mappingparam, & ! (in)
834  rotc_cos(i,j), rotc_sin(i,j) ) ! (out)
835  end do
836  end do
837 
838  return
839  end subroutine mapprojection_rotcoef_initialized
840 
841  subroutine mapprojection_rotcoef_param( &
842  IA, IS, IE, JA, JS, JE, &
843  lon, lat, &
844  mapping_name, &
845  mapping_param, &
846  rotc_cos, rotc_sin )
847  use scale_prc, only: &
848  prc_abort
849  implicit none
850  integer, intent(in) :: IA, IS, IE
851  integer, intent(in) :: JA, JS, JE
852 
853  real(RP), intent(in) :: lon(IA,JA) ! [rad]
854  real(RP), intent(in) :: lat(IA,JA) ! [rad]
855  character(len=*), intent(in) :: mapping_name
856  type(mappingparam), intent(in) :: mapping_param
857 
858  real(RP), intent(out) :: rotc_cos(IA,JA)
859  real(RP), intent(out) :: rotc_sin(IA,JA)
860 
861  integer :: i, j
862  !---------------------------------------------------------------------------
863 
864  select case( mapping_name )
865  case( "" )
866  !$omp parallel do
867  do j = js, je
868  do i = is, ie
869  call mapprojection_rotcoef_none( lon(i,j), lat(i,j), & ! (in)
870  mapping_param, & ! (in)
871  rotc_cos(i,j), rotc_sin(i,j) ) ! (out)
872  end do
873  end do
874  case( "lambert_conformal_conic" )
875  !$omp parallel do
876  do j = js, je
877  do i = is, ie
878  call mapprojection_rotcoef_lambertconformal( lon(i,j), lat(i,j), & ! (in)
879  mapping_param, & ! (in)
880  rotc_cos(i,j), rotc_sin(i,j) ) ! (out)
881  end do
882  end do
883  case( "polar_stereographic" )
884  !$omp parallel do
885  do j = js, je
886  do i = is, ie
887  call mapprojection_rotcoef_polarstereographic( lon(i,j), lat(i,j), & ! (in)
888  mapping_param, & ! (in)
889  rotc_cos(i,j), rotc_sin(i,j) ) ! (out)
890  end do
891  end do
892  case( "mercator" )
893  !$omp parallel do
894  do j = js, je
895  do i = is, ie
896  call mapprojection_rotcoef_mercator( lon(i,j), lat(i,j), & ! (in)
897  mapping_param, & ! (in)
898  rotc_cos(i,j), rotc_sin(i,j) ) ! (out)
899  end do
900  end do
901  case( "equirectangular" )
902  !$omp parallel do
903  do j = js, je
904  do i = is, ie
905  call mapprojection_rotcoef_equidistantcylindrical( lon(i,j), lat(i,j), & ! (in)
906  mapping_param, & ! (in)
907  rotc_cos(i,j), rotc_sin(i,j) ) ! (out)
908  end do
909  end do
910  case default
911  log_error("MAPPROJECTION_rotcoef_param",*) 'Unsupported mapping type. STOP: ', trim(mapping_name)
912  call prc_abort
913  endselect
914 
915  return
916  end subroutine mapprojection_rotcoef_param
917 
918  !-----------------------------------------------------------------------------
921  implicit none
922  !---------------------------------------------------------------------------
923 
924  return
925  end subroutine mapprojection_get_param_none
926 
927  !-----------------------------------------------------------------------------
929  subroutine mapprojection_xy2lonlat_none( &
930  x, y, &
931  param, &
932  lon, lat )
933  use scale_const, only: &
934  undef => const_undef
935  implicit none
936  real(rp), intent(in) :: x, y
937  type(mappingparam), intent(in) :: param
938  real(rp), intent(out) :: lon, lat ! [rad]
939 
940  real(dp) :: xx, yy
941  real(dp) :: gno1, gno2
942  real(dp) :: rho, gmm
943  real(dp) :: gmmc, gmms
944  real(dp) :: latc, lats
945  !---------------------------------------------------------------------------
946 
947  if ( x==undef .or. y==undef ) then
948  lon = undef
949  lat = undef
950 
951  return
952  end if
953 
954  call xy_rotate( x, y, & ! (in)
955  param, & ! (in)
956  xx, yy ) ! (out)
957 
958  gno1 = ( x - param%basepoint_x ) / radius
959  gno2 = ( y - param%basepoint_y ) / radius
960 
961  rho = sqrt( gno1 * gno1 + gno2 * gno2 )
962  gmm = atan( rho )
963 
964  if ( rho == 0.0_dp ) then
965  lon = param%basepoint_lon
966  lat = param%basepoint_lat
967  else
968  gmmc = cos(gmm)
969  gmms = sin(gmm)
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 )
977  endif
978 
979  return
980  end subroutine mapprojection_xy2lonlat_none
981 
982  subroutine mapprojection_lonlat2xy_none( &
983  lon, lat, &
984  param, &
985  x, y )
986  use scale_const, only: &
987  undef => const_undef
988  implicit none
989  real(rp), intent(in) :: lon, lat ! [rad]
990  type(mappingparam), intent(in) :: param
991  real(rp), intent(out) :: x, y
992 
993  real(dp) :: xx, yy
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
997  real(dp) :: cos_gmm
998  !---------------------------------------------------------------------------
999  ! http://mathworld.wolfram.com/GnomonicProjection.html
1000 
1001  if ( lon==undef .or. lat==undef ) then
1002  x = undef
1003  y = undef
1004 
1005  return
1006  end if
1007 
1008  lat_d = real(lat,kind=dp)
1009  lat0_d = param%basepoint_lat
1010 
1011  dlon = lon - param%basepoint_lon
1012 
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)
1017  dlon_c = cos(dlon)
1018  dlon_s = sin(dlon)
1019 
1020  cos_gmm = lat0_d_s * lat_d_s &
1021  + lat0_d_c * lat_d_c * dlon_c
1022 
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
1026 
1027  xx = gno1 * radius + param%basepoint_x
1028  yy = gno2 * radius + param%basepoint_y
1029 
1030  call xy_unrotate( xx, yy, & ! (in)
1031  param, & ! (in)
1032  x, y ) ! (out)
1033 
1034  return
1035  end subroutine mapprojection_lonlat2xy_none
1036 
1037  !-----------------------------------------------------------------------------
1039  subroutine mapprojection_mapfactor_none( &
1040  lat, &
1041  param, &
1042  m1, m2 )
1043  use scale_const, only: &
1044  undef => const_undef
1045  implicit none
1046  real(rp), intent(in) :: lat ! [rad]
1047  type(mappingparam), intent(in) :: param
1048  real(rp), intent(out) :: m1, m2
1049  !---------------------------------------------------------------------------
1050 
1051  m1 = 1.0_rp
1052  m2 = m1
1053 
1054  return
1055  end subroutine mapprojection_mapfactor_none
1056 
1057  !-----------------------------------------------------------------------------
1059  subroutine mapprojection_rotcoef_none( &
1060  lon, lat, &
1061  param, &
1062  rotc_cos, rotc_sin )
1063  implicit none
1064  real(rp), intent(in) :: lon, lat
1065  type(mappingparam), intent(in) :: param
1066  real(rp), intent(out) :: rotc_cos, rotc_sin
1067  !---------------------------------------------------------------------------
1068 
1069  rotc_cos = param%rot_fact_cos
1070  rotc_sin = param%rot_fact_sin
1071 
1072  return
1073  end subroutine mapprojection_rotcoef_none
1074 
1075  !-----------------------------------------------------------------------------
1078  info, &
1079  param )
1080  use scale_prc, only: &
1081  prc_abort
1082  implicit none
1083  type(mappinginfo), intent(in) :: info
1084  type(mappingparam), intent(out) :: param
1085 
1086  real(dp) :: lc_lat1, lc_lat2
1087  real(dp) :: basepoint_lat
1088  real(dp) :: basepoint_x, basepoint_y
1089 
1090  real(dp) :: lat1rot, lat2rot
1091  real(dp) :: dlon, latrot, dist
1092  !---------------------------------------------------------------------------
1093 
1094  lc_lat1 = info%standard_parallel(1)
1095  lc_lat2 = info%standard_parallel(2)
1096 
1097  basepoint_lat = info%latitude_of_projection_origin
1098  basepoint_x = info%false_easting
1099  basepoint_y = info%false_northing
1100 
1101  if ( lc_lat1 >= lc_lat2 ) then
1102  log_error("MAPPROJECTION_get_param_LambertConformal",*) 'Please set LC_lat1 < LC_lat2 in degree. STOP'
1103  call prc_abort
1104  endif
1105 
1106  ! check hemisphere: 1=north, -1=south
1107  param%hemisphere = sign(1.0_dp,lc_lat1+lc_lat2)
1108 
1109  lat1rot = 0.5_dp*pi - param%hemisphere * lc_lat1 * d2r
1110  lat2rot = 0.5_dp*pi - param%hemisphere * lc_lat2 * d2r
1111 
1112  ! calc conformal factor c
1113  param%c = ( log( sin(lat1rot) ) - log( sin(lat2rot) ) ) &
1114  / ( log( tan(0.5_dp*lat1rot) ) - log( tan(0.5_dp*lat2rot) ) )
1115 
1116  ! pre-calc factor
1117  param%fact = sin(lat1rot) / param%c / tan(0.5_dp*lat1rot)**param%c
1118 
1119  ! calc (x,y) at pole point
1120  dlon = 0.0_dp
1121 
1122  latrot = 0.5_dp*pi - param%hemisphere * basepoint_lat * d2r
1123 
1124  dist = param%fact * radius * tan(0.5_dp*latrot)**param%c
1125 
1126  param%x = basepoint_x - dist * sin(param%c*dlon)
1127  param%y = basepoint_y + param%hemisphere * dist * cos(param%c*dlon)
1128 
1129  log_newline
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
1138 
1139  return
1141 
1142  !-----------------------------------------------------------------------------
1145  x, y, &
1146  param, &
1147  lon, lat )
1148  use scale_const, only: &
1149  undef => const_undef
1150  implicit none
1151  real(rp), intent(in) :: x, y
1152  type(mappingparam), intent(in) :: param
1153  real(rp), intent(out) :: lon, lat ! [rad]
1154 
1155  real(dp) :: xx, yy, dist
1156  !---------------------------------------------------------------------------
1157 
1158  if ( x==undef .or. y==undef ) then
1159  lon = undef
1160  lat = undef
1161 
1162  return
1163  end if
1164 
1165  call xy_rotate( x, y, & ! (in)
1166  param, & ! (in)
1167  xx, yy ) ! (out)
1168  xx = ( xx - param%x ) / radius / param%fact
1169  yy = - param%hemisphere * ( yy - param%y ) / radius / param%fact
1170 
1171  dist = sqrt( xx*xx + yy*yy )
1172 
1173  lon = param%basepoint_lon + atan2(xx,yy) / param%c
1174 
1175  ! check hemisphere: 1=north, -1=south
1176  lat = param%hemisphere * ( 0.5_dp*pi - 2.0_dp*atan( dist**(1.0_dp/param%c) ) )
1177 
1178  return
1180 
1181  !-----------------------------------------------------------------------------
1184  lon, lat, &
1185  param, &
1186  x, y )
1187  use scale_const, only: &
1188  undef => const_undef
1189  implicit none
1190  real(rp), intent(in) :: lon, lat ! [rad]
1191  type(mappingparam), intent(in) :: param
1192  real(rp), intent(out) :: x, y
1193 
1194  real(dp) :: xx, yy
1195  real(dp) :: dlon, latrot, dist
1196  !---------------------------------------------------------------------------
1197 
1198  if ( lon==undef .or. lat==undef ) then
1199  x = undef
1200  y = undef
1201 
1202  return
1203  end if
1204 
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
1208 
1209  latrot = 0.5_dp*pi - param%hemisphere * lat
1210 
1211  dist = param%fact * radius * tan(0.5_dp*latrot)**param%c
1212 
1213  xx = param%x + dist * sin(param%c*dlon)
1214  yy = param%y - param%hemisphere * dist * cos(param%c*dlon)
1215 
1216  call xy_unrotate( xx, yy, & ! (in)
1217  param, & ! (in)
1218  x, y ) ! (out)
1219 
1220  return
1222 
1223  !-----------------------------------------------------------------------------
1226  lat, &
1227  param, &
1228  m1, m2 )
1229  implicit none
1230  real(rp), intent(in) :: lat ! [rad]
1231  type(mappingparam), intent(in) :: param
1232  real(rp), intent(out) :: m1, m2
1233 
1234  real(dp) :: latrot
1235  !---------------------------------------------------------------------------
1236 
1237  latrot = 0.5_dp*pi - param%hemisphere * lat
1238 
1239  m1 = param%fact / sin(latrot) * param%c * tan(0.5_dp*latrot)**param%c
1240  m2 = m1
1241 
1242  return
1244 
1245  !-----------------------------------------------------------------------------
1247  lon, lat, &
1248  param, &
1249  rotc_cos, rotc_sin )
1250  implicit none
1251  real(rp), intent(in) :: lon, lat
1252  type(mappingparam), intent(in) :: param
1253  real(rp), intent(out) :: rotc_cos, rotc_sin
1254 
1255  real(dp) :: dlon
1256  real(dp) :: alpha
1257  !---------------------------------------------------------------------------
1258 
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
1262 
1263  alpha = - param%c * dlon * param%hemisphere &
1264  + param%rotation
1265 
1266  rotc_cos = cos( alpha )
1267  rotc_sin = sin( alpha )
1268 
1269  return
1271 
1272  !-----------------------------------------------------------------------------
1275  info, &
1276  param )
1277  implicit none
1278  type(mappinginfo), intent(in) :: info
1279  type(mappingparam), intent(out) :: param
1280 
1281  real(dp) :: ps_lat
1282  real(dp) :: basepoint_lat
1283  real(dp) :: basepoint_x, basepoint_y
1284 
1285  real(dp) :: lat0
1286  real(dp) :: dlon, latrot, dist
1287  !---------------------------------------------------------------------------
1288 
1289  ps_lat = info%standard_parallel(1)
1290 
1291  basepoint_lat = info%latitude_of_projection_origin
1292  basepoint_x = info%false_easting
1293  basepoint_y = info%false_northing
1294 
1295  ! check hemisphere: 1=north, -1=south
1296  param%hemisphere = sign(1.0_dp,ps_lat)
1297 
1298  lat0 = param%hemisphere * ps_lat * d2r
1299 
1300  ! pre-calc factor
1301  param%fact = 1.0_dp + sin(lat0)
1302 
1303  ! calc (x,y) at pole point
1304  dlon = 0.0_dp
1305 
1306  latrot = 0.5_dp*pi - param%hemisphere * basepoint_lat * d2r
1307 
1308  dist = param%fact * radius * tan(0.5_dp*latrot)
1309 
1310  param%x = basepoint_x - dist * sin(dlon)
1311  param%y = basepoint_y + param%hemisphere * dist * cos(dlon)
1312 
1313  log_newline
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
1319 
1320  return
1322 
1323  !-----------------------------------------------------------------------------
1326  x, y, &
1327  param, &
1328  lon, lat )
1329  use scale_const, only: &
1330  undef => const_undef
1331  implicit none
1332  real(rp), intent(in) :: x, y
1333  type(mappingparam), intent(in) :: param
1334  real(rp), intent(out) :: lon, lat ! [rad]
1335 
1336  real(dp) :: xx, yy, dist
1337  !---------------------------------------------------------------------------
1338 
1339  if ( x==undef .or. y==undef ) then
1340  lon = undef
1341  lat = undef
1342 
1343  return
1344  end if
1345 
1346  call xy_rotate( x, y, & ! (in)
1347  param, & ! (in)
1348  xx, yy ) ! (out)
1349 
1350  xx = ( xx - param%x ) / radius / param%fact
1351  yy = - param%hemisphere * ( yy - param%y ) / radius / param%fact
1352 
1353  dist = sqrt( xx*xx + yy*yy )
1354 
1355  lon = param%basepoint_lon + atan2(xx,yy)
1356  lat = param%hemisphere * ( 0.5_dp*pi - 2.0_dp*atan(dist) )
1357 
1358  return
1360 
1361  !-----------------------------------------------------------------------------
1364  lon, lat, &
1365  param, &
1366  x, y )
1367  use scale_const, only: &
1368  undef => const_undef
1369  implicit none
1370  real(rp), intent(in) :: lon, lat ! [rad]
1371  type(mappingparam), intent(in) :: param
1372  real(rp), intent(out) :: x, y
1373 
1374  real(dp) :: dlon, latrot, dist
1375  real(dp) :: xx, yy
1376  !---------------------------------------------------------------------------
1377 
1378  if ( lon==undef .or. lat==undef ) then
1379  x = undef
1380  y = undef
1381 
1382  return
1383  end if
1384 
1385  dlon = lon - param%basepoint_lon
1386 
1387  latrot = 0.5_dp*pi - param%hemisphere * lat
1388 
1389  dist = param%fact * radius * tan(0.5_dp*latrot)
1390 
1391  xx = param%x + dist * sin(dlon)
1392  yy = param%y - param%hemisphere * dist * cos(dlon)
1393 
1394  call xy_unrotate( xx, yy, & ! (in)
1395  param, & ! (in)
1396  x, y ) ! (out)
1397 
1398  return
1400 
1401  !-----------------------------------------------------------------------------
1404  lat, &
1405  param, &
1406  m1, m2 )
1407  implicit none
1408  real(rp), intent(in) :: lat ! [rad]
1409  type(mappingparam), intent(in) :: param
1410  real(rp), intent(out) :: m1, m2
1411  !---------------------------------------------------------------------------
1412 
1413  m1 = param%fact / ( 1.0_dp + sin(param%hemisphere*lat) )
1414  m2 = m1
1415 
1416  return
1418 
1419  !-----------------------------------------------------------------------------
1421  lon, lat, &
1422  param, &
1423  rotc_cos, rotc_sin )
1424  implicit none
1425  real(rp), intent(in) :: lon, lat
1426  type(mappingparam), intent(in) :: param
1427  real(rp), intent(out) :: rotc_cos, rotc_sin
1428 
1429  real(dp) :: dlon
1430  real(dp) :: alpha
1431  !---------------------------------------------------------------------------
1432 
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
1436 
1437  alpha = - dlon * param%hemisphere &
1438  + param%rotation
1439 
1440  rotc_cos = cos( alpha )
1441  rotc_sin = sin( alpha )
1442 
1443  return
1445 
1446  !-----------------------------------------------------------------------------
1448  subroutine mapprojection_get_param_mercator( &
1449  info, &
1450  param )
1451  use scale_prc, only: &
1452  prc_abort
1453  implicit none
1454  type(mappinginfo), intent(in) :: info
1455  type(mappingparam), intent(out) :: param
1456 
1457  real(dp) :: m_lat
1458  real(dp) :: basepoint_lat
1459  real(dp) :: basepoint_x, basepoint_y
1460 
1461  real(dp) :: lat0
1462  real(dp) :: latrot, dist
1463  !---------------------------------------------------------------------------
1464 
1465  m_lat = info%standard_parallel(1)
1466 
1467  basepoint_lat = info%latitude_of_projection_origin
1468  basepoint_x = info%false_easting
1469  basepoint_y = info%false_northing
1470 
1471  lat0 = m_lat * d2r
1472 
1473  ! pre-calc factor
1474  param%fact = cos(lat0)
1475 
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
1478  call prc_abort
1479  endif
1480 
1481  ! calc (x,y) at (lon,lat) = (base,0)
1482  latrot = 0.5_dp*pi - basepoint_lat * d2r
1483 
1484  dist = 1.0_dp / tan(0.5_dp*latrot)
1485 
1486  param%x = basepoint_x
1487  param%y = basepoint_y - radius * param%fact * log(dist)
1488 
1489  log_newline
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
1494 
1495  return
1496  end subroutine mapprojection_get_param_mercator
1497 
1498  !-----------------------------------------------------------------------------
1500  subroutine mapprojection_xy2lonlat_mercator( &
1501  x, y, &
1502  param, &
1503  lon, lat )
1504  use scale_const, only: &
1505  undef => const_undef
1506  implicit none
1507  real(rp), intent(in) :: x, y
1508  type(mappingparam), intent(in) :: param
1509  real(rp), intent(out) :: lon, lat ! [rad]
1510 
1511  real(dp) :: xx, yy
1512  !---------------------------------------------------------------------------
1513 
1514  if ( x==undef .or. y==undef ) then
1515  lon = undef
1516  lat = undef
1517 
1518  return
1519  end if
1520 
1521  call xy_rotate( x, y, & ! (in)
1522  param, & ! (in)
1523  xx, yy ) ! (out)
1524 
1525  xx = ( xx - param%x ) / radius / param%fact
1526  yy = ( yy - param%y ) / radius / param%fact
1527 
1528  lon = xx + param%basepoint_lon
1529 
1530  lat = 0.5_dp*pi - 2.0_dp*atan( 1.0_dp/exp(yy) )
1531 
1532  return
1533  end subroutine mapprojection_xy2lonlat_mercator
1534 
1535  !-----------------------------------------------------------------------------
1537  subroutine mapprojection_lonlat2xy_mercator( &
1538  lon, lat, &
1539  param, &
1540  x, y )
1541  use scale_const, only: &
1542  undef => const_undef
1543  implicit none
1544  real(rp), intent(in) :: lon, lat ! [rad]
1545  type(mappingparam), intent(in) :: param
1546  real(rp), intent(out) :: x, y
1547 
1548  real(dp) :: dlon, latrot, dist
1549  !---------------------------------------------------------------------------
1550 
1551  if ( lon==undef .or. lat==undef ) then
1552  x = undef
1553  y = undef
1554 
1555  return
1556  end if
1557 
1558  dlon = lon - param%basepoint_lon
1559 
1560  latrot = 0.5_dp*pi - lat
1561 
1562  dist = 1.0_dp / tan(0.5_dp*latrot)
1563 
1564  x = param%x + radius * param%fact * dlon
1565  y = param%y + radius * param%fact * log(dist)
1566 
1567  return
1568  end subroutine mapprojection_lonlat2xy_mercator
1569 
1570  !-----------------------------------------------------------------------------
1572  subroutine mapprojection_mapfactor_mercator( &
1573  lat, &
1574  param, &
1575  m1, m2 )
1576  implicit none
1577  real(rp), intent(in) :: lat ! [rad]
1578  type(mappingparam), intent(in) :: param
1579  real(rp), intent(out) :: m1, m2
1580  !---------------------------------------------------------------------------
1581 
1582  m1 = param%fact / cos(real(lat,kind=dp))
1583  m2 = m1
1584 
1585  return
1586  end subroutine mapprojection_mapfactor_mercator
1587 
1588  !-----------------------------------------------------------------------------
1589  subroutine mapprojection_rotcoef_mercator( &
1590  lon, lat, &
1591  param, &
1592  rotc_cos, rotc_sin )
1593  implicit none
1594  real(rp), intent(in) :: lon, lat
1595  type(mappingparam), intent(in) :: param
1596  real(rp), intent(out) :: rotc_cos, rotc_sin
1597  !---------------------------------------------------------------------------
1598 
1599  rotc_cos = param%rot_fact_cos
1600  rotc_sin = param%rot_fact_sin
1601 
1602  return
1603  end subroutine mapprojection_rotcoef_mercator
1604 
1605  !-----------------------------------------------------------------------------
1608  info, &
1609  param )
1610  use scale_prc, only: &
1611  prc_abort
1612  implicit none
1613  type(mappinginfo), intent(in) :: info
1614  type(mappingparam), intent(out) :: param
1615 
1616  real(dp) :: ec_lat
1617  real(dp) :: basepoint_lat
1618  real(dp) :: basepoint_x, basepoint_y
1619 
1620  real(dp) :: lat0
1621  !---------------------------------------------------------------------------
1622 
1623  ec_lat = info%standard_parallel(1)
1624 
1625  basepoint_lat = info%latitude_of_projection_origin
1626  basepoint_x = info%false_easting
1627  basepoint_y = info%false_northing
1628 
1629  lat0 = ec_lat * d2r
1630 
1631  ! pre-calc factor
1632  param%fact = cos(lat0)
1633 
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
1636  call prc_abort
1637  endif
1638 
1639  param%x = basepoint_x
1640  param%y = basepoint_y - radius * basepoint_lat * d2r
1641 
1642  log_newline
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
1647 
1648  return
1650 
1651  !-----------------------------------------------------------------------------
1654  x, y, &
1655  param, &
1656  lon, lat )
1657  use scale_prc, only: &
1658  prc_abort
1659  use scale_const, only: &
1660  undef => const_undef
1661  implicit none
1662  real(rp), intent(in) :: x, y
1663  type(mappingparam), intent(in) :: param
1664  real(rp), intent(out) :: lon, lat ! [rad]
1665 
1666  real(dp) :: xx, yy
1667  !---------------------------------------------------------------------------
1668 
1669  if ( x==undef .or. y==undef ) then
1670  lon = undef
1671  lat = undef
1672 
1673  return
1674  end if
1675 
1676  call xy_rotate( x, y, & ! (in)
1677  param, & ! (in)
1678  xx, yy ) ! (out)
1679 
1680  xx = ( xx - param%x ) / radius / param%fact
1681  yy = ( yy - param%y ) / radius
1682 
1683  lon = xx + param%basepoint_lon
1684  lat = yy
1685 
1686  if ( abs(lat) > 0.5_dp*pi ) then
1687  log_error("MAPPROJECTION_EquidistantCylindrical_xy2lonlat",*) 'Invalid latitude range! value=', lat
1688  call prc_abort
1689  endif
1690 
1691  return
1693 
1694  !-----------------------------------------------------------------------------
1697  lon, lat, &
1698  param, &
1699  x, y )
1700  use scale_const, only: &
1701  undef => const_undef
1702  implicit none
1703  real(rp), intent(in) :: lon, lat ! [rad]
1704  type(mappingparam), intent(in) :: param
1705  real(rp), intent(out) :: x, y
1706 
1707  real(dp) :: dlon
1708  real(dp) :: xx, yy
1709  !---------------------------------------------------------------------------
1710 
1711  if ( lon==undef .or. lat==undef ) then
1712  x = undef
1713  y = undef
1714 
1715  return
1716  end if
1717 
1718  dlon = lon - param%basepoint_lon
1719 
1720  xx = param%x + radius * param%fact * dlon
1721  yy = param%y + radius * lat
1722 
1723  call xy_unrotate( xx, yy, & ! (in)
1724  param, & ! (in)
1725  x, y ) ! (out)
1726 
1727  return
1729 
1730  !-----------------------------------------------------------------------------
1733  lat, &
1734  param, &
1735  m1, m2 )
1736  implicit none
1737  real(rp), intent(in) :: lat ! [rad]
1738  type(mappingparam), intent(in) :: param
1739  real(rp), intent(out) :: m1, m2
1740 
1741  real(dp) :: mm1, mm2
1742  !---------------------------------------------------------------------------
1743 
1744  mm1 = param%fact / cos(real(lat,kind=dp))
1745  mm2 = 1.0_dp
1746 
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 )
1749 
1750  return
1752 
1753  !-----------------------------------------------------------------------------
1755  lon, lat, &
1756  param, &
1757  rotc_cos, rotc_sin )
1758  implicit none
1759  real(rp), intent(in) :: lon, lat
1760  type(mappingparam), intent(in) :: param
1761  real(rp), intent(out) :: rotc_cos, rotc_sin
1762  !---------------------------------------------------------------------------
1763 
1764  rotc_cos = param%rot_fact_cos
1765  rotc_sin = param%rot_fact_sin
1766 
1767  return
1769 
1770 ! private
1771 
1772  subroutine xy_rotate( &
1773  x, y, &
1774  param, &
1775  xx, yy )
1776  implicit none
1777  real(rp), intent(in) :: x, y
1778  type(mappingparam), intent(in) :: param
1779  real(dp), intent(out) :: xx, yy
1780 
1781  real(dp) :: xd, yd
1782 
1783  xd = x - param%basepoint_x
1784  yd = y - param%basepoint_y
1785 
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
1790 
1791  return
1792  end subroutine xy_rotate
1793 
1794  subroutine xy_unrotate( &
1795  xx, yy, &
1796  param, &
1797  x, y )
1798  implicit none
1799  real(dp), intent(in) :: xx, yy
1800  type(mappingparam), intent(in) :: param
1801  real(rp), intent(out) :: x, y
1802 
1803  real(dp) :: xxd, yyd
1804 
1805  xxd = xx - param%basepoint_x
1806  yyd = yy - param%basepoint_y
1807 
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 )
1812 
1813  return
1814  end subroutine xy_unrotate
1815 
1816 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:382
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:342
scale_mapprojection::mapprojection_mapfactor_mercator
subroutine, public mapprojection_mapfactor_mercator(lat, param, m1, m2)
Mercator projection: (lon,lat) -> (m1=m2)
Definition: scale_mapprojection.F90:1576
scale_mapprojection::mapprojection_get_param_lambertconformal
subroutine, public mapprojection_get_param_lambertconformal(info, param)
Lambert Conformal projection.
Definition: scale_mapprojection.F90:1080
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:933
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:921
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:1700
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:470
scale_mapprojection::mapprojection_get_param_mercator
subroutine, public mapprojection_get_param_mercator(info, param)
Mercator projection.
Definition: scale_mapprojection.F90:1451
scale_mapprojection::mapprojection_lonlat2xy_none
subroutine, public mapprojection_lonlat2xy_none(lon, lat, param, x, y)
Definition: scale_mapprojection.F90:986
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:1148
scale_mapprojection::mapprojection_rotcoef_equidistantcylindrical
subroutine, public mapprojection_rotcoef_equidistantcylindrical(lon, lat, param, rotc_cos, rotc_sin)
Definition: scale_mapprojection.F90:1758
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_const::const_pi
real(rp), public const_pi
pi
Definition: scale_const.F90:31
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:1541
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:1229
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:1329
scale_mapprojection::mapprojection_mapfactor_none
subroutine, public mapprojection_mapfactor_none(lat, param, m1, m2)
No projection: m1=m2=1.
Definition: scale_mapprojection.F90:1043
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:817
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:1250
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:1407
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_mapprojection::mapprojection_mapfactor_equidistantcylindrical
subroutine, public mapprojection_mapfactor_equidistantcylindrical(lat, param, m1, m2)
Equidistant Cylindrical projection: (lon,lat) -> (m1,m2)
Definition: scale_mapprojection.F90:1736
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:1504
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:711
scale_mapprojection::mapprojection_get_param
subroutine, public mapprojection_get_param(info, param)
Definition: scale_mapprojection.F90:333
scale_const::const_radius
real(rp), public const_radius
radius of the planet [m]
Definition: scale_const.F90:44
scale_mapprojection::mapprojection_rotcoef_polarstereographic
subroutine, public mapprojection_rotcoef_polarstereographic(lon, lat, param, rotc_cos, rotc_sin)
Definition: scale_mapprojection.F90:1424
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:1367
scale_const::const_d2r
real(rp), public const_d2r
degree to radian
Definition: scale_const.F90:32
scale_mapprojection::mapprojection_get_param_polarstereographic
subroutine, public mapprojection_get_param_polarstereographic(info, param)
Polar Stereographic projection.
Definition: scale_mapprojection.F90:1277
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:1187
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:1657
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:634
scale_mapprojection::mapprojection_get_param_equidistantcylindrical
subroutine, public mapprojection_get_param_equidistantcylindrical(info, param)
Equidistant Cylindrical projection.
Definition: scale_mapprojection.F90:1610
scale_mapprojection::mapprojection_rotcoef_none
subroutine, public mapprojection_rotcoef_none(lon, lat, param, rotc_cos, rotc_sin)
No projection:
Definition: scale_mapprojection.F90:1063
scale_const::const_undef
real(rp), public const_undef
Definition: scale_const.F90:41
scale_io::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:56
scale_mapprojection::mappingparam
Definition: scale_mapprojection.F90:106
scale_prc_cartesc::prc_twod
logical, public prc_twod
2D experiment
Definition: scale_prc_cartesC.F90:55
scale_mapprojection::mapprojection_mappinginfo
type(mappinginfo), public mapprojection_mappinginfo
Definition: scale_mapprojection.F90:104
scale_mapprojection::mapprojection_rotcoef_mercator
subroutine, public mapprojection_rotcoef_mercator(lon, lat, param, rotc_cos, rotc_sin)
Definition: scale_mapprojection.F90:1593
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:546