SCALE-RM
Functions/Subroutines | Variables
scale_vector Module Reference

module vector More...

Functions/Subroutines

subroutine, public vectr_xyz2latlon (x, y, z, lat, lon)
 
subroutine, public vectr_latlon2xyz (lat, lon, x, y, z, radius)
 
subroutine, public vectr_cross (nv, a, b, c, d)
 exterior product of vector a->b and c->d More...
 
subroutine, public vectr_dot (l, a, b, c, d)
 interior product of vector a->b and c->d More...
 
subroutine, public vectr_abs (l, a)
 length of vector o->a More...
 
subroutine, public vectr_angle (angle, a, b, c)
 calc angle between two vector(b->a,b->c) More...
 
subroutine, public vectr_intersec (ifcross, p, a, b, c, d)
 judge intersection of two vector More...
 
subroutine, public vectr_anticlockwise (vertex, nvert)
 bubble sort anticlockwise by angle More...
 
real(rp) function, public vectr_triangle (a, b, c, polygon_type, radius)
 calc triangle area More...
 
real(rp) function, public vectr_triangle_plane (a, b, c)
 calc triangle area on plane More...
 
subroutine, public vectr_rotation (a, angle, iaxis)
 Apply rotation matrix. More...
 
subroutine, public vectr_distance (r, lon1, lat1, lon2, lat2, dist)
 Get horizontal distance on the sphere. More...
 

Variables

integer, parameter, public i_xaxis = 1
 
integer, parameter, public i_yaxis = 2
 
integer, parameter, public i_zaxis = 3
 

Detailed Description

module vector

Description
module for 3D vector on the sphere
Author
NICAM developers, Team SCALE

Function/Subroutine Documentation

◆ vectr_xyz2latlon()

subroutine, public scale_vector::vectr_xyz2latlon ( real(rp), intent(in)  x,
real(rp), intent(in)  y,
real(rp), intent(in)  z,
real(rp), intent(out)  lat,
real(rp), intent(out)  lon 
)

Definition at line 63 of file scale_vector.F90.

References scale_const::const_eps.

63  use scale_const, only: &
64  eps => const_eps
65  implicit none
66 
67  real(RP), intent(in) :: x
68  real(RP), intent(in) :: y
69  real(RP), intent(in) :: z
70  real(RP), intent(out) :: lat
71  real(RP), intent(out) :: lon
72 
73  real(RP) :: length, length_h
74  !---------------------------------------------------------------------------
75 
76  length = sqrt( x*x + y*y + z*z )
77 
78  if ( length < eps ) then ! 3D vector length is
79  lat = 0.0_rp
80  lon = 0.0_rp
81  return
82  endif
83 
84  if ( z / length >= 1.0_rp ) then ! vector is parallele to z axis.
85  lat = asin(1.0_rp)
86  lon = 0.0_rp
87  return
88  elseif( z / length <= -1.0_rp ) then ! vector is parallele to z axis.
89  lat = asin(-1.0_rp)
90  lon = 0.0_rp
91  return
92  else
93  lat = asin( z / length )
94  endif
95 
96  length_h = sqrt( x*x + y*y )
97 
98  if ( length_h < eps ) then
99  lon = 0.0_rp
100  return
101  endif
102 
103  if ( x / length_h >= 1.0_rp ) then
104  lon = acos(1.0_rp)
105  elseif( x / length_h <= -1.0_rp ) then
106  lon = acos(-1.0_rp)
107  else
108  lon = acos( x / length_h )
109  endif
110 
111  if( y < 0.0_rp ) lon = -lon
112 
113  return
scale_const
module CONSTANT
Definition: scale_const.F90:14

◆ vectr_latlon2xyz()

subroutine, public scale_vector::vectr_latlon2xyz ( real(rp), intent(in)  lat,
real(rp), intent(in)  lon,
real(rp), intent(out)  x,
real(rp), intent(out)  y,
real(rp), intent(out)  z,
real(rp), intent(in)  radius 
)

Definition at line 124 of file scale_vector.F90.

