Changeset 153

Timestamp:

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

Author:

dubos

Message:

Schmidt transform

Location:

codes/icosagcm/trunk/src

Files:

• geometry.f90 (modified) (7 diffs)
• spherical_geom.f90 (modified) (1 diff)

codes/icosagcm/trunk/src/geometry.f90

@@ -55 +55 @@
   INTEGER, PARAMETER :: ne_ldown=1
   INTEGER, PARAMETER :: ne_rdown=-1
-  
+
+  REAL(rstd),private :: schmidt_factor, schmidt_lon, schmidt_lat
+
 CONTAINS
 
@@ -115 +117 @@
 !$OMP BARRIER    
   END SUBROUTINE swap_geometry
-  
+
+  SUBROUTINE update_circumcenters 
+    USE domain_mod
+    USE dimensions
+    USE spherical_geom_mod
+    USE vector
+    USE transfert_mod
+
+    IMPLICIT NONE
+    REAL(rstd) :: x1(3),x2(3)
+    REAL(rstd) :: vect(3,6)
+    REAL(rstd) :: centr(3)
+
+    INTEGER :: ind,i,j,n,k
+
+    CALL transfert_request(geom%xyz_i,req_i1)
+
+    DO ind=1,ndomain
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          n=(j-1)*iim+i
+          DO k=0,5
+            x1(:) = xyz_i(n+t_pos(k+1),:)
+            x2(:) = xyz_i(n+t_pos(MOD(k+1,6)+1),:)
+            if (norm(x1-x2)<1e-16) x2(:) = xyz_i(n+t_pos(MOD(k+2,6)+1),:)
+            CALL circumcenter(xyz_i(n,:), x1, x2, xyz_v(n+z_pos(k+1),:))
+          ENDDO
+        ENDDO
+      ENDDO
+    ENDDO
+
+  END SUBROUTINE update_circumcenters
+
   SUBROUTINE optimize_geometry
   USE metric
@@ -149 +185 @@
     ENDDO
     
-    CALL transfert_request(geom%xyz_i,req_i1)
-    
-    DO ind=1,ndomain
-      CALL swap_dimensions(ind)
-      CALL swap_geometry(ind)
-      DO j=jj_begin,jj_end
-        DO i=ii_begin,ii_end
-          n=(j-1)*iim+i
-          DO k=0,5
-            x1(:) = xyz_i(n+t_pos(k+1),:)
-            x2(:) = xyz_i(n+t_pos(MOD(k+1,6)+1),:)
-            if (norm(x1-x2)<1e-16) x2(:) = xyz_i(n+t_pos(MOD(k+2,6)+1),:)
-            CALL circumcenter(xyz_i(n,:), x1, x2, xyz_v(n+z_pos(k+1),:))
-          ENDDO
-        ENDDO
-      ENDDO
-    ENDDO
-    
+    CALL update_circumcenters
+
     DO ind=1,ndomain
       d=>domain(ind)
@@ -221 +241 @@
         PRINT *,"it = ",it,"  diff centroid circumcenter ",sum
       ENDIF 
+
+     CALL update_circumcenters
+
+    ENDDO
+    
+  END SUBROUTINE optimize_geometry
+
+  SUBROUTINE schmidt_transform(xyz,cc)
+    ! Based on formula (12) from Guo & Drake, JCP 2005
+    USE spherical_geom_mod
+    IMPLICIT NONE
+    REAL(rstd),INTENT(INOUT) :: xyz(3)
+    REAL(rstd), INTENT(IN) :: cc  ! stretching factor>0
+    REAL(rstd) :: lat,lon,mu
+    CALL xyz2lonlat_relative(xyz,schmidt_lon,schmidt_lat, lon,lat)
+    mu = sin(lat)
+    mu = ((cc-1)+mu*(cc+1)) / ((cc+1)+mu*(cc-1))
+    lat = asin(mu)
+    CALL lonlat2xyz_relative(lon,lat,schmidt_lon, schmidt_lat, xyz)
+  END SUBROUTINE schmidt_transform
       
-     CALL transfert_request(geom%xyz_i,req_i1)
-
-    DO ind=1,ndomain
-      CALL swap_dimensions(ind)
-      CALL swap_geometry(ind)
-      DO j=jj_begin,jj_end
-        DO i=ii_begin,ii_end
-          n=(j-1)*iim+i
-          DO k=0,5
-            x1(:) = xyz_i(n+t_pos(k+1),:)
-            x2(:) = xyz_i(n+t_pos(MOD(k+1,6)+1),:)
-            if (norm(x1-x2)<1e-16) x2(:) = xyz_i(n+t_pos(MOD(k+2,6)+1),:)
-            CALL circumcenter(xyz_i(n,:), x1, x2, xyz_v(n+z_pos(k+1),:))
-          ENDDO
-        ENDDO
-      ENDDO
-    ENDDO
-
-      
-    ENDDO
-    
-  END SUBROUTINE optimize_geometry
-  
-    
   SUBROUTINE set_geometry
   USE metric
@@ -253 +269 @@
   USE dimensions
   USE transfert_mod
+  USE ioipsl
   IMPLICIT NONE
 
@@ -259 +276 @@
     REAL(rstd) :: vect(3,6)
     REAL(rstd) :: centr(3)
-    REAL(rstd) :: vet(3),vep(3)
+    REAL(rstd) :: vet(3),vep(3), vertex(3)
     INTEGER :: ind,i,j,k,n
     TYPE(t_domain),POINTER :: d
@@ -270 +287 @@
       
     CALL optimize_geometry
-    
+
+    ! Schmidt transform
+    schmidt_factor = 1.
+    CALL getin('schmidt_factor', schmidt_factor)
+
+    schmidt_lon = 0.
+    CALL getin('schmidt_lon', schmidt_lon)
+    schmidt_lon = schmidt_lon * pi/180.
+
+    schmidt_lat = 45.
+    CALL getin('schmidt_lat', schmidt_lat)
+    schmidt_lat = schmidt_lat * pi/180.
+
+    DO ind=1,ndomain
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      DO j=jj_begin-1,jj_end+1
+        DO i=ii_begin-1,ii_end+1
+          n=(j-1)*iim+i
+          vertex(:)=xyz_i(n,:)
+          CALL schmidt_transform(vertex, schmidt_factor**2.)
+          xyz_i(n,:)=vertex(:)
+       END DO
+     END DO
+   END DO
+
+    ! at this point the mesh is defined on the unit sphere by xyz_i
+    ! all other mesh quantities are now deduced from xyz_i
+
+    CALL update_circumcenters
+
+!    DO ind=1,ndomain
+!      d=>domain(ind)
+!      CALL swap_dimensions(ind)
+!      CALL swap_geometry(ind)
+!      DO j=1,jjm
+!        DO i=1,iim
+!          n=(j-1)*iim+i
+!          d%xyz(:,i,j)=xyz_i(n,:)
+!        ENDDO
+!      ENDDO
+!    ENDDO
+!
+!    CALL compute_boundary
+ 
     DO ind=1,ndomain
       d=>domain(ind)

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
+