SCALE-RM
scale_mapproj.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
13 !-------------------------------------------------------------------------------
14 module scale_mapproj
15  !-----------------------------------------------------------------------------
16  !
17  !++ used modules
18  !
19  use scale_precision
20  use scale_stdio
21  use scale_prof
22  use scale_grid_index
23 
24  use scale_const, only: &
25  undef => const_undef
26  !-----------------------------------------------------------------------------
27  implicit none
28  private
29  !-----------------------------------------------------------------------------
30  !
31  !++ Public procedure
32  !
33  public :: mprj_setup
34  public :: mprj_xy2lonlat
35  public :: mprj_lonlat2xy
36  public :: mprj_mapfactor
37  public :: mprj_rotcoef
38 
39  interface mprj_rotcoef
40  module procedure mprj_rotcoef_0d
41  module procedure mprj_rotcoef_2d
42  end interface mprj_rotcoef
43 
44  !-----------------------------------------------------------------------------
45  !
46  !++ Public parameters & variables
47  !
48  real(RP), public :: mprj_basepoint_lon = 135.221_rp ! position of base point (domain center) in real world [deg]
49  real(RP), public :: mprj_basepoint_lat = 34.653_rp ! position of base point (domain center) in real world [deg]
50 
51  !-----------------------------------------------------------------------------
52  !
53  !++ Private procedure
54  !
55  private :: mprj_none_setup
56  private :: mprj_lambertconformal_setup
57  private :: mprj_polarstereographic_setup
58  private :: mprj_mercator_setup
59  private :: mprj_equidistantcylindrical_setup
60 
61  private :: mprj_none_xy2lonlat
62  private :: mprj_lambertconformal_xy2lonlat
63  private :: mprj_polarstereographic_xy2lonlat
64  private :: mprj_mercator_xy2lonlat
65  private :: mprj_equidistantcylindrical_xy2lonlat
66 
67  private :: mprj_none_lonlat2xy
68  private :: mprj_lambertconformal_lonlat2xy
69  private :: mprj_polarstereographic_lonlat2xy
70  private :: mprj_mercator_lonlat2xy
71  private :: mprj_equidistantcylindrical_lonlat2xy
72 
73  private :: mprj_none_mapfactor
74  private :: mprj_lambertconformal_mapfactor
75  private :: mprj_polarstereographic_mapfactor
76  private :: mprj_mercator_mapfactor
77  private :: mprj_equidistantcylindrical_mapfactor
78 
79  private :: mprj_none_rotcoef_2d
80  private :: mprj_lambertconformal_rotcoef_2d
81  private :: mprj_polarstereographic_rotcoef_2d
82  private :: mprj_mercator_rotcoef_2d
83  private :: mprj_equidistantcylindrical_rotcoef_2d
84 
85  private :: mprj_none_rotcoef_0d
86  private :: mprj_lambertconformal_rotcoef_0d
87  private :: mprj_polarstereographic_rotcoef_0d
88  private :: mprj_mercator_rotcoef_0d
89  private :: mprj_equidistantcylindrical_rotcoef_0d
90 
91  !-----------------------------------------------------------------------------
92  !
93  !++ Private parameters & variables
94  !
95  character(len=H_SHORT), private :: mprj_type = 'NONE'
96  ! 'NONE'
97  ! 'LC'
98  ! 'PS'
99  ! 'MER'
100  ! 'EC'
101 
102  real(DP), private :: mprj_hemisphere ! hemisphere flag: 1=north, -1=south
103 
104  real(DP), private :: mprj_basepoint_x ! position of base point in the model [m]
105  real(DP), private :: mprj_basepoint_y ! position of base point in the model [m]
106 
107  real(DP), private :: mprj_pole_x ! position of north/south pole in the model [m]
108  real(DP), private :: mprj_pole_y ! position of north/south pole in the model [m]
109  real(DP), private :: mprj_eq_x ! position of equator at the base lon. in the model [m]
110  real(DP), private :: mprj_eq_y ! position of equator at the base lon. in the model [m]
111 
112  real(DP), private :: mprj_rotation = 0.0_dp ! rotation factor (only for 'NONE' type)
113 
114  real(DP), private :: mprj_lc_lat1 = 30.0_dp ! standard latitude1 for L.C. projection [deg]
115  real(DP), private :: mprj_lc_lat2 = 60.0_dp ! standard latitude2 for L.C. projection [deg]
116  real(DP), private :: mprj_lc_c ! conformal factor
117  real(DP), private :: mprj_lc_fact ! pre-calc factor
118 
119  real(DP), private :: mprj_ps_lat ! standard latitude1 for P.S. projection [deg]
120  real(DP), private :: mprj_ps_fact ! pre-calc factor
121 
122  real(DP), private :: mprj_m_lat ! standard latitude1 for Mer. projection [deg]
123  real(DP), private :: mprj_m_fact ! pre-calc factor
124 
125  real(DP), private :: mprj_ec_lat ! standard latitude1 for E.C. projection [deg]
126  real(DP), private :: mprj_ec_fact ! pre-calc factor
127 
128  real(DP), private :: pi
129  real(DP), private :: d2r
130  real(DP), private :: radius
131 
132  !-----------------------------------------------------------------------------
133 contains
134  !-----------------------------------------------------------------------------
136  subroutine mprj_setup( DOMAIN_CENTER_X, DOMAIN_CENTER_Y )
137  use scale_process, only: &
138  prc_mpistop
139  use scale_const, only: &
140  pi_rp => const_pi, &
141  d2r_rp => const_d2r, &
142  radius_rp => const_radius
143  implicit none
144 
145  real(RP), intent(in) :: DOMAIN_CENTER_X
146  real(RP), intent(in) :: DOMAIN_CENTER_Y
147 
148  namelist / param_mapproj / &
149  mprj_basepoint_lon, &
150  mprj_basepoint_lat, &
151  mprj_basepoint_x, &
152  mprj_basepoint_y, &
153  mprj_type, &
154  mprj_rotation, &
155  mprj_lc_lat1, &
156  mprj_lc_lat2, &
157  mprj_ps_lat, &
158  mprj_m_lat, &
159  mprj_ec_lat
160 
161  integer :: ierr
162  !---------------------------------------------------------------------------
163 
164  if( io_l ) write(io_fid_log,*)
165  if( io_l ) write(io_fid_log,*) '+++ Module[MAPPROJ]/Categ[GRID]'
166 
167  pi = real(pi_rp, kind=dp)
168  d2r = real(d2r_rp, kind=dp)
169  radius = real(radius_rp, kind=dp)
170 
171  mprj_basepoint_x = undef
172  mprj_basepoint_y = undef
173  mprj_ps_lat = undef
174  mprj_m_lat = undef
175  mprj_ec_lat = undef
176 
177  mprj_basepoint_x = domain_center_x
178  mprj_basepoint_y = domain_center_y
179 
180  !--- read namelist
181  rewind(io_fid_conf)
182  read(io_fid_conf,nml=param_mapproj,iostat=ierr)
183 
184  if( ierr < 0 ) then !--- missing
185  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
186  elseif( ierr > 0 ) then !--- fatal error
187  write(*,*) 'xxx Not appropriate names in namelist PARAM_MAPPROJ. Check!'
188  call prc_mpistop
189  endif
190  if( io_l ) write(io_fid_log,nml=param_mapproj)
191 
192 
193  if( mprj_ps_lat == undef ) mprj_ps_lat = mprj_basepoint_lat
194  if( mprj_m_lat == undef ) mprj_m_lat = mprj_basepoint_lat
195  if( mprj_ec_lat == undef ) mprj_ec_lat = mprj_basepoint_lat
196 
197  if( io_l ) write(io_fid_log,*)
198  if( io_l ) write(io_fid_log,*) '+++ MPRJ_basepoint_lon:', mprj_basepoint_lon
199  if( io_l ) write(io_fid_log,*) '+++ MPRJ_basepoint_lat:', mprj_basepoint_lat
200  if( io_l ) write(io_fid_log,*) '+++ MPRJ_basepoint_x :', mprj_basepoint_x
201  if( io_l ) write(io_fid_log,*) '+++ MPRJ_basepoint_y :', mprj_basepoint_y
202 
203  if( io_l ) write(io_fid_log,*)
204  if( io_l ) write(io_fid_log,*) '*** Projection type:', trim(mprj_type)
205  select case(trim(mprj_type))
206  case('NONE')
207  if( io_l ) write(io_fid_log,*) '=> NO map projection'
208  call mprj_none_setup
209  case('LC')
210  if( io_l ) write(io_fid_log,*) '=> Lambert Conformal projection'
211  call mprj_lambertconformal_setup
212  case('PS')
213  if( io_l ) write(io_fid_log,*) '=> Polar Stereographic projection'
214  call mprj_polarstereographic_setup
215  case('MER')
216  if( io_l ) write(io_fid_log,*) '=> Mercator projection'
217  call mprj_mercator_setup
218  case('EC')
219  if( io_l ) write(io_fid_log,*) '=> Equidistant Cylindrical projection'
220  call mprj_equidistantcylindrical_setup
221  case default
222  write(*,*) ' xxx Unsupported TYPE. STOP'
223  call prc_mpistop
224  endselect
225 
226  return
227  end subroutine mprj_setup
228 
229  !-----------------------------------------------------------------------------
231  subroutine mprj_xy2lonlat( &
232  x, &
233  y, &
234  lon, &
235  lat )
236  use scale_process, only: &
237  prc_mpistop
238  implicit none
239 
240  real(RP), intent(in) :: x
241  real(RP), intent(in) :: y
242  real(RP), intent(out) :: lon ! [rad]
243  real(RP), intent(out) :: lat ! [rad]
244  !---------------------------------------------------------------------------
245 
246  select case(mprj_type)
247  case('NONE')
248  call mprj_none_xy2lonlat( x, y, lon, lat )
249  case('LC')
250  call mprj_lambertconformal_xy2lonlat( x, y, lon, lat )
251  case('PS')
252  call mprj_polarstereographic_xy2lonlat( x, y, lon, lat )
253  case('MER')
254  call mprj_mercator_xy2lonlat( x, y, lon, lat )
255  case('EC')
256  call mprj_equidistantcylindrical_xy2lonlat( x, y, lon, lat )
257  case default
258  write(*,*) ' xxx Unsupported TYPE. STOP'
259  call prc_mpistop
260  endselect
261 
262  return
263  end subroutine mprj_xy2lonlat
264 
265  !-----------------------------------------------------------------------------
267  subroutine mprj_lonlat2xy( &
268  lon, &
269  lat, &
270  x, &
271  y )
272  use scale_process, only: &
273  prc_mpistop
274  implicit none
275 
276  real(RP), intent(in) :: lon ! [rad]
277  real(RP), intent(in) :: lat ! [rad]
278  real(RP), intent(out) :: x
279  real(RP), intent(out) :: y
280  !---------------------------------------------------------------------------
281 
282  select case(mprj_type)
283  case('NONE')
284  call mprj_none_lonlat2xy( lon, lat, x, y )
285  case('LC')
286  call mprj_lambertconformal_lonlat2xy( lon, lat, x, y )
287  case('PS')
288  call mprj_polarstereographic_lonlat2xy( lon, lat, x, y )
289  case('MER')
290  call mprj_mercator_lonlat2xy( lon, lat, x, y )
291  case('EC')
292  call mprj_equidistantcylindrical_lonlat2xy( lon, lat, x, y )
293  case default
294  write(*,*) ' xxx Unsupported TYPE. STOP'
295  call prc_mpistop
296  endselect
297 
298  return
299  end subroutine mprj_lonlat2xy
300 
301  !-----------------------------------------------------------------------------
303  subroutine mprj_mapfactor( &
304  lat, &
305  m1, &
306  m2 )
307  use scale_process, only: &
308  prc_mpistop
309  implicit none
310 
311  real(RP), intent(in) :: lat(ia,ja) ! [rad]
312  real(RP), intent(out) :: m1 (ia,ja)
313  real(RP), intent(out) :: m2 (ia,ja)
314  !---------------------------------------------------------------------------
315 
316  select case(mprj_type)
317  case('NONE')
318  call mprj_none_mapfactor( lat, m1, m2 )
319  case('LC')
320  call mprj_lambertconformal_mapfactor( lat, m1, m2 )
321  case('PS')
322  call mprj_polarstereographic_mapfactor( lat, m1, m2 )
323  case('MER')
324  call mprj_mercator_mapfactor( lat, m1, m2 )
325  case('EC')
326  call mprj_equidistantcylindrical_mapfactor( lat, m1, m2 )
327  case default
328  write(*,*) ' xxx Unsupported TYPE. STOP'
329  call prc_mpistop
330  endselect
331 
332  return
333  end subroutine mprj_mapfactor
334 
335  !-----------------------------------------------------------------------------
338  subroutine mprj_rotcoef_0d( &
339  rotc, &
340  lon, &
341  lat )
342  use scale_process, only: &
343  prc_mpistop
344  implicit none
345 
346  real(RP), intent(out) :: rotc(2)
347  real(RP), intent(in) :: lon ! [rad]
348  real(RP), intent(in) :: lat ! [rad]
349  !---------------------------------------------------------------------------
350 
351  select case(mprj_type)
352  case('NONE')
353  call mprj_none_rotcoef_0d( rotc )
354  case('LC')
355  call mprj_lambertconformal_rotcoef_0d( rotc, lon, lat )
356  case('PS')
357  call mprj_polarstereographic_rotcoef_0d( rotc, lon, lat )
358  case('MER')
359  call mprj_mercator_rotcoef_0d( rotc )
360  case('EC')
361  call mprj_equidistantcylindrical_rotcoef_0d( rotc )
362  case default
363  write(*,*) ' xxx Unsupported TYPE. STOP'
364  call prc_mpistop
365  endselect
366 
367  return
368  end subroutine mprj_rotcoef_0d
369 
370  !-----------------------------------------------------------------------------
373  subroutine mprj_rotcoef_2d( &
374  rotc, &
375  lon, &
376  lat )
377  use scale_process, only: &
378  prc_mpistop
379  implicit none
380 
381  real(RP), intent(out) :: rotc(ia,ja,2)
382  real(RP), intent(in) :: lon (ia,ja) ! [rad]
383  real(RP), intent(in) :: lat (ia,ja) ! [rad]
384  !---------------------------------------------------------------------------
385 
386  select case(mprj_type)
387  case('NONE')
388  call mprj_none_rotcoef_2d( rotc )
389  case('LC')
390  call mprj_lambertconformal_rotcoef_2d( rotc, lon, lat )
391  case('PS')
392  call mprj_polarstereographic_rotcoef_2d( rotc, lon, lat )
393  case('MER')
394  call mprj_mercator_rotcoef_2d( rotc )
395  case('EC')
396  call mprj_equidistantcylindrical_rotcoef_2d( rotc )
397  case default
398  write(*,*) ' xxx Unsupported TYPE. STOP'
399  call prc_mpistop
400  endselect
401 
402  return
403  end subroutine mprj_rotcoef_2d
404 
405  !-----------------------------------------------------------------------------
407  subroutine mprj_none_setup
408  implicit none
409  !---------------------------------------------------------------------------
410 
411  if( io_l ) write(io_fid_log,*)
412  if( io_l ) write(io_fid_log,*) '+++ MPRJ_rotation:', mprj_rotation
413 
414  return
415  end subroutine mprj_none_setup
416 
417  !-----------------------------------------------------------------------------
419  subroutine mprj_none_xy2lonlat( &
420  x, &
421  y, &
422  lon, &
423  lat )
424  implicit none
425 
426  real(RP), intent(in) :: x
427  real(RP), intent(in) :: y
428  real(RP), intent(out) :: lon ! [rad]
429  real(RP), intent(out) :: lat ! [rad]
430 
431  real(DP) :: gno(2)
432  real(DP) :: rho, gmm
433  !---------------------------------------------------------------------------
434 
435  gno(1) = ( (y-mprj_basepoint_y) * sin(mprj_rotation*d2r) &
436  + (x-mprj_basepoint_x) * cos(mprj_rotation*d2r) ) / radius
437  gno(2) = ( (y-mprj_basepoint_y) * cos(mprj_rotation*d2r) &
438  - (x-mprj_basepoint_x) * sin(mprj_rotation*d2r) ) / radius
439 
440  rho = sqrt( gno(1) * gno(1) + gno(2) * gno(2) )
441  gmm = atan( rho )
442 
443  if ( rho == 0.0_dp ) then
444  lon = mprj_basepoint_lon * d2r
445  lat = mprj_basepoint_lat * d2r
446  else
447  lon = mprj_basepoint_lon * d2r &
448  + atan( gno(1)*sin(gmm) / ( rho *cos(mprj_basepoint_lat*d2r)*cos(gmm) &
449  - gno(2)*sin(mprj_basepoint_lat*d2r)*sin(gmm) ) )
450  lat = asin( sin(mprj_basepoint_lat*d2r)*cos(gmm) &
451  + gno(2)*cos(mprj_basepoint_lat*d2r)*sin(gmm) / rho )
452  endif
453 
454  return
455  end subroutine mprj_none_xy2lonlat
456 
457  !-----------------------------------------------------------------------------
459  subroutine mprj_none_lonlat2xy( &
460  lon, &
461  lat, &
462  x, &
463  y )
464  use scale_process, only: &
465  prc_mpistop
466  implicit none
467 
468  real(RP), intent(in) :: lon ! [rad]
469  real(RP), intent(in) :: lat ! [rad]
470  real(RP), intent(out) :: x
471  real(RP), intent(out) :: y
472 
473  real(DP) :: lat_d
474  real(DP) :: gno(2)
475  real(DP) :: cos_gmm
476  !---------------------------------------------------------------------------
477  ! http://mathworld.wolfram.com/GnomonicProjection.html
478 
479  lat_d = real(lat,kind=dp)
480  cos_gmm = sin(mprj_basepoint_lat*d2r) * sin(lat_d) &
481  + cos(mprj_basepoint_lat*d2r) * cos(lat_d) * cos(lon - mprj_basepoint_lon*d2r)
482 
483  gno(1) = (cos(lat_d) * sin(lon - mprj_basepoint_lon*d2r)) / cos_gmm
484  gno(2) = (cos(mprj_basepoint_lat*d2r) * sin(lat_d) &
485  - sin(mprj_basepoint_lat*d2r) * cos(lat_d) * cos(lon - mprj_basepoint_lon*d2r)) / cos_gmm
486 
487  x = mprj_basepoint_x + (gno(1) * cos(mprj_rotation*d2r) &
488  - gno(2) * sin(mprj_rotation*d2r)) * radius
489  y = mprj_basepoint_y + (gno(1) * sin(mprj_rotation*d2r) &
490  + gno(2) * cos(mprj_rotation*d2r)) * radius
491 
492  return
493  end subroutine mprj_none_lonlat2xy
494 
495  !-----------------------------------------------------------------------------
497  subroutine mprj_none_mapfactor( &
498  lat, &
499  m1, &
500  m2 )
501  implicit none
502 
503  real(RP), intent(in) :: lat(ia,ja) ! [rad]
504  real(RP), intent(out) :: m1 (ia,ja)
505  real(RP), intent(out) :: m2 (ia,ja)
506  !---------------------------------------------------------------------------
507 
508  m1 = 1.0_rp
509  m2 = m1
510 
511  return
512  end subroutine mprj_none_mapfactor
513 
514  !-----------------------------------------------------------------------------
516  subroutine mprj_none_rotcoef_0d( &
517  rotc )
518  implicit none
519 
520  real(RP), intent(out) :: rotc(2)
521  !---------------------------------------------------------------------------
522 
523  rotc(1) = 1.0_rp
524  rotc(2) = 0.0_rp
525 
526  return
527  end subroutine mprj_none_rotcoef_0d
528 
529  !-----------------------------------------------------------------------------
531  subroutine mprj_none_rotcoef_2d( &
532  rotc )
533  implicit none
534 
535  real(RP), intent(out) :: rotc(ia,ja,2)
536  !---------------------------------------------------------------------------
537 
538  rotc(:,:,1) = 1.0_rp
539  rotc(:,:,2) = 0.0_rp
540 
541  return
542  end subroutine mprj_none_rotcoef_2d
543 
544  !-----------------------------------------------------------------------------
546  subroutine mprj_lambertconformal_setup
547  use scale_process, only: &
548  prc_mpistop
549  implicit none
550 
551  real(DP) :: lat1rot, lat2rot
552  real(DP) :: dlon, latrot, dist
553  !---------------------------------------------------------------------------
554 
555  if ( mprj_lc_lat1 >= mprj_lc_lat2 ) then
556  write(*,*) ' xxx Please set MPRJ_LC_lat1 < MPRJ_LC_lat2 in degree. STOP'
557  call prc_mpistop
558  endif
559 
560  ! check hemisphere: 1=north, -1=south
561  mprj_hemisphere = sign(1.0_dp,mprj_lc_lat1+mprj_lc_lat2)
562 
563  lat1rot = 0.5_dp*pi - mprj_hemisphere * mprj_lc_lat1 * d2r
564  lat2rot = 0.5_dp*pi - mprj_hemisphere * mprj_lc_lat2 * d2r
565 
566  ! calc conformal factor c
567  mprj_lc_c = ( log( sin(lat1rot) ) - log( sin(lat2rot) ) ) &
568  / ( log( tan(0.5_dp*lat1rot) ) - log( tan(0.5_dp*lat2rot) ) )
569 
570  ! pre-calc factor
571  mprj_lc_fact = sin(lat1rot) / mprj_lc_c / tan(0.5_dp*lat1rot)**mprj_lc_c
572 
573  ! calc (x,y) at pole point
574  dlon = 0.0_dp
575 
576  latrot = 0.5_dp*pi - mprj_hemisphere * mprj_basepoint_lat * d2r
577 
578  dist = mprj_lc_fact * radius * tan(0.5_dp*latrot)**mprj_lc_c
579 
580  mprj_pole_x = mprj_basepoint_x - dist * sin(mprj_lc_c*dlon)
581  mprj_pole_y = mprj_basepoint_y + mprj_hemisphere * dist * cos(mprj_lc_c*dlon)
582 
583  if( io_l ) write(io_fid_log,*)
584  if( io_l ) write(io_fid_log,*) '+++ MPRJ_LC_lat1 :', mprj_lc_lat1
585  if( io_l ) write(io_fid_log,*) '+++ MPRJ_LC_lat2 :', mprj_lc_lat2
586  if( io_l ) write(io_fid_log,*) '+++ MPRJ_hemisphere:', mprj_hemisphere
587  if( io_l ) write(io_fid_log,*) '+++ MPRJ_LC_c :', mprj_lc_c
588  if( io_l ) write(io_fid_log,*) '+++ MPRJ_LC_fact :', mprj_lc_fact
589  if( io_l ) write(io_fid_log,*) '+++ MPRJ_pole_x :', mprj_pole_x
590  if( io_l ) write(io_fid_log,*) '+++ MPRJ_pole_y :', mprj_pole_y
591 
592  return
593  end subroutine mprj_lambertconformal_setup
594 
595  !-----------------------------------------------------------------------------
597  subroutine mprj_lambertconformal_xy2lonlat( &
598  x, &
599  y, &
600  lon, &
601  lat )
602  implicit none
603 
604  real(RP), intent(in) :: x
605  real(RP), intent(in) :: y
606  real(RP), intent(out) :: lon ! [rad]
607  real(RP), intent(out) :: lat ! [rad]
608 
609  real(DP) :: xx, yy, dist
610  !---------------------------------------------------------------------------
611 
612  xx = ( x - mprj_pole_x ) / radius / mprj_lc_fact
613  yy = -mprj_hemisphere * ( y - mprj_pole_y ) / radius / mprj_lc_fact
614 
615  dist = sqrt( xx*xx + yy*yy )
616 
617  lon = mprj_basepoint_lon * d2r + atan2(xx,yy) / mprj_lc_c
618  lon = mod( lon+2.0_dp*pi, 2.0_dp*pi )
619 
620  ! check hemisphere: 1=north, -1=south
621  lat = mprj_hemisphere * ( 0.5_dp*pi - 2.0_dp*atan( dist**(1.0_dp/mprj_lc_c) ) )
622 
623  return
624  end subroutine mprj_lambertconformal_xy2lonlat
625 
626  !-----------------------------------------------------------------------------
628  subroutine mprj_lambertconformal_lonlat2xy( &
629  lon, &
630  lat, &
631  x, &
632  y )
633  implicit none
634 
635  real(RP), intent(in) :: lon ! [rad]
636  real(RP), intent(in) :: lat ! [rad]
637  real(RP), intent(out) :: x
638  real(RP), intent(out) :: y
639 
640  real(DP) :: dlon, latrot, dist
641  !---------------------------------------------------------------------------
642 
643  dlon = lon - mprj_basepoint_lon * d2r
644  if ( dlon > pi ) dlon = dlon - pi*2.0_dp
645  if ( dlon < -pi ) dlon = dlon + pi*2.0_dp
646 
647  latrot = 0.5_dp*pi - lat
648 
649  dist = mprj_lc_fact * radius * tan(0.5_dp*latrot)**mprj_lc_c
650 
651  x = mprj_pole_x + dist * sin(mprj_lc_c*dlon)
652  y = mprj_pole_y - mprj_hemisphere * dist * cos(mprj_lc_c*dlon)
653 
654  return
655  end subroutine mprj_lambertconformal_lonlat2xy
656 
657  !-----------------------------------------------------------------------------
659  subroutine mprj_lambertconformal_mapfactor( &
660  lat, &
661  m1, &
662  m2 )
663  implicit none
664 
665  real(RP), intent(in) :: lat(ia,ja) ! [rad]
666  real(RP), intent(out) :: m1 (ia,ja)
667  real(RP), intent(out) :: m2 (ia,ja)
668 
669  real(DP) :: latrot
670  integer :: i, j
671  !---------------------------------------------------------------------------
672 
673  do j = 1, ja
674  do i = 1, ia
675  latrot = 0.5_dp*pi - mprj_hemisphere * lat(i,j)
676 
677  m1(i,j) = mprj_lc_fact / sin(latrot) * mprj_lc_c * tan(0.5_dp*latrot)**mprj_lc_c
678  m2(i,j) = m1(i,j)
679  enddo
680  enddo
681 
682  return
683  end subroutine mprj_lambertconformal_mapfactor
684 
685  !-----------------------------------------------------------------------------
686  subroutine mprj_lambertconformal_rotcoef_0d( &
687  rotc, &
688  lon, &
689  lat )
690  implicit none
691 
692  real(RP), intent(out) :: rotc(2)
693  real(RP), intent(in) :: lon ! [rad]
694  real(RP), intent(in) :: lat ! [rad]
695 
696  real(DP) :: dlon
697  real(DP) :: alpha
698  !---------------------------------------------------------------------------
699 
700  dlon = lon - mprj_basepoint_lon * d2r
701  if( dlon > pi ) dlon = dlon - pi*2.0_dp
702  if( dlon < -pi ) dlon = dlon + pi*2.0_dp
703  alpha = - mprj_lc_c * dlon * mprj_hemisphere
704  rotc(1) = cos( alpha )
705  rotc(2) = sin( alpha )
706 
707  return
708  end subroutine mprj_lambertconformal_rotcoef_0d
709 
710  !-----------------------------------------------------------------------------
711  subroutine mprj_lambertconformal_rotcoef_2d( &
712  rotc, &
713  lon, &
714  lat )
715  implicit none
716 
717  real(RP), intent(out) :: rotc(ia,ja,2)
718  real(RP), intent(in) :: lon (ia,ja) ! [rad]
719  real(RP), intent(in) :: lat (ia,ja) ! [rad]
720 
721  real(DP) :: dlon
722  real(DP) :: alpha
723 
724  integer :: i, j
725  !---------------------------------------------------------------------------
726 
727  do j = 1, ja
728  do i = 1, ia
729  dlon = lon(i,j) - mprj_basepoint_lon * d2r
730  if( dlon > pi ) dlon = dlon - pi*2.0_dp
731  if( dlon < -pi ) dlon = dlon + pi*2.0_dp
732  alpha = - mprj_lc_c * dlon * mprj_hemisphere
733  rotc(i,j,1) = cos( alpha )
734  rotc(i,j,2) = sin( alpha )
735  enddo
736  enddo
737 
738  return
739  end subroutine mprj_lambertconformal_rotcoef_2d
740 
741  !-----------------------------------------------------------------------------
743  subroutine mprj_polarstereographic_setup
744  implicit none
745 
746  real(DP) :: lat0
747  real(DP) :: dlon, latrot, dist
748  !---------------------------------------------------------------------------
749 
750  ! check hemisphere: 1=north, -1=south
751  mprj_hemisphere = sign(1.0_dp,mprj_ps_lat)
752 
753  lat0 = mprj_hemisphere * mprj_ps_lat * d2r
754 
755  ! pre-calc factor
756  mprj_ps_fact = 1.0_dp + sin(lat0)
757 
758  ! calc (x,y) at pole point
759  dlon = 0.0_dp
760 
761  latrot = 0.5_dp*pi - mprj_hemisphere * mprj_basepoint_lat * d2r
762 
763  dist = mprj_ps_fact * radius * tan(0.5_dp*latrot)
764 
765  mprj_pole_x = mprj_basepoint_x - dist * sin(dlon)
766  mprj_pole_y = mprj_basepoint_y + mprj_hemisphere * dist * cos(dlon)
767 
768  if( io_l ) write(io_fid_log,*)
769  if( io_l ) write(io_fid_log,*) '+++ MPRJ_PS_lat1 :', mprj_ps_lat
770  if( io_l ) write(io_fid_log,*) '+++ MPRJ_hemisphere:', mprj_hemisphere
771  if( io_l ) write(io_fid_log,*) '+++ MPRJ_PS_fact :', mprj_ps_fact
772  if( io_l ) write(io_fid_log,*) '+++ MPRJ_pole_x :', mprj_pole_x
773  if( io_l ) write(io_fid_log,*) '+++ MPRJ_pole_y :', mprj_pole_y
774 
775  return
776  end subroutine mprj_polarstereographic_setup
777 
778  !-----------------------------------------------------------------------------
780  subroutine mprj_polarstereographic_xy2lonlat( &
781  x, &
782  y, &
783  lon, &
784  lat )
785  implicit none
786 
787  real(RP), intent(in) :: x
788  real(RP), intent(in) :: y
789  real(RP), intent(out) :: lon ! [rad]
790  real(RP), intent(out) :: lat ! [rad]
791 
792  real(DP) :: xx, yy, dist
793  !---------------------------------------------------------------------------
794 
795  xx = ( x - mprj_pole_x ) / radius / mprj_ps_fact
796  yy = -mprj_hemisphere * ( y - mprj_pole_y ) / radius / mprj_ps_fact
797 
798  dist = sqrt( xx*xx + yy*yy )
799 
800  lon = mprj_basepoint_lon * d2r + atan2(xx,yy)
801  lon = mod( lon+2.0_dp*pi, 2.0_dp*pi )
802  lat = mprj_hemisphere * ( 0.5_dp*pi - 2.0_dp*atan(dist) )
803 
804  return
805  end subroutine mprj_polarstereographic_xy2lonlat
806 
807  !-----------------------------------------------------------------------------
809  subroutine mprj_polarstereographic_lonlat2xy( &
810  lon, &
811  lat, &
812  x, &
813  y )
814  implicit none
815 
816  real(RP), intent(in) :: lon ! [rad]
817  real(RP), intent(in) :: lat ! [rad]
818  real(RP), intent(out) :: x
819  real(RP), intent(out) :: y
820 
821  real(DP) :: dlon, latrot, dist
822  !---------------------------------------------------------------------------
823 
824  dlon = lon - mprj_basepoint_lon * d2r
825 
826  latrot = 0.5_dp*pi - lat
827 
828  dist = mprj_ps_fact * radius * tan(0.5_dp*latrot)
829 
830  x = mprj_pole_x + dist * sin(dlon)
831  y = mprj_pole_y - mprj_hemisphere * dist * cos(dlon)
832 
833  return
834  end subroutine mprj_polarstereographic_lonlat2xy
835 
836  !-----------------------------------------------------------------------------
838  subroutine mprj_polarstereographic_mapfactor( &
839  lat, &
840  m1, &
841  m2 )
842  implicit none
843 
844  real(RP), intent(in) :: lat(ia,ja) ! [rad]
845  real(RP), intent(out) :: m1 (ia,ja)
846  real(RP), intent(out) :: m2 (ia,ja)
847 
848  integer :: i, j
849  !---------------------------------------------------------------------------
850 
851  do j = 1, ja
852  do i = 1, ia
853  m1(i,j) = mprj_ps_fact / ( 1.0_dp + sin(mprj_hemisphere*lat(i,j)) )
854  m2(i,j) = m1(i,j)
855  enddo
856  enddo
857 
858  return
859  end subroutine mprj_polarstereographic_mapfactor
860 
861  !-----------------------------------------------------------------------------
862  subroutine mprj_polarstereographic_rotcoef_0d( &
863  rotc, &
864  lon, &
865  lat )
866  implicit none
867 
868  real(RP), intent(out) :: rotc(2)
869  real(RP), intent(in) :: lon ! [rad]
870  real(RP), intent(in) :: lat ! [rad]
871 
872  real(DP) :: dlon
873  real(DP) :: alpha
874  !---------------------------------------------------------------------------
875 
876  dlon = lon - mprj_basepoint_lon * d2r
877  if( dlon > pi ) dlon = dlon - pi*2.0_dp
878  if( dlon < -pi ) dlon = dlon + pi*2.0_dp
879  alpha = - dlon * mprj_hemisphere
880  rotc(1) = cos( alpha )
881  rotc(2) = sin( alpha )
882 
883  return
884  end subroutine mprj_polarstereographic_rotcoef_0d
885 
886  !-----------------------------------------------------------------------------
887  subroutine mprj_polarstereographic_rotcoef_2d( &
888  rotc, &
889  lon, &
890  lat )
891  implicit none
892 
893  real(RP), intent(out) :: rotc(ia,ja,2)
894  real(RP), intent(in) :: lon (ia,ja) ! [rad]
895  real(RP), intent(in) :: lat (ia,ja) ! [rad]
896 
897  real(DP) :: dlon
898  real(DP) :: alpha
899 
900  integer :: i, j
901  !---------------------------------------------------------------------------
902 
903  do j = 1, ja
904  do i = 1, ia
905  dlon = lon(i,j) - mprj_basepoint_lon * d2r
906  if( dlon > pi ) dlon = dlon - pi*2.0_dp
907  if( dlon < -pi ) dlon = dlon + pi*2.0_dp
908  alpha = - dlon * mprj_hemisphere
909  rotc(i,j,1) = cos( alpha )
910  rotc(i,j,2) = sin( alpha )
911  enddo
912  enddo
913 
914  return
915  end subroutine mprj_polarstereographic_rotcoef_2d
916 
917  !-----------------------------------------------------------------------------
919  subroutine mprj_mercator_setup
920  implicit none
921 
922  real(DP) :: lat0
923  real(DP) :: latrot, dist
924  !---------------------------------------------------------------------------
925 
926  lat0 = mprj_m_lat * d2r
927 
928  ! pre-calc factor
929  mprj_m_fact = cos(lat0)
930 
931  ! calc (x,y) at (lon,lat) = (base,0)
932  latrot = 0.5_dp*pi - mprj_basepoint_lat * d2r
933 
934  dist = 1.0_dp / tan(0.5_dp*latrot)
935 
936  mprj_eq_x = mprj_basepoint_x
937  mprj_eq_y = mprj_basepoint_y - radius * mprj_m_fact * log(dist)
938 
939  if( io_l ) write(io_fid_log,*)
940  if( io_l ) write(io_fid_log,*) '+++ MPRJ_M_lat :', mprj_m_lat
941  if( io_l ) write(io_fid_log,*) '+++ MPRJ_M_fact:', mprj_m_fact
942  if( io_l ) write(io_fid_log,*) '+++ MPRJ_eq_x :', mprj_eq_x
943  if( io_l ) write(io_fid_log,*) '+++ MPRJ_eq_y :', mprj_eq_y
944 
945  return
946  end subroutine mprj_mercator_setup
947 
948  !-----------------------------------------------------------------------------
950  subroutine mprj_mercator_xy2lonlat( &
951  x, &
952  y, &
953  lon, &
954  lat )
955  implicit none
956 
957  real(RP), intent(in) :: x
958  real(RP), intent(in) :: y
959  real(RP), intent(out) :: lon ! [rad]
960  real(RP), intent(out) :: lat ! [rad]
961 
962  real(DP) :: xx, yy
963  !---------------------------------------------------------------------------
964 
965  xx = ( x - mprj_eq_x ) / radius / mprj_m_fact
966  yy = ( y - mprj_eq_y ) / radius / mprj_m_fact
967 
968  lon = xx + mprj_basepoint_lon * d2r
969  lon = mod( lon+2.0_dp*pi, 2.0_dp*pi )
970  lat = 0.5_dp*pi - 2.0_dp*atan( 1.0_dp/exp(yy) )
971 
972  return
973  end subroutine mprj_mercator_xy2lonlat
974 
975  !-----------------------------------------------------------------------------
977  subroutine mprj_mercator_lonlat2xy( &
978  lon, &
979  lat, &
980  x, &
981  y )
982  implicit none
983 
984  real(RP), intent(in) :: lon ! [rad]
985  real(RP), intent(in) :: lat ! [rad]
986  real(RP), intent(out) :: x
987  real(RP), intent(out) :: y
988 
989  real(DP) :: dlon, latrot, dist
990  !---------------------------------------------------------------------------
991 
992  dlon = lon - mprj_basepoint_lon * d2r
993 
994  latrot = 0.5_dp*pi - lat
995 
996  dist = 1.0_dp / tan(0.5_dp*latrot)
997 
998  x = mprj_eq_x + radius * mprj_m_fact * dlon
999  y = mprj_eq_y + radius * mprj_m_fact * log(dist)
1000 
1001  return
1002  end subroutine mprj_mercator_lonlat2xy
1003 
1004  !-----------------------------------------------------------------------------
1006  subroutine mprj_mercator_mapfactor( &
1007  lat, &
1008  m1, &
1009  m2 )
1010  implicit none
1011 
1012  real(RP), intent(in) :: lat(ia,ja) ! [rad]
1013  real(RP), intent(out) :: m1 (ia,ja)
1014  real(RP), intent(out) :: m2 (ia,ja)
1015  integer :: i, j
1016  !---------------------------------------------------------------------------
1017 
1018  do j = 1, ja
1019  do i = 1, ia
1020  m1(i,j) = mprj_m_fact / cos(real(lat(i,j),kind=dp))
1021  m2(i,j) = m1(i,j)
1022  enddo
1023  enddo
1024 
1025  return
1026  end subroutine mprj_mercator_mapfactor
1027 
1028  !-----------------------------------------------------------------------------
1029  subroutine mprj_mercator_rotcoef_0d( &
1030  rotc )
1031  implicit none
1032 
1033  real(RP), intent(out) :: rotc(2)
1034  !---------------------------------------------------------------------------
1035 
1036  rotc(1) = 1.0_rp
1037  rotc(2) = 0.0_rp
1038 
1039  return
1040  end subroutine mprj_mercator_rotcoef_0d
1041 
1042  !-----------------------------------------------------------------------------
1043  subroutine mprj_mercator_rotcoef_2d( &
1044  rotc )
1045  implicit none
1046 
1047  real(RP), intent(out) :: rotc(ia,ja,2)
1048  !---------------------------------------------------------------------------
1049 
1050  rotc(:,:,1) = 1.0_rp
1051  rotc(:,:,2) = 0.0_rp
1052 
1053  return
1054  end subroutine mprj_mercator_rotcoef_2d
1055 
1056  !-----------------------------------------------------------------------------
1058  subroutine mprj_equidistantcylindrical_setup
1059  implicit none
1060 
1061  real(DP) :: lat0
1062  !---------------------------------------------------------------------------
1063 
1064  lat0 = mprj_ec_lat * d2r
1065 
1066  ! pre-calc factor
1067  mprj_ec_fact = cos(lat0)
1068 
1069  mprj_eq_x = mprj_basepoint_x
1070  mprj_eq_y = mprj_basepoint_y - radius * mprj_basepoint_lat * d2r
1071 
1072  if( io_l ) write(io_fid_log,*)
1073  if( io_l ) write(io_fid_log,*) '+++ MPRJ_EC_lat :', mprj_ec_lat
1074  if( io_l ) write(io_fid_log,*) '+++ MPRJ_EC_fact:', mprj_ec_fact
1075  if( io_l ) write(io_fid_log,*) '+++ MPRJ_eq_x :', mprj_eq_x
1076  if( io_l ) write(io_fid_log,*) '+++ MPRJ_eq_y :', mprj_eq_y
1077 
1078  return
1079  end subroutine mprj_equidistantcylindrical_setup
1080 
1081  !-----------------------------------------------------------------------------
1083  subroutine mprj_equidistantcylindrical_xy2lonlat( &
1084  x, &
1085  y, &
1086  lon, &
1087  lat )
1088  implicit none
1089 
1090  real(RP), intent(in) :: x
1091  real(RP), intent(in) :: y
1092  real(RP), intent(out) :: lon ! [rad]
1093  real(RP), intent(out) :: lat ! [rad]
1094 
1095  real(DP) :: xx, yy
1096  !---------------------------------------------------------------------------
1097 
1098  xx = ( x - mprj_eq_x ) / radius / mprj_ec_fact
1099  yy = ( y - mprj_eq_y ) / radius
1100 
1101  lon = xx + mprj_basepoint_lon * d2r
1102  lon = mod( lon+2.0_dp*pi, 2.0_dp*pi )
1103  lat = yy
1104 
1105  return
1106  end subroutine mprj_equidistantcylindrical_xy2lonlat
1107 
1108  !-----------------------------------------------------------------------------
1110  subroutine mprj_equidistantcylindrical_lonlat2xy( &
1111  lon, &
1112  lat, &
1113  x, &
1114  y )
1115  implicit none
1116 
1117  real(RP), intent(in) :: lon ! [rad]
1118  real(RP), intent(in) :: lat ! [rad]
1119  real(RP), intent(out) :: x
1120  real(RP), intent(out) :: y
1121 
1122  real(DP) :: dlon
1123  !---------------------------------------------------------------------------
1124 
1125  dlon = lon - mprj_basepoint_lon * d2r
1126 
1127  x = mprj_eq_x + radius * mprj_ec_fact * dlon
1128  y = mprj_eq_y + radius * lat
1129 
1130  return
1131  end subroutine mprj_equidistantcylindrical_lonlat2xy
1132 
1133  !-----------------------------------------------------------------------------
1135  subroutine mprj_equidistantcylindrical_mapfactor( &
1136  lat, &
1137  m1, &
1138  m2 )
1139  implicit none
1140 
1141  real(RP), intent(in) :: lat(ia,ja) ! [rad]
1142  real(RP), intent(out) :: m1 (ia,ja)
1143  real(RP), intent(out) :: m2 (ia,ja)
1144  integer :: i, j
1145  !---------------------------------------------------------------------------
1146 
1147  do j = 1, ja
1148  do i = 1, ia
1149  m1(i,j) = mprj_ec_fact / cos(real(lat(i,j),kind=dp))
1150  m2(i,j) = 1.0_rp
1151  enddo
1152  enddo
1153 
1154  return
1155  end subroutine mprj_equidistantcylindrical_mapfactor
1156 
1157  !-----------------------------------------------------------------------------
1158  subroutine mprj_equidistantcylindrical_rotcoef_0d( &
1159  rotc )
1160  implicit none
1161 
1162  real(RP), intent(out) :: rotc(2)
1163  !---------------------------------------------------------------------------
1164 
1165  rotc(1) = 1.0_rp
1166  rotc(2) = 0.0_rp
1167 
1168  return
1169  end subroutine mprj_equidistantcylindrical_rotcoef_0d
1170 
1171  !-----------------------------------------------------------------------------
1172  subroutine mprj_equidistantcylindrical_rotcoef_2d( &
1173  rotc )
1174  implicit none
1175 
1176  real(RP), intent(out) :: rotc(ia,ja,2)
1177  !---------------------------------------------------------------------------
1178 
1179  rotc(:,:,1) = 1.0_rp
1180  rotc(:,:,2) = 0.0_rp
1181 
1182  return
1183  end subroutine mprj_equidistantcylindrical_rotcoef_2d
1184 
1185 end module scale_mapproj
scale_process::prc_mpistop
subroutine, public prc_mpistop
Abort MPI.
Definition: scale_process.F90:234
scale_mapproj::mprj_basepoint_lat
real(rp), public mprj_basepoint_lat
Definition: scale_mapproj.F90:49
scale_stdio::io_l
logical, public io_l
output log or not? (this process)
Definition: scale_stdio.F90:59
scale_mapproj::mprj_xy2lonlat
subroutine, public mprj_xy2lonlat(x, y, lon, lat)
(x,y) -> (lon,lat)
Definition: scale_mapproj.F90:236
scale_const::const_radius
real(rp), public const_radius
radius of the planet [m]
Definition: scale_const.F90:46
scale_stdio
module STDIO
Definition: scale_stdio.F90:12
scale_const::const_d2r
real(rp), public const_d2r
degree to radian
Definition: scale_const.F90:35
scale_mapproj::mprj_basepoint_lon
real(rp), public mprj_basepoint_lon
Definition: scale_mapproj.F90:48
scale_const::const_undef
real(rp), public const_undef
Definition: scale_const.F90:43
scale_mapproj::mprj_rotcoef_2d
subroutine mprj_rotcoef_2d(rotc, lon, lat)
u(lat,lon) = cos u(x,y) - sin v(x,y) v(lat,lon) = sin u(x,y) + cos v(x,y)
Definition: scale_mapproj.F90:377
scale_grid_index
module grid index
Definition: scale_grid_index.F90:14
scale_grid_index::ia
integer, public ia
of x whole cells (local, with HALO)
Definition: scale_grid_index.F90:55
dc_types::dp
integer, parameter, public dp
Definition: dc_types.f90:27
scale_process
module PROCESS
Definition: scale_process.F90:14
scale_mapproj::mprj_lonlat2xy
subroutine, public mprj_lonlat2xy(lon, lat, x, y)
(lon,lat) -> (x,y)
Definition: scale_mapproj.F90:272
dc_log::log
subroutine, public log(type, message)
Definition: dc_log.f90:133
scale_const
module CONSTANT
Definition: scale_const.F90:14
scale_mapproj::mprj_mapfactor
subroutine, public mprj_mapfactor(lat, m1, m2)
(x,y) -> (lon,lat)
Definition: scale_mapproj.F90:307
scale_prof
module profiler
Definition: scale_prof.F90:10
scale_precision
module PRECISION
Definition: scale_precision.F90:13
scale_const::const_pi
real(rp), public const_pi
pi
Definition: scale_const.F90:34
scale_mapproj::mprj_rotcoef_0d
subroutine mprj_rotcoef_0d(rotc, lon, lat)
u(lat,lon) = cos u(x,y) - sin v(x,y) v(lat,lon) = sin u(x,y) + cos v(x,y)
Definition: scale_mapproj.F90:342
scale_stdio::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_stdio.F90:55
scale_stdio::io_fid_log
integer, public io_fid_log
Log file ID.
Definition: scale_stdio.F90:56
scale_mapproj::mprj_setup
subroutine, public mprj_setup(DOMAIN_CENTER_X, DOMAIN_CENTER_Y)
Setup.
Definition: scale_mapproj.F90:137
scale_mapproj
module Map projection
Definition: scale_mapproj.F90:14
scale_grid_index::ja
integer, public ja
of y whole cells (local, with HALO)
Definition: scale_grid_index.F90:56
Generated by   doxygen 1.8.13