Changeset 427 for codes

Timestamp:

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

Author:

dubos

Message:

New kinetic energy (diagnostic)

Location:

codes/icosagcm/trunk/src

Files:

• geometry.f90 (modified) (6 diffs)
• kinetic.f90 (modified) (3 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
 

codes/icosagcm/trunk/src/kinetic.f90

@@ -1 +1 @@
 MODULE kinetic_mod
+IMPLICIT NONE
+PRIVATE
 
+PUBLIC :: kinetic, kinetic_v, kinetic_new, gradient
 
 CONTAINS
@@ -24 +27 @@
       Ki=f_Ki(ind)
       CALL compute_kinetic(ue, Ki)
+    ENDDO  
+  END SUBROUTINE kinetic
+  
+  SUBROUTINE kinetic_new(f_ue,f_Ki)
+    USE icosa
+    IMPLICIT NONE
+    TYPE(t_field), POINTER :: f_ue(:)
+    TYPE(t_field), POINTER :: f_Ki(:)
+    
+    REAL(rstd), POINTER :: ue(:,:)
+    REAL(rstd), POINTER :: Ki(:,:)
+    INTEGER :: ind
+    
+    CALL transfert_request(f_ue,req_e1_vect)
+    CALL transfert_request(f_ue,req_e1_vect)
+    
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       ue=f_ue(ind)
+       Ki=f_Ki(ind)
+       CALL compute_Ki_new(ue, Ki)
     ENDDO
+  END SUBROUTINE kinetic_new
   
-  END SUBROUTINE kinetic
+  SUBROUTINE kinetic_v(f_ue,f_Kv)
+  USE icosa
+  IMPLICIT NONE
+    TYPE(t_field), POINTER :: f_ue(:)
+    TYPE(t_field), POINTER :: f_Kv(:)
+  
+    REAL(rstd), POINTER :: ue(:,:)
+    REAL(rstd), POINTER :: Kv(:,:)
+    INTEGER :: ind
+  
+    CALL transfert_request(f_ue,req_e1_vect)
+    CALL transfert_request(f_ue,req_e1_vect)
+
+    DO ind=1,ndomain
+      IF (.NOT. assigned_domain(ind)) CYCLE
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      ue=f_ue(ind)
+      Kv=f_Kv(ind)
+      CALL compute_kv(ue, Kv)
+    ENDDO  
+  END SUBROUTINE kinetic_v
   
   SUBROUTINE compute_kinetic(ue, Ki)
@@ -51 +99 @@
       ENDDO
     ENDDO
+  END SUBROUTINE compute_kinetic
+  
+  SUBROUTINE compute_kv(ue, Kv)
+  USE icosa
+  USE omp_para
+  IMPLICIT NONE
+    REAL(rstd),INTENT(IN) :: ue(3*iim*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: Kv(2*iim*jjm,llm)
+    INTEGER :: ij,l, u_up, u_down
+
+    u_up    = t_lup + u_right
+    u_down  = t_rdown + u_left
     
+    DO l=ll_begin,ll_end
+      DO ij=ij_begin,ij_end
+          Kv(ij+z_up,l) = (radius**2/Av(ij+z_up))*(                         &
+                               S1(ij,vup)*ue(ij+u_rup,l)**2 +        &
+                               S2(ij,vup)*ue(ij+u_lup,l)**2 +        &
+                               S2(ij+t_lup,vrdown)*ue(ij+u_up,l)**2)
+
+          Kv(ij+z_down,l) = (radius**2/Av(ij+z_down))*(                         &
+                               S1(ij,vdown)*ue(ij+u_ldown,l)**2 +       &
+                               S2(ij,vdown)*ue(ij+u_rdown,l)**2 +       &
+                               S2(ij+t_rdown,vlup)*ue(ij+u_down,l)**2 )
+       ENDDO
+    ENDDO
+  END SUBROUTINE compute_kv
   
-  END SUBROUTINE compute_kinetic
-        
+  SUBROUTINE compute_Ki_new(ue, Ki)
+    USE icosa
+    USE omp_para
+    IMPLICIT NONE
+    REAL(rstd),INTENT(IN) :: ue(3*iim*jjm,llm)
+    REAL(rstd),INTENT(OUT):: Ki(iim*jjm,llm)
+    REAL(rstd) :: Kv(2*iim*jjm,llm)
+    INTEGER :: ij,l, u_up, u_down
+    
+    CALL compute_kv(ue,Kv)
+    
+    DO l=ll_begin,ll_end
+       DO ij=ij_begin,ij_end
+          Ki(ij,l) = Riv(ij,vup)*Kv(ij+z_up,l) + &
+               Riv(ij,vlup)  *Kv(ij+z_lup,l) + &
+               Riv(ij,vldown)*Kv(ij+z_ldown,l) + &
+               Riv(ij,vdown) *Kv(ij+z_down,l) + &
+               Riv(ij,vrdown)*Kv(ij+z_rdown,l) + &
+               Riv(ij,vrup)  *Kv(ij+z_rup,l)
+       END DO
+    END DO
+  END SUBROUTINE compute_Ki_new
+  
+  SUBROUTINE gradient(f_berni, f_du)
+  USE icosa
+  IMPLICIT NONE
+    TYPE(t_field), POINTER :: f_berni(:)
+    TYPE(t_field), POINTER :: f_du(:)
+  
+    REAL(rstd), POINTER :: du(:,:)
+    REAL(rstd), POINTER :: berni(:,:)
+    INTEGER :: ind
+  
+    CALL transfert_request(f_du,req_e1_vect)
+    CALL transfert_request(f_du,req_e1_vect)
+
+    DO ind=1,ndomain
+      IF (.NOT. assigned_domain(ind)) CYCLE
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      berni=f_berni(ind)
+      du=f_du(ind)
+      CALL compute_grad(berni, du)
+    ENDDO
+  
+  END SUBROUTINE gradient
+  
+  SUBROUTINE compute_grad(berni, du)
+  USE icosa
+  USE omp_para
+  IMPLICIT NONE
+    REAL(rstd),INTENT(IN) :: berni(iim*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: du(3*iim*jjm,llm)
+    INTEGER :: i,j,ij,l
+    
+    DO l=ll_begin,ll_end
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+            du(ij+u_right,l) = ne_right*(berni(ij,l)-berni(ij+t_right,l))/de(ij+u_right)
+            du(ij+u_lup,l)   = ne_lup*(berni(ij,l)-berni(ij+t_lup,l))/de(ij+u_right)
+            du(ij+u_ldown,l) = ne_ldown*(berni(ij,l)-berni(ij+t_ldown,l))/de(ij+u_right)
+         ENDDO
+      ENDDO
+    ENDDO
+
+  END SUBROUTINE compute_grad
+
+
 END MODULE kinetic_mod