Changeset 15 for codes/icosagcm/trunk/src/geometry.f90

Timestamp:

07/09/12 15:23:38 (12 years ago)

Author:

ymipsl

Message:

Update on 3D dynamic

YM

File:

• codes/icosagcm/trunk/src/geometry.f90 (modified) (8 diffs)

codes/icosagcm/trunk/src/geometry.f90

@@ -4 +4 @@
   TYPE t_geometry
     TYPE(t_field),POINTER :: xyz_i(:)
+    TYPE(t_field),POINTER :: centroid(:)
     TYPE(t_field),POINTER :: xyz_e(:)
     TYPE(t_field),POINTER :: xyz_v(:)
+    TYPE(t_field),POINTER :: ep_e(:)
+    TYPE(t_field),POINTER :: et_e(:)
+    TYPE(t_field),POINTER :: elon_i(:)
+    TYPE(t_field),POINTER :: elat_i(:)
+    TYPE(t_field),POINTER :: elon_e(:)
+    TYPE(t_field),POINTER :: elat_e(:)
     TYPE(t_field),POINTER :: Ai(:)
     TYPE(t_field),POINTER :: Av(:)
@@ -20 +27 @@
   TYPE(t_geometry),TARGET :: geom
   
-  REAL(rstd),POINTER :: Ai(:)
-  REAL(rstd),POINTER :: xyz_i(:,:)
-  REAL(rstd),POINTER :: xyz_e(:,:)
-  REAL(rstd),POINTER :: xyz_v(:,:)
-  REAL(rstd),POINTER :: Av(:)
-  REAL(rstd),POINTER :: de(:)
-  REAL(rstd),POINTER :: le(:)
-  REAL(rstd),POINTER :: Riv(:,:)
-  REAL(rstd),POINTER :: Riv2(:,:)
-  INTEGER,POINTER :: ne(:,:)
-  REAL(rstd),POINTER :: Wee(:,:,:)
-  REAL(rstd),POINTER :: bi(:)
-  REAL(rstd),POINTER :: fv(:)
+  REAL(rstd),POINTER :: Ai(:)          ! area of a cell
+  REAL(rstd),POINTER :: xyz_i(:,:)     ! coordinate of the center of the cell (voronoi)
+  REAL(rstd),POINTER :: centroid(:,:)  ! coordinate of the centroid of the cell
+  REAL(rstd),POINTER :: xyz_e(:,:)     ! coordinate of a wind point on the cell on a edge
+  REAL(rstd),POINTER :: ep_e(:,:)      ! perpendicular unit vector of a edge (outsider)
+  REAL(rstd),POINTER :: et_e(:,:)      ! tangeantial unit vector of a edge
+  REAL(rstd),POINTER :: elon_i(:,:)    ! unit longitude vector on the center 
+  REAL(rstd),POINTER :: elat_i(:,:)    ! unit latitude vector on the center 
+  REAL(rstd),POINTER :: elon_e(:,:)    ! unit longitude vector on a wind point 
+  REAL(rstd),POINTER :: elat_e(:,:)    ! unit latitude vector on a wind point
+  REAL(rstd),POINTER :: xyz_v(:,:)     ! coordinate of a vertex (center of the dual mesh)
+  REAL(rstd),POINTER :: Av(:)          ! area of dual mesk cell
+  REAL(rstd),POINTER :: de(:)          ! distance from a neighbour == lenght of an edge of the dual mesh
+  REAL(rstd),POINTER :: le(:)          ! lenght of a edge
+  REAL(rstd),POINTER :: Riv(:,:)       ! weight
+  REAL(rstd),POINTER :: Riv2(:,:)      ! weight
+  INTEGER,POINTER    :: ne(:,:)        ! convention for the way on the normal wind on an edge 
+  REAL(rstd),POINTER :: Wee(:,:,:)     ! weight
+  REAL(rstd),POINTER :: bi(:)          ! orographie
+  REAL(rstd),POINTER :: fv(:)          ! coriolis (evaluted on a vertex)
      
 
@@ -44 +58 @@
     CALL allocate_field(geom%Ai,field_t,type_real)
     CALL allocate_field(geom%xyz_i,field_t,type_real,3)
+    CALL allocate_field(geom%centroid,field_t,type_real,3)
     CALL allocate_field(geom%xyz_e,field_u,type_real,3)
