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