124  implicit none
125 
126  real(RP), intent(in) :: lat
127  real(RP), intent(in) :: lon
128  real(RP), intent(out) :: x
129  real(RP), intent(out) :: y
130  real(RP), intent(out) :: z
131  real(RP), intent(in) :: radius
132  !---------------------------------------------------------------------------
133 
134  x = radius * cos(lat) * cos(lon)
135  y = radius * cos(lat) * sin(lon)
136  z = radius * sin(lat)
137 
138  return

◆ vectr_cross()

subroutine, public scale_vector::vectr_cross ( real(rp), dimension(3), intent(out)  nv,
real(rp), dimension(3), intent(in)  a,
real(rp), dimension(3), intent(in)  b,
real(rp), dimension(3), intent(in)  c,
real(rp), dimension(3), intent(in)  d 
)

exterior product of vector a->b and c->d

Definition at line 144 of file scale_vector.F90.

Referenced by vectr_angle(), vectr_anticlockwise(), vectr_intersec(), and vectr_triangle().

144  implicit none
145 
146  real(RP), intent(out) :: nv(3) ! normal vector
147  real(RP), intent(in) :: a(3), b(3), c(3), d(3) ! x,y,z(cartesian)
148  !---------------------------------------------------------------------------
149 
150  nv(1) = ( b(2)-a(2) ) * ( d(3)-c(3) ) &
151  - ( b(3)-a(3) ) * ( d(2)-c(2) )
152  nv(2) = ( b(3)-a(3) ) * ( d(1)-c(1) ) &
153  - ( b(1)-a(1) ) * ( d(3)-c(3) )
154  nv(3) = ( b(1)-a(1) ) * ( d(2)-c(2) ) &
155  - ( b(2)-a(2) ) * ( d(1)-c(1) )
156 
157  return
Here is the caller graph for this function:

◆ vectr_dot()

subroutine, public scale_vector::vectr_dot ( real(rp), intent(out)  l,
real(rp), dimension(3), intent(in)  a,
real(rp), dimension(3), intent(in)  b,
real(rp), dimension(3), intent(in)  c,
real(rp), dimension(3), intent(in)  d 
)

interior product of vector a->b and c->d

Definition at line 163 of file scale_vector.F90.

Referenced by vectr_angle(), vectr_anticlockwise(), vectr_intersec(), and vectr_triangle_plane().

163  implicit none
164 
165  real(RP), intent(out) :: l
166  real(RP), intent(in) :: a(3), b(3), c(3), d(3) ! x,y,z(cartesian)
167  !---------------------------------------------------------------------------
168  ! if a=c=zero-vector and b=d, result is abs|a|^2
169 
170  l = ( b(1)-a(1) ) * ( d(1)-c(1) ) &
171  + ( b(2)-a(2) ) * ( d(2)-c(2) ) &
172  + ( b(3)-a(3) ) * ( d(3)-c(3) )
173 
174  return
Here is the caller graph for this function:

◆ vectr_abs()

subroutine, public scale_vector::vectr_abs ( real(rp), intent(out)  l,
real(rp), dimension(3), intent(in)  a 
)

length of vector o->a

Definition at line 180 of file scale_vector.F90.

Referenced by vectr_angle(), vectr_intersec(), and vectr_triangle().

180  implicit none
181 
182  real(RP), intent(out) :: l
183  real(RP), intent(in) :: a(3) ! x,y,z(cartesian)
184  !---------------------------------------------------------------------------
185 
186  l = a(1)*a(1) + a(2)*a(2) + a(3)*a(3)
187  l = sqrt(l)
188 
189  return
Here is the caller graph for this function:

◆ vectr_angle()

subroutine, public scale_vector::vectr_angle ( real(rp), intent(out)  angle,
real(rp), dimension(3), intent(in)  a,
real(rp), dimension(3), intent(in)  b,
real(rp), dimension(3), intent(in)  c 
)

calc angle between two vector(b->a,b->c)

Definition at line 195 of file scale_vector.F90.

References vectr_abs(), vectr_cross(), and vectr_dot().

Referenced by vectr_anticlockwise(), vectr_intersec(), and vectr_triangle().