+    CALL allocate_field(geom%ep_e,field_u,type_real,3)
+    CALL allocate_field(geom%et_e,field_u,type_real,3)
+    CALL allocate_field(geom%elon_i,field_t,type_real,3)
+    CALL allocate_field(geom%elat_i,field_t,type_real,3)
+    CALL allocate_field(geom%elon_e,field_u,type_real,3)
+    CALL allocate_field(geom%elat_e,field_u,type_real,3)
     CALL allocate_field(geom%xyz_v,field_z,type_real,3)
     CALL allocate_field(geom%de,field_u,type_real)
@@ -66 +87 @@
     Ai=geom%Ai(ind)
     xyz_i=geom%xyz_i(ind)
+    centroid=geom%centroid(ind)
     xyz_e=geom%xyz_e(ind)
+    ep_e=geom%ep_e(ind)
+    et_e=geom%et_e(ind)
+    elon_i=geom%elon_i(ind)
+    elat_i=geom%elat_i(ind)
+    elon_e=geom%elon_e(ind)
+    elat_e=geom%elat_e(ind)
     xyz_v=geom%xyz_v(ind)
     de=geom%de(ind)
@@ -79 +107 @@
     
   END SUBROUTINE swap_geometry
-    
-  SUBROUTINE set_geometry
+  
+  SUBROUTINE optimize_geometry
   USE metric
   USE spherical_geom_mod
   USE domain_mod
   USE dimensions
+  USE transfert_mod
+  USE vector
+  IMPLICIT NONE
+    INTEGER,PARAMETER :: nb_it=3000
+    TYPE(t_domain),POINTER :: d
+    INTEGER :: ind,it,i,j,n,k
+    REAL(rstd) :: x1(3),x2(3)
+    REAL(rstd) :: vect(3,6)
+    REAL(rstd) :: centr(3)
+    REAL(rstd) :: sum
+    LOGICAL    :: check
+    
+    DO ind=1,ndomain
+      d=>domain(ind)
+      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
+          xyz_i(n,:)=d%xyz(:,i,j) 
+        ENDDO
+      ENDDO
+    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
+    
+    DO ind=1,ndomain
+      d=>domain(ind)
+      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_v(n+z_pos(k+1),:)
+            x2(:) = d%vertex(:,k,i,j) 
+            IF (norm(x1-x2)>1e-10) THEN
+              PRINT*,"vertex diff ",ind,i,j,k
+              PRINT*,x1
+              PRINT*,x2
+            ENDIF
+          ENDDO
+        ENDDO
+      ENDDO
+    ENDDO
+    
+    
+    DO it=1,nb_it
+      IF (MOD(it,100)==0) THEN 
+        check=.TRUE.
+      ELSE
+        check=.FALSE.
+      ENDIF
+      
+      sum=0
+      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
+            vect(:,1)=xyz_v(n+z_rup,:)
+            vect(:,2)=xyz_v(n+z_up,:)
+            vect(:,3)=xyz_v(n+z_lup,:)
+            vect(:,4)=xyz_v(n+z_ldown,:)
+            vect(:,5)=xyz_v(n+z_down,:)
+            vect(:,6)=xyz_v(n+z_rdown,:)
+            CALL compute_centroid(vect,6,centr)
+            IF (check) THEN
+              sum=MAX(sum,norm(xyz_i(n,:)-centr(:)))
+            ENDIF
+            xyz_i(n,:)=centr(:)
+          ENDDO
+        ENDDO
+!        i=ii_begin ; j=jj_begin ; n=(j-1)*iim+i ; xyz_i(n,:)=domain(ind)%xyz(:,i,j) 
+!        i=ii_begin ; j=jj_end   ; n=(j-1)*iim+i ; xyz_i(n,:)=domain(ind)%xyz(:,i,j)
+!        i=ii_end   ; j=jj_begin ; n=(j-1)*iim+i ;   xyz_i(n,:)=domain(ind)%xyz(:,i,j)
+!        PRINT *,"Pb ?? : " 
+!        PRINT *, xyz_i(n,:), domain(ind)%xyz(:,i,j), norm(xyz_i(n,:)- domain(ind)%xyz(:,i,j))
+!        i=ii_end   ; j=jj_end   ; n=(j-1)*iim+i ; xyz_i(n,:)=domain(ind)%xyz(:,i,j)
+        
+      ENDDO
+
+      IF (check) THEN
+!        sum=sum/(ndomain*ii_nb*jj_nb)
+        PRINT *,"it = ",it,"  diff centroid circumcenter ",sum
+      ENDIF 
+      
+     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
+  USE vector
+  USE spherical_geom_mod
+  USE domain_mod
+  USE dimensions
+  USE transfert_mod
   IMPLICIT NONE
 
     REAL(rstd) :: surf(6)
     REAL(rstd) :: surf_v(6)
