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

Timestamp:

06/14/16 07:04:15 (8 years ago)

Author:

dubos

Message:

New kinetic energy (diagnostic)

File:

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

codes/icosagcm/trunk/src/geometry.f90

@@ -24 +24 @@
     TYPE(t_field),POINTER :: le_de(:) ! le/de, 0. if de=0.
     TYPE(t_field),POINTER :: Riv(:)
+    TYPE(t_field),POINTER :: S1(:)
+    TYPE(t_field),POINTER :: S2(:)
     TYPE(t_field),POINTER :: Riv2(:)
     TYPE(t_field),POINTER :: ne(:)
@@ -70 +72 @@
   REAL(rstd),POINTER :: le(:)          ! lenght of a edge
 !$OMP THREADPRIVATE(le)
-  REAL(rstd),POINTER :: le_de(:)          ! le/de
+  REAL(rstd),POINTER :: le_de(:)       ! le/de
 !$OMP THREADPRIVATE(le_de)
+  REAL(rstd),POINTER :: S1(:,:)        ! area of sub-triangle
+!$OMP THREADPRIVATE(S1)
+  REAL(rstd),POINTER :: S2(:,:)        ! area of sub-tirangle
+!$OMP THREADPRIVATE(S2)
   REAL(rstd),POINTER :: Riv(:,:)       ! weight
 !$OMP THREADPRIVATE(Riv)
@@ -122 +128 @@
     CALL allocate_field(geom%bi,field_t,type_real)
     CALL allocate_field(geom%Av,field_z,type_real)
+    CALL allocate_field(geom%S1,field_t,type_real,6)
+    CALL allocate_field(geom%S2,field_t,type_real,6)
     CALL allocate_field(geom%Riv,field_t,type_real,6)
     CALL allocate_field(geom%Riv2,field_t,type_real,6)
@@ -156 +164 @@
     le_de=geom%le_de(ind)
     Av=geom%Av(ind)
+    S1=geom%S1(ind)
+    S2=geom%S2(ind)
     Riv=geom%Riv(ind)
     Riv2=geom%Riv2(ind)
@@ -420 +430 @@
     INTEGER :: ind,i,j,k,n
     TYPE(t_domain),POINTER :: d
-    REAL(rstd) ::  S
+    REAL(rstd) ::  S12
     REAL(rstd) :: w(6)
     REAL(rstd) :: lon,lat
     INTEGER :: ii_glo,jj_glo
-    REAL(rstd) :: S1,S2
+    REAL(rstd) :: S
           
       
@@ -569 +579 @@
           
           DO k=0,5
-            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)*(radius**2) ! Definition modified for DEC
-            Riv(n,k+1)=0.5*(S1+S2)/Ai(n)
-            Riv2(n,k+1)=0.5*(S1+S2)/surf_v(k+1)
+            CALL surf_triangle(xyz_i(n,:), xyz_v(n+z_pos(k+1),:), xyz_i(n+t_pos(k+1),:),S1(n,k+1) )
+            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(n,k+1) )
+            S12 = .5*(S1(n,k+1)+S2(n,k+1))
+            Riv(n,k+1)=S12/Ai(n)
+            Riv2(n,k+1)=S12/surf_v(k+1)
           ENDDO