195  implicit none
196 
197  real(RP), intent(out) :: angle
198  real(RP), intent(in) :: a(3), b(3), c(3)
199 
200  real(RP) :: nv(3), nvlenS, nvlenC
201  !---------------------------------------------------------------------
202 
203  call vectr_dot ( nvlenc, b, a, b, c )
204  call vectr_cross( nv(:), b, a, b, c )
205  call vectr_abs ( nvlens, nv(:) )
206  angle = atan2( nvlens, nvlenc )
207 
208  return
Here is the call graph for this function:
Here is the caller graph for this function:

◆ vectr_intersec()

subroutine, public scale_vector::vectr_intersec ( logical, intent(out)  ifcross,
real(rp), dimension(3), intent(out)  p,
real(rp), dimension(3), intent(in)  a,
real(rp), dimension(3), intent(in)  b,
real(rp), dimension(3), intent(in)  c,
real(rp), dimension(3), intent(in)  d 
)

judge intersection of two vector

Definition at line 214 of file scale_vector.F90.

References scale_const::const_eps, vectr_abs(), vectr_angle(), vectr_cross(), and vectr_dot().

214  use scale_const, only: &
215  eps => const_eps
216  implicit none
217 
218  logical, intent(out) :: ifcross
219  ! .true. : line a->b and c->d intersect
220  ! .false.: line a->b and c->d do not intersect and p = (0,0)
221  real(RP), intent(out) :: p(3) ! intersection point
222  real(RP), intent(in) :: a(3), b(3), c(3), d(3)
223 
224  real(RP), parameter :: o(3) = 0.0_rp
225 
226  real(RP) :: oaob(3), ocod(3), cdab(3)
227  real(RP) :: ip, length
228  real(RP) :: angle_aop, angle_pob, angle_aob
229  real(RP) :: angle_cop, angle_pod, angle_cod
230  !---------------------------------------------------------------------
231 
232  call vectr_cross( oaob, o, a, o, b )
233  call vectr_cross( ocod, o, c, o, d )
234  call vectr_cross( cdab, o, ocod, o, oaob )
235 
236  call vectr_abs ( length, cdab )
237  call vectr_dot ( ip, o, cdab, o, a )
238 
239  p(:) = cdab(:) / sign(length,ip)
240 ! if( IO_L ) write(IO_FID_LOG,*), "p:", p(:)
241 
242  call vectr_angle( angle_aop, a, o, p )
243  call vectr_angle( angle_pob, p, o, b )
244  call vectr_angle( angle_aob, a, o, b )
245 ! if( IO_L ) write(IO_FID_LOG,*), "angle a-p-b:", angle_aop, angle_pob, angle_aob
246 
247  call vectr_angle( angle_cop, c, o, p )
248  call vectr_angle( angle_pod, p, o, d )
249  call vectr_angle( angle_cod, c, o, d )
250 ! if( IO_L ) write(IO_FID_LOG,*), "angle c-p-d:", angle_cop, angle_pod, angle_cod
251 
252 ! if( IO_L ) write(IO_FID_LOG,*), "judge:", angle_aob-(angle_aop+angle_pob), angle_cod-(angle_cop+angle_pod)
253 
254  ! --- judge intersection
255  if ( abs(angle_aob-(angle_aop+angle_pob)) < eps &
256  .AND. abs(angle_cod-(angle_cop+angle_pod)) < eps &
257  .AND. abs(angle_aop) > eps &
258  .AND. abs(angle_pob) > eps &
259  .AND. abs(angle_cop) > eps &
260  .AND. abs(angle_pod) > eps ) then
261  ifcross = .true.
262  else
263  ifcross = .false.
264  p(:) = 0.0_rp
265  endif
266 
267  return
scale_const
module CONSTANT
Definition: scale_const.F90:14
Here is the call graph for this function:

◆ vectr_anticlockwise()

subroutine, public scale_vector::vectr_anticlockwise ( real(rp), dimension(nvert,3), intent(inout)  vertex,
integer, intent(in)  nvert 
)

bubble sort anticlockwise by angle

Definition at line 273 of file scale_vector.F90.

References scale_const::const_eps, vectr_angle(), vectr_cross(), and vectr_dot().