+    REAL(rstd) :: vect(3,6)
+    REAL(rstd) :: centr(6)
+    REAL(rstd) :: vet(3),vep(3)
     INTEGER :: ind,i,j,k,n
     TYPE(t_domain),POINTER :: d
     REAL(rstd) ::  S
     REAL(rstd) :: w(6)
-    REAl(rstd) :: lon,lat
+    REAL(rstd) :: lon,lat
     INTEGER :: ii_glo,jj_glo
     REAL(rstd) :: S1,S2
+          
       
+    CALL optimize_geometry
+    
     DO ind=1,ndomain
       d=>domain(ind)
@@ -104 +273 @@
         DO i=ii_begin-1,ii_end+1
           n=(j-1)*iim+i
-        
-          xyz_i(n,:)=d%xyz(:,i,j) 
-          xyz_v(n+z_up,:)=d%vertex(:,vup-1,i,j) 
-          xyz_v(n+z_down,:)=d%vertex(:,vdown-1,i,j) 
-          
+
+          DO k=0,5
+            ne(n,k+1)=d%ne(k,i,j)
+          ENDDO
+                  
+!          xyz_i(n,:)=d%xyz(:,i,j) 
+!          xyz_v(n+z_up,:)=d%vertex(:,vup-1,i,j) 
+!          xyz_v(n+z_down,:)=d%vertex(:,vdown-1,i,j) 
+
+          vect(:,1)=xyz_v(n+z_rup,:)
+          vect(:,2)=xyz_v(n+z_up,:)
+          vect(:,3)=xyz_v(n+z_lup,:)
+          vect(:,4)=xyz_v(n+z_ldown,:)
+          vect(:,5)=xyz_v(n+z_down,:)
+          vect(:,6)=xyz_v(n+z_rdown,:)
+          CALL compute_centroid(vect,6,centr)
+          centroid(n,:)=centr(:)
+
+               
           CALL xyz2lonlat(xyz_v(n+z_up,:),lon,lat)
           fv(n+z_up)=2*sin(lat)*omega
@@ -116 +299 @@
           bi(n)=0. 
         
