Changeset 155

Timestamp:

05/20/13 16:35:23 (11 years ago)

Author:

dubos

Message:

Schmidt transform in metric.f90

Location:

codes/icosagcm/trunk/src

Files:

• geometry.f90 (modified) (3 diffs)
• metric.f90 (modified) (2 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
@@ -248 +246 @@
   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
-      
   SUBROUTINE set_geometry
   USE metric
@@ -287 +271 @@
       
     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

codes/icosagcm/trunk/src/metric.f90

@@ -52 +52 @@
   INTEGER,PARAMETER :: vdown=5
   INTEGER,PARAMETER :: vrdown=6
-
   
 CONTAINS
-  
+
+  SUBROUTINE remap_schmidt
+    USE spherical_geom_mod
+    USE ioipsl
+    IMPLICIT NONE
+    INTEGER :: nf,i,j
+    REAL(rstd) :: schmidt_factor, schmidt_lon, schmidt_lat
+
+    ! Schmidt transform parameters
+    schmidt_factor = 1.
+    CALL getin('schmidt_factor', schmidt_factor)
+    schmidt_factor =  schmidt_factor**2.
+    
+    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 nf=1,nb_face
+       DO j=1,jjm_glo
+          DO i=1,iim_glo
+             CALL schmidt_transform(vertex_glo(i,j,nf)%xyz, schmidt_factor, schmidt_lon, schmidt_lat)
+          END DO
+       END DO
+    END DO
+  END SUBROUTINE remap_schmidt
+
   SUBROUTINE allocate_metric
   IMPLICIT NONE
@@ -830 +858 @@
 !    CALL compute_face
     CALL compute_face_projection
+    CALL remap_schmidt
+
     CALL set_index
     CALL set_cell

codes/icosagcm/trunk/src/spherical_geom.f90

@@ -67 +67 @@
 END SUBROUTINE xyz2lonlat_relative
 
-
+SUBROUTINE schmidt_transform(xyz,cc, lon0, lat0)
+  ! Based on formula (12) from Guo & Drake, JCP 2005
+  IMPLICIT NONE
+  REAL(rstd),INTENT(INOUT) :: xyz(3)
+  REAL(rstd), INTENT(IN) :: cc, lon0, lat0  ! stretching factor>0, lon/lat of zoomed area
+  REAL(rstd) :: lat,lon,mu
+  CALL xyz2lonlat_relative(xyz,lon0,lat0, lon,lat)
+  mu = sin(lat)
+  mu = ((cc-1)+mu*(cc+1)) / ((cc+1)+mu*(cc-1))
+  lat = asin(mu)
+  CALL lonlat2xyz_relative(lon,lat, lon0,lat0, xyz)
+END SUBROUTINE schmidt_transform
 
 SUBROUTINE dist_cart(A,B,d)