273  use scale_const, only: &
274  eps => const_eps
275  implicit none
276 
277  integer, intent(in) :: nvert
278  real(RP), intent(inout) :: vertex(nvert,3)
279 
280  real(RP), parameter :: o(3) = 0.0_rp
281 
282  real(RP) :: v1(3), v2(3), v3(3)
283  real(RP) :: xp(3), ip
284  real(RP) :: angle1, angle2
285 
286  integer :: i, j
287  !---------------------------------------------------------------------
288 
289  do j = 2 , nvert-1
290  do i = j+1, nvert
291  v1(:) = vertex(1,:)
292  v2(:) = vertex(j,:)
293  v3(:) = vertex(i,:)
294 
295  call vectr_cross( xp(:), v1(:), v2(:), v1(:), v3(:) )
296  call vectr_dot ( ip, o(:), v1(:), o(:), xp(:) )
297 
298  if ( ip < -eps ) then ! right hand : exchange
299 ! if( IO_L ) write(IO_FID_LOG,*) 'exchange by ip', i, '<->',j
300  vertex(i,:) = v2(:)
301  vertex(j,:) = v3(:)
302  endif
303 
304  enddo
305  enddo
306 
307  v1(:) = vertex(1,:)
308  v2(:) = vertex(2,:)
309  v3(:) = vertex(3,:)
310  ! if 1->2->3 is on the line
311  call vectr_cross( xp(:), v1(:), v2(:), v1(:), v3(:) )
312  call vectr_dot ( ip, o(:), v1(:), o(:), xp(:) )
313  call vectr_angle( angle1, v1(:), o, v2(:) )
314  call vectr_angle( angle2, v1(:), o, v3(:) )
315 ! if( IO_L ) write(IO_FID_LOG,*) ip, angle1, angle2, abs(angle1)-abs(angle2)
316 
317  if ( abs(ip) < eps & ! on the same line
318  .AND. abs(angle2)-abs(angle1) < 0.0_rp ) then ! which is far?
319 ! if( IO_L ) write(IO_FID_LOG,*) 'exchange by angle', 2, '<->', 3
320  vertex(2,:) = v3(:)
321  vertex(3,:) = v2(:)
322  endif
323 
324  v2(:) = vertex(nvert ,:)
325  v3(:) = vertex(nvert-1,:)
326  ! if 1->nvert->nvert-1 is on the line
327  call vectr_cross( xp(:), v1(:), v2(:), v1(:), v3(:) )
328  call vectr_dot ( ip, o(:), v1(:), o(:), xp(:) )
329  call vectr_angle( angle1, v1(:), o, v2(:) )
330  call vectr_angle( angle2, v1(:), o, v3(:) )
331 ! if( IO_L ) write(IO_FID_LOG,*) ip, angle1, angle2, abs(angle1)-abs(angle2)
332 
333  if ( abs(ip) < eps & ! on the same line
334  .AND. abs(angle2)-abs(angle1) < 0.0_rp ) then ! which is far?
335 ! if( IO_L ) write(IO_FID_LOG,*) 'exchange by angle', nvert, '<->', nvert-1
336  vertex(nvert, :) = v3(:)
337  vertex(nvert-1,:) = v2(:)
338  endif
339 
340  return
scale_const
module CONSTANT
Definition: scale_const.F90:14
Here is the call graph for this function:

◆ vectr_triangle()

real(rp) function, public scale_vector::vectr_triangle ( real(rp), dimension(3), intent(in)  a,
real(rp), dimension(3), intent(in)  b,
real(rp), dimension(3), intent(in)  c,
character(len=*), intent(in)  polygon_type,
real(rp), intent(in)  radius 
)

calc triangle area

Returns
area

Definition at line 351 of file scale_vector.F90.

References scale_const::const_eps, scale_const::const_pi, vectr_abs(), vectr_angle(), and vectr_cross().