-          CALL dist_cart(d%xyz(:,i,j),d%neighbour(:,right-1,i,j),de(n+u_right))
-          CALL dist_cart(d%xyz(:,i,j),d%neighbour(:,lup-1,i,j),de(n+u_lup))
-          CALL dist_cart(d%xyz(:,i,j),d%neighbour(:,ldown-1,i,j),de(n+u_ldown))
-          
-          CALL div_arc_bis(d%xyz(:,i,j),d%neighbour(:,right-1,i,j),0.5,xyz_e(n+u_right,:))
-          CALL div_arc_bis(d%xyz(:,i,j),d%neighbour(:,lup-1,i,j),0.5,xyz_e(n+u_lup,:))
-          CALL div_arc_bis(d%xyz(:,i,j),d%neighbour(:,ldown-1,i,j),0.5,xyz_e(n+u_ldown,:))
-          
-          CALL dist_cart(d%vertex(:,vrdown-1,i,j),d%vertex(:,vrup-1,i,j),le(n+u_right))
-          CALL dist_cart(d%vertex(:,vup-1,i,j),d%vertex(:,vlup-1,i,j),le(n+u_lup))
-          CALL dist_cart(d%vertex(:,vldown-1,i,j),d%vertex(:,vdown-1,i,j),le(n+u_ldown))
-
+!          CALL dist_cart(d%xyz(:,i,j),d%neighbour(:,right-1,i,j),de(n+u_right))
+!          CALL dist_cart(d%xyz(:,i,j),d%neighbour(:,lup-1,i,j),de(n+u_lup))
+!          CALL dist_cart(d%xyz(:,i,j),d%neighbour(:,ldown-1,i,j),de(n+u_ldown))
+
+          CALL dist_cart(xyz_i(n,:),xyz_i(n+t_right,:),de(n+u_right))
+          CALL dist_cart(xyz_i(n,:),xyz_i(n+t_lup,:),de(n+u_lup))
+          CALL dist_cart(xyz_i(n,:),xyz_i(n+t_ldown,:),de(n+u_ldown))
+          
+!          CALL div_arc_bis(d%xyz(:,i,j),d%neighbour(:,right-1,i,j),0.5,xyz_e(n+u_right,:))
+!          CALL div_arc_bis(d%xyz(:,i,j),d%neighbour(:,lup-1,i,j),0.5,xyz_e(n+u_lup,:))
+!          CALL div_arc_bis(d%xyz(:,i,j),d%neighbour(:,ldown-1,i,j),0.5,xyz_e(n+u_ldown,:))
+
+          CALL div_arc_bis(xyz_i(n,:),xyz_i(n+t_right,:),0.5,xyz_e(n+u_right,:))
+          CALL div_arc_bis(xyz_i(n,:),xyz_i(n+t_lup,:),0.5,xyz_e(n+u_lup,:))
+          CALL div_arc_bis(xyz_i(n,:),xyz_i(n+t_ldown,:),0.5,xyz_e(n+u_ldown,:))
+         
+!          CALL dist_cart(d%vertex(:,vrdown-1,i,j),d%vertex(:,vrup-1,i,j),le(n+u_right))
+!          CALL dist_cart(d%vertex(:,vup-1,i,j),d%vertex(:,vlup-1,i,j),le(n+u_lup))
+!          CALL dist_cart(d%vertex(:,vldown-1,i,j),d%vertex(:,vdown-1,i,j),le(n+u_ldown))
+
+          CALL dist_cart(xyz_v(n+z_rdown,:), xyz_v(n+z_rup,:),le(n+u_right))
+          CALL dist_cart(xyz_v(n+z_up,:), xyz_v(n+z_lup,:),le(n+u_lup))
+          CALL dist_cart(xyz_v(n+z_ldown,:), xyz_v(n+z_down,:),le(n+u_ldown))
          
           Ai(n)=0
           DO k=0,5
-            CALL surf_triangle(d%xyz(:,i,j),d%neighbour(:,k,i,j),d%neighbour(:,MOD((k+1+6),6),i,j),surf_v(k+1))
-            CALL surf_triangle(d%xyz(:,i,j),d%vertex(:,MOD((k-1+6),6),i,j),d%vertex(:,k,i,j),surf(k+1))
-            ne(n,k+1)=d%ne(k,i,j)
+!            CALL surf_triangle(d%xyz(:,i,j),d%neighbour(:,k,i,j),d%neighbour(:,MOD((k+1+6),6),i,j),surf_v(k+1))
+!            CALL surf_triangle(d%xyz(:,i,j),d%vertex(:,MOD((k-1+6),6),i,j),d%vertex(:,k,i,j),surf(k+1))
+            CALL surf_triangle(xyz_i(n,:),xyz_i(n+t_pos(k+1),:),xyz_i(n+t_pos(MOD((k+1+6),6)+1),:),surf_v(k+1))
+            CALL surf_triangle(xyz_i(n,:),xyz_v(n+z_pos(MOD((k-1+6),6)+1),:),xyz_v(n+z_pos(k+1),:),surf(k+1))
             Ai(n)=Ai(n)+surf(k+1)
