Changeset 153 for codes/icosagcm/trunk/src/spherical_geom.f90

Timestamp:

05/17/13 18:21:50 (11 years ago)

Author:

dubos

Message:

Schmidt transform

File:

• codes/icosagcm/trunk/src/spherical_geom.f90 (modified) (1 diff)

codes/icosagcm/trunk/src/spherical_geom.f90

@@ -34 +34 @@
   
   lat=asin(xyzn(3))
-  coslat=cos(lat)
-  IF (abs(coslat)<1e-15) THEN
-    lon=0.
-  ELSE
-    lon=acos(MAX(MIN((xyzn(1)/coslat),1.),-1.))
-    IF (xyzn(2)/coslat<0) lon=-lon
-  ENDIF
-    
+  lon=atan2(xyzn(2),xyzn(1))
 END SUBROUTINE xyz2lonlat
+
+! lat/lon with respect to a displaced pole (rotated basis) :
+!  ( cos(lon0)*sin(lat0), sin(lon0)*sin(lat0), -cos(lat0)) 
+!  (-sin(lon0),           cos(lon0),           0)
+!  ( cos(lon0)*cos(lat0), sin(lon0)*cos(lat0), sin(lat0))
+
+SUBROUTINE lonlat2xyz_relative(lon,lat,lon0,lat0, xyz)
+IMPLICIT NONE
+  REAL(rstd),INTENT(IN) :: lon0, lat0, lon,lat
+  REAL(rstd),INTENT(OUT) :: xyz(3)
+  REAL(rstd) :: xx,yy,zz
+  xx = cos(lon)*cos(lat)
+  yy = sin(lon)*cos(lat)
+  zz = sin(lat)
+  xyz(1) = cos(lon0)*(sin(lat0)*xx+cos(lat0)*zz)-sin(lon0)*yy
+  xyz(2) = sin(lon0)*(sin(lat0)*yy+cos(lat0)*zz)+cos(lon0)*yy
+  xyz(3) = sin(lat0)*zz-cos(lat0)*xx
+END SUBROUTINE lonlat2xyz_relative
+
+SUBROUTINE xyz2lonlat_relative(xyz,lon0,lat0, lon,lat)
+IMPLICIT NONE
+  REAL(rstd),INTENT(IN) :: xyz(3), lon0, lat0
+  REAL(rstd),INTENT(OUT) :: lon,lat
+  REAL(rstd) :: xx,yy,zz
+  xx = sin(lat0)*(xyz(1)*cos(lon0)+xyz(2)*sin(lon0))-cos(lat0)*xyz(3)
+  yy = xyz(2)*cos(lon0)-xyz(1)*sin(lon0)
+  zz = cos(lat0)*(xyz(1)*cos(lon0)+xyz(2)*sin(lon0))+sin(lat0)*xyz(3)
+  lon = atan2(yy,xx)
+  lat = asin(zz)
+END SUBROUTINE xyz2lonlat_relative
+