351  use scale_const, only: &
352  pi => const_pi, &
353  eps => const_eps
354  implicit none
355 
356  real(RP), intent(in) :: a(3), b(3), c(3)
357  character(len=*), intent(in) :: polygon_type
358  real(RP), intent(in) :: radius
359  real(RP) :: area
360 
361  real(RP), parameter :: o(3) = 0.0_rp
362 
363  ! ON_PLANE
364  real(RP) :: abc(3)
365  real(RP) :: prd, r
366 
367  ! ON_SPHERE
368  real(RP) :: angle(3)
369  real(RP) :: oaob(3), oaoc(3)
370  real(RP) :: oboc(3), oboa(3)
371  real(RP) :: ocoa(3), ocob(3)
372  real(RP) :: abab, acac
373  real(RP) :: bcbc, baba
374  real(RP) :: caca, cbcb
375  !---------------------------------------------------------------------------
376 
377  area = 0.0_rp
378 
379  if ( polygon_type == 'ON_PLANE' ) then ! Note : On a plane, area = | ourter product of two vectors |.
380 
381  call vectr_cross( abc(:), a(:), b(:), a(:), c(:) )
382  call vectr_abs( prd, abc(:) )
383  call vectr_abs( r , a(:) )
384 
385  prd = 0.5_rp * prd !! triangle area
386  if ( r < eps ) then
387  print *, "zero length?", a(:)
388  else
389  r = 1.0_rp / r !! 1 / length
390  endif
391 
392  area = prd * r*r * radius*radius
393 
394  elseif( polygon_type == 'ON_SPHERE' ) then ! On a unit sphere, area = sum of angles - pi.
395 
396  ! angle 1
397  call vectr_cross( oaob(:), o(:), a(:), o(:), b(:) )
398  call vectr_cross( oaoc(:), o(:), a(:), o(:), c(:) )
399  call vectr_abs( abab, oaob(:) )
400  call vectr_abs( acac, oaoc(:) )
401 
402  if ( abab < eps .OR. acac < eps ) then
403  !write(*,'(A,3(ES20.10))') "zero length abab or acac:", abab, acac
404  return
405  endif
406 
407  call vectr_angle( angle(1), oaob(:), o(:), oaoc(:) )
408 
409  ! angle 2
410  call vectr_cross( oboc(:), o(:), b(:), o(:), c(:) )
411  oboa(:) = -oaob(:)
412  call vectr_abs( bcbc, oboc(:) )
413  baba = abab
414 
415  if ( bcbc < eps .OR. baba < eps ) then
416  !write(*,'(A,3(ES20.10))') "zero length bcbc or baba:", bcbc, baba
417  return
418  endif
419 
420  call vectr_angle( angle(2), oboc(:), o(:), oboa(:) )
421 
422  ! angle 3
423  ocoa(:) = -oaoc(:)
424  ocob(:) = -oboc(:)
425  caca = acac
426  cbcb = bcbc
427 
428  if ( caca < eps .OR. cbcb < eps ) then
429  !write(*,'(A,3(ES20.10))') "zero length caca or cbcb:", caca, cbcb
430  return
431  endif
432 
433  call vectr_angle( angle(3), ocoa(:), o(:), ocob(:) )
434 
435  ! calc area
436  area = ( angle(1)+angle(2)+angle(3) - pi ) * radius*radius
437 
438  endif
439 
440  return
scale_const
module CONSTANT
Definition: scale_const.F90:14
Here is the call graph for this function:

◆ vectr_triangle_plane()

real(rp) function, public scale_vector::vectr_triangle_plane ( real(rp), dimension(3), intent(in)  a,
real(rp), dimension(3), intent(in)  b,
real(rp), dimension(3), intent(in)  c 
)

calc triangle area on plane

Returns
area

Definition at line 449 of file scale_vector.F90.

References vectr_dot().

449  implicit none
450 
451  real(RP), intent(in) :: a(3), b(3), c(3)
452  real(RP) :: area
453  !
454  real(RP) :: len_ab, len_ac, prd
455  !---------------------------------------------------------------------------
456 
457  call vectr_dot( len_ab, a, b, a, b )
458  call vectr_dot( len_ac, a, c, a, c )
459  call vectr_dot( prd , a, b, a, c )
460 
461  area = 0.5_rp * sqrt( len_ab * len_ac - prd * prd )
462 
Here is the call graph for this function:

◆ vectr_rotation()

subroutine, public scale_vector::vectr_rotation ( real(rp), dimension(3), intent(inout)  a,
real(rp), intent(in)  angle,
integer, intent(in)  iaxis 
)

Apply rotation matrix.

Definition at line 471 of file scale_vector.F90.

References i_xaxis, i_yaxis, and i_zaxis.

