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

Timestamp:

09/25/15 14:36:36 (9 years ago)

Author:

dubos

Message:

Introduced DEC convention for velocity - HEVI only - cleanup to follow

File:

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

codes/icosagcm/trunk/src/geometry.f90

@@ -22 +22 @@
     TYPE(t_field),POINTER :: de(:)
     TYPE(t_field),POINTER :: le(:)
+    TYPE(t_field),POINTER :: le_de(:) ! le/de, 0. if de=0.
     TYPE(t_field),POINTER :: Riv(:)
     TYPE(t_field),POINTER :: Riv2(:)
@@ -69 +70 @@
   REAL(rstd),POINTER :: le(:)          ! lenght of a edge
 !$OMP THREADPRIVATE(le)
+  REAL(rstd),POINTER :: le_de(:)          ! le/de
+!$OMP THREADPRIVATE(le_de)
   REAL(rstd),POINTER :: Riv(:,:)       ! weight
 !$OMP THREADPRIVATE(Riv)
@@ -116 +119 @@
     CALL allocate_field(geom%de,field_u,type_real)
     CALL allocate_field(geom%le,field_u,type_real)
+    CALL allocate_field(geom%le_de,field_u,type_real)
     CALL allocate_field(geom%bi,field_t,type_real)
     CALL allocate_field(geom%Av,field_z,type_real)
@@ -150 +154 @@
     de=geom%de(ind)
     le=geom%le(ind)
+    le_de=geom%le_de(ind)
     Av=geom%Av(ind)
     Riv=geom%Riv(ind)
@@ -438 +443 @@
           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))
-         
+
+          le_de(n+u_right)=le(n+u_right)/de(n+u_right) ! NaN possible but should be harmless
+          le_de(n+u_lup)  =le(n+u_lup)  /de(n+u_lup)
+          le_de(n+u_ldown)=le(n+u_ldown)/de(n+u_ldown)
+
           Ai(n)=0
           DO k=0,5
@@ -528 +537 @@
             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)