-          ENDDO
-
-
+            IF (i==ii_end .AND. j==jj_begin) THEN
+              IF (Ai(n)<1e20) THEN
+              ELSE
+                PRINT *,"PB !!",Ai(n),k,surf(k+1)
+                PRINT*,xyz_i(n,:),xyz_v(n+z_pos(MOD((k-1+6),6)+1),:),xyz_v(n+z_pos(k+1),:)
+              ENDIF
+            ENDIF
+          ENDDO
+
+
+          CALL cross_product2(xyz_v(n+z_rdown,:),xyz_v(n+z_rup,:),vep)
+          IF (norm(vep)>1e-30) THEN
+            vep(:)=vep(:)/norm(vep)
+            CALL cross_product2(vep,xyz_e(n+u_right,:),vet)
+            vet(:)=vet(:)/norm(vet)
+            ep_e(n+u_right,:)=vep(:)*ne(n,right)
+            et_e(n+u_right,:)=vet(:)*ne(n,right)
+          ENDIF
+
+          CALL cross_product2(xyz_v(n+z_up,:),xyz_v(n+z_lup,:),vep)
+          IF (norm(vep)>1e-30) THEN
+            vep(:)=vep(:)/norm(vep)
+            CALL cross_product2(vep,xyz_e(n+u_lup,:),vet)
+            vet(:)=vet(:)/norm(vet)
+            ep_e(n+u_lup,:)=vep(:)*ne(n,lup)
+            et_e(n+u_lup,:)=vet(:)*ne(n,lup)
+          ENDIF
+
+          CALL cross_product2(xyz_v(n+z_ldown,:),xyz_v(n+z_down,:),vep)
+          IF (norm(vep)>1e-30) THEN
+            vep(:)=vep(:)/norm(vep)
+            CALL cross_product2(vep,xyz_e(n+u_ldown,:),vet)
+            vet(:)=vet(:)/norm(vet)
+            ep_e(n+u_ldown,:)=vep(:)*ne(n,ldown)
+            et_e(n+u_ldown,:)=vet(:)*ne(n,ldown)
+          ENDIF
+
+          CALL xyz2lonlat(xyz_i(n,:),lon,lat)
+          elon_i(n,1) = -sin(lon)
+          elon_i(n,2) = cos(lon)
+          elon_i(n,3) = 0
+          elat_i(n,1) = -cos(lon)*sin(lat)
+          elat_i(n,2) = -sin(lon)*sin(lat)
+          elat_i(n,3) = cos(lat)
+
+          
+          CALL xyz2lonlat(xyz_e(n+u_right,:),lon,lat)
+          elon_e(n+u_right,1) = -sin(lon)
+          elon_e(n+u_right,2) = cos(lon)
+          elon_e(n+u_right,3) = 0
+          elat_e(n+u_right,1) = -cos(lon)*sin(lat)
+          elat_e(n+u_right,2) = -sin(lon)*sin(lat)
+          elat_e(n+u_right,3) = cos(lat)
+          
+          CALL xyz2lonlat(xyz_e(n+u_lup,:),lon,lat)
+          elon_e(n+u_lup,1) = -sin(lon)
+          elon_e(n+u_lup,2) = cos(lon)
+          elon_e(n+u_lup,3) = 0
+          elat_e(n+u_lup,1) = -cos(lon)*sin(lat)
+          elat_e(n+u_lup,2) = -sin(lon)*sin(lat)
+          elat_e(n+u_lup,3) = cos(lat)
+ 
+          CALL xyz2lonlat(xyz_e(n+u_ldown,:),lon,lat)
+          elon_e(n+u_ldown,1) = -sin(lon)
+          elon_e(n+u_ldown,2) = cos(lon)
+          elon_e(n+u_ldown,3) = 0
+          elat_e(n+u_ldown,1) = -cos(lon)*sin(lat)
+          elat_e(n+u_ldown,2) = -sin(lon)*sin(lat)
+          elat_e(n+u_ldown,3) = cos(lat)
+ 
+          
           DO k=0,5
-            CALL surf_triangle(d%xyz(:,i,j),d%vertex(:,k,i,j),d%neighbour(:,k,i,j),S1)
-            CALL surf_triangle(d%xyz(:,i,j),d%vertex(:,k,i,j),d%neighbour(:,MOD(k+1+6,6),i,j),S2)
+!            CALL surf_triangle(d%xyz(:,i,j),d%vertex(:,k,i,j),d%neighbour(:,k,i,j),S1)
+!            CALL surf_triangle(d%xyz(:,i,j),d%vertex(:,k,i,j),d%neighbour(:,MOD(k+1+6,6),i,j),S2)
+            CALL surf_triangle(xyz_i(n,:), xyz_v(n+z_pos(k+1),:), xyz_i(n+t_pos(k+1),:),S1)
+            CALL surf_triangle(xyz_i(n,:), xyz_v(n+z_pos(k+1),:), xyz_i(n+t_pos(MOD(k+1+6,6)+1),:),S2)
             Riv(n,k+1)=0.5*(S1+S2)/Ai(n)
             Riv2(n,k+1)=0.5*(S1+S2)/surf_v(k+1)
@@ -236 +501 @@
                  
     ENDDO
+    CALL transfert_request(geom%Ai,req_i1)
+    CALL transfert_request(geom%centroid,req_i1)
     CALL surf_triangle(d%xyz(:,ii_begin,jj_begin),d%xyz(:,ii_begin,jj_end),d%xyz(:,ii_end,jj_begin),S)
 !    PRINT *,"Surf triangle : ",S*20/(4*Pi)