471  implicit none
472 
473  real(RP), intent(inout) :: a(3)
474  real(RP), intent(in) :: angle
475  integer, intent(in) :: iaxis
476 
477  real(RP) :: m(3,3), b(3)
478  !---------------------------------------------------------------------------
479 
480  if ( iaxis == i_xaxis ) then
481  m(1,1) = 1.0_rp
482  m(1,2) = 0.0_rp
483  m(1,3) = 0.0_rp
484 
485  m(2,1) = 0.0_rp
486  m(2,2) = cos(angle)
487  m(2,3) = sin(angle)
488 
489  m(3,1) = 0.0_rp
490  m(3,2) = -sin(angle)
491  m(3,3) = cos(angle)
492  elseif( iaxis == i_yaxis ) then
493  m(1,1) = cos(angle)
494  m(1,2) = 0.0_rp
495  m(1,3) = -sin(angle)
496 
497  m(2,1) = 0.0_rp
498  m(2,2) = 1.0_rp
499  m(2,3) = 0.0_rp
500 
501  m(3,1) = sin(angle)
502  m(3,2) = 0.0_rp
503  m(3,3) = cos(angle)
504  elseif( iaxis == i_zaxis ) then
505  m(1,1) = cos(angle)
506  m(1,2) = sin(angle)
507  m(1,3) = 0.0_rp
508 
509  m(2,1) = -sin(angle)
510  m(2,2) = cos(angle)
511  m(2,3) = 0.0_rp
512 
513  m(3,1) = 0.0_rp
514  m(3,2) = 0.0_rp
515  m(3,3) = 1.0_rp
516  else
517  return
518  endif
519 
520  b(1) = m(1,1) * a(1) + m(1,2) * a(2) + m(1,3) * a(3)
521  b(2) = m(2,1) * a(1) + m(2,2) * a(2) + m(2,3) * a(3)
522  b(3) = m(3,1) * a(1) + m(3,2) * a(2) + m(3,3) * a(3)
523 
524  a(:) = b(:)
525 
526  return

◆ vectr_distance()

subroutine, public scale_vector::vectr_distance ( real(rp), intent(in)  r,
real(rp), intent(in)  lon1,
real(rp), intent(in)  lat1,
real(rp), intent(in)  lon2,
real(rp), intent(in)  lat2,
real(rp), intent(out)  dist 
)

Get horizontal distance on the sphere.

Definition at line 538 of file scale_vector.F90.

Referenced by scale_gridtrans::gtrans_rotcoef().

538  implicit none
539 
540  real(RP), intent(in) :: r ! radius in meter
541  real(RP), intent(in) :: lon1, lat1 ! in radian
542  real(RP), intent(in) :: lon2, lat2 ! in radian
543  real(RP), intent(out) :: dist ! distance of the two points in meter
544 
545  real(RP) :: gmm, gno_x, gno_y
546  !-----------------------------------------------------------------------
547 
548  gmm = sin(lat1) * sin(lat2) &
549  + cos(lat1) * cos(lat2) * cos(lon2-lon1)
550 
551  gno_x = gmm * ( cos(lat2) * sin(lon2-lon1) )
552  gno_y = gmm * ( cos(lat1) * sin(lat2) &
553  - sin(lat1) * cos(lat2) * cos(lon2-lon1) )
554 
555  dist = r * atan2( sqrt(gno_x*gno_x+gno_y*gno_y), gmm )
556 
557  return
Here is the caller graph for this function:

Variable Documentation

◆ i_xaxis

integer, parameter, public scale_vector::i_xaxis = 1

Definition at line 42 of file scale_vector.F90.

Referenced by vectr_rotation().

42  integer, public, parameter :: I_Xaxis = 1

◆ i_yaxis

integer, parameter, public scale_vector::i_yaxis = 2

Definition at line 43 of file scale_vector.F90.

Referenced by vectr_rotation().

43  integer, public, parameter :: I_Yaxis = 2

◆ i_zaxis

integer, parameter, public scale_vector::i_zaxis = 3

Definition at line 44 of file scale_vector.F90.

Referenced by vectr_rotation().

44  integer, public, parameter :: I_Zaxis = 3
Generated by   doxygen 1.8.14