doc-doxygen/html/sqdis_8f_source.html

00001       FUNCTION sqdis (px1, py1, px2, py2)
00002 C****
00003 C               ******************************
00004 C               * OASIS FUNCTION  -  LEVEL T *
00005 C               * --------------     ------- *
00006 C               ******************************
00007 C
00008 C**** *sqdis*  - Arithmetic function
00009 C
00010 C     Purpose:
00011 C     -------
00012 C     Calculate the distance squared between 2 points on a spheric grid
00013 C
00014 C**   Interface:
00015 C     ---------
00016 C       *zs =*  *sqdis (px1, py1, px2, py2)*
00017 C
00018 C     Input:
00019 C     -----
00020 C                px1    : longitude of first point
00021 C                py1    : latitude of first point
00022 C                px2    : longitude of second point
00023 C                py2    : latitude of second point
00024 C
00025 C     NB: the coordinates must be in degrees
00026 C
00027 C     Output:
00028 C     ------
00029 C     None
00030 C
00031 C     Workspace:
00032 C     ---------
00033 C     None
00034 C
00035 C     Externals:
00036 C     ---------
00037 C     None
00038 C
00039 C     Reference:
00040 C     ---------
00041 C     See OASIS manual (1995)
00042 C
00043 C     History:
00044 C     -------
00045 C       Version   Programmer     Date      Description
00046 C       -------   ----------     ----      -----------
00047 C       2.0       L. Terray      95/09/01  created
00048 C
00049 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00050 C
00051 C* ---------------------------- Include files ---------------------------
00052 C
00053       USE mod_unit
00054 C
00055 C* ---------------------------- Function declarations -------------------
00056 C
00057       REAL (kind=ip_realwp_p) sqdis, fast
00058       fast(x,y,z) = cos(x)*cos(y)*z + sin(x)*sin(y)
00059 C
00060 C* ---------------------------- Poema verses ----------------------------
00061 C
00062 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00063 C
00064 C*    1. Calculate distance
00065 C        ------------------
00066 C
00067 C* Degree to radians --> zcon = 2. pi / 360.
00068 C
00069       zconv = 1.74532925199432957692e-2
00070       zfi1 = px1 * zconv
00071       zth1 = py1 * zconv
00072       zfi2 = px2 * zconv
00073       zth2 = py2 * zconv
00074 C
00075 C* Scalar product
00076 C
00077       ztp1 = fast(zfi1, zfi2, 1.)
00078       zsca = fast(zth1, zth2, ztp1)
00079 C
00080 C* Angular distance
00081 C
00082       IF (zsca .GT. 1.) THEN
00083           WRITE(UNIT = nulan,FMT = *)
00084      $        ' Scalar product is greater than 1 '
00085           WRITE(UNIT = nulan,FMT = *)
00086      $        ' for point 1 -->  longitude = ', px1
00087           WRITE(UNIT = nulan,FMT = *)
00088      $        '             -->  latitude  = ', py1
00089           WRITE(UNIT = nulan,FMT = *)
00090      $        ' and point 2 -->  longitude = ', px2
00091           WRITE(UNIT = nulan,FMT = *)
00092      $        '             -->  latitude  = ', py2
00093           WRITE(UNIT = nulan,FMT = *)
00094      $        ' Scalar product  zsca = ', zsca
00095           zsca =1.
00096       ENDIF
00097       IF (zsca .LT. -1.) THEN
00098           WRITE(UNIT = nulan,FMT = *)
00099      $        ' Scalar product is less than -1 '
00100           WRITE(UNIT = nulan,FMT = *)
00101      $        ' for point 1 -->  longitude = ', px1
00102           WRITE(UNIT = nulan,FMT = *)
00103      $        '             -->  latitude  = ', py1
00104           WRITE(UNIT = nulan,FMT = *)
00105      $        ' and point 2 -->  longitude = ', px2
00106           WRITE(UNIT = nulan,FMT = *)
00107      $        '             -->  latitude  = ', py2
00108           WRITE(UNIT = nulan,FMT = *)
00109      $        ' Scalar product  zsca = ', zsca
00110           zsca = -1.
00111       ENDIF
00112 C
00113 C* Get the distance
00114 C
00115       zalp = acos(zsca)
00116       sqdis = zalp * zalp
00117 C
00118 C* End of function
00119 C
00120       RETURN
00121       END