Changeset 13 for codes/icosagcm/trunk/src

Timestamp:

05/22/12 18:06:53 (12 years ago)

Author:

ymipsl

Message:

some update

YM

Location:

codes/icosagcm/trunk/src

Files:

• dissip.f90 (modified) (15 diffs)
• etat0_academic.f90 (modified) (1 diff)
• icosa_sw.f90 (modified) (1 diff)
• metric.f90 (modified) (2 diffs)
• spherical_geom.f90 (modified) (1 diff)
• timeloop_sw.f90 (modified) (3 diffs)

codes/icosagcm/trunk/src/dissip.f90

@@ -6 +6 @@
   TYPE(t_field),POINTER,SAVE :: f_gradrot(:)
   TYPE(t_request),POINTER :: req_dissip(:) 
-
+  TYPE(t_field),POINTER,SAVE :: f_due(:)
+ 
   INTEGER,PARAMETER :: nitergdiv=1
   INTEGER,PARAMETER :: nitergrot=1
@@ -23 +24 @@
   IMPLICIT NONE  
     CALL allocate_field(f_gradrot,field_u,type_real)
+    CALL allocate_field(f_due,field_u,type_real)
   END SUBROUTINE allocate_dissip
  
@@ -42 +44 @@
   REAL(rstd)            :: r
   INTEGER :: i,j,n,ind,it
+  REAL :: sum1,sum2
 
     tetagdiv=5000
@@ -184 +187 @@
       dumaxm1=dumax 
       dumax=0
-
+      sum1=0 ; sum2=0
       DO ind=1,ndomain
         u=f_u(ind)
         du=f_du(ind)
-        CALL gradrot(u,du)
+        CALL gradrot(u,du,ind,sum1,sum2)
       ENDDO
 
@@ -225 +228 @@
  
   
-  SUBROUTINE dissip(f_ue,f_due)
+  SUBROUTINE dissip(f_ue)
   USE domain_mod
   USE dimensions
@@ -232 +235 @@
   IMPLICIT NONE
     TYPE(t_field),POINTER :: f_ue(:)
-    TYPE(t_field),POINTER :: f_due(:)
+!    TYPE(t_field),POINTER :: f_due(:)
     REAL(rstd),POINTER         :: ue(:)
     REAL(rstd),POINTER         :: due(:)
@@ -238 +241 @@
     INTEGER :: ind
     REAL :: v
+    REAL :: sum1,sum2
     
     CALL transfert_request(f_ue,req_dissip)
     CALL transfert_request(f_ue,req_dissip)
-
+    sum1=0 ; sum2=0
     DO ind=1,ndomain
       CALL swap_dimensions(ind)
@@ -250 +254 @@
       CALL gradiv(ue,due)
       due=due*dtdissip/tetagdiv
-      CALL gradrot(ue,gradrot_e)
+      CALL gradrot(ue,gradrot_e,ind,sum1,sum2)
       due=due+gradrot_e*dtdissip/tetagrot 
     ENDDO
-  
+    PRINT *,"SUM1 == SUM2 ??  ", sum1,sum2
+    DO ind=1,ndomain
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      ue=f_ue(ind)
+      due=f_due(ind) 
+      
+      ue(:)=ue(:)+0.5*due(:)
+    ENDDO
+      
   END SUBROUTINE dissip
   
@@ -271 +284 @@
       DO i=ii_begin,ii_end
         n=(j-1)*iim+i
-        divu_i(n)=-1./Ai(n)*(ne(n,right)*ue(n+u_right)*le(n+u_right)  +  &
+        divu_i(n)=1./Ai(n)*(ne(n,right)*ue(n+u_right)*le(n+u_right)  +  &
                              ne(n,rup)*ue(n+u_rup)*le(n+u_rup)        +  &  
                              ne(n,lup)*ue(n+u_lup)*le(n+u_lup)        +  &  
@@ -285 +298 @@
         n=(j-1)*iim+i
  
-        gradivu_e(n+u_right)=1/de(n+u_right)*(ne(n,right)*divu_i(n)+ ne(n+t_right,left)*divu_i(n+t_right) )       
-
-        gradivu_e(n+u_lup)=1/de(n+u_lup)*(ne(n,lup)*divu_i(n)+ ne(n+t_lup,rdown)*divu_i(n+t_lup ))       
-    
-        gradivu_e(n+u_ldown)=1/de(n+u_ldown)*(ne(n,ldown)*divu_i(n)+ne(n+t_ldown,rup)*divu_i(n+t_ldown) )
+        gradivu_e(n+u_right)=-1/de(n+u_right)*(ne(n,right)*divu_i(n)+ ne(n+t_right,left)*divu_i(n+t_right) )       
+
+        gradivu_e(n+u_lup)=-1/de(n+u_lup)*(ne(n,lup)*divu_i(n)+ ne(n+t_lup,rdown)*divu_i(n+t_lup ))       
+    
+        gradivu_e(n+u_ldown)=-1/de(n+u_ldown)*(ne(n,ldown)*divu_i(n)+ne(n+t_ldown,rup)*divu_i(n+t_ldown) )
 
       ENDDO
@@ -297 +310 @@
       DO i=ii_begin,ii_end
         n=(j-1)*iim+i
-        gradivu_e(n+u_right)=-gradivu_e(n+u_right)*cdivu
-        gradivu_e(n+u_lup)=-gradivu_e(n+u_lup)*cdivu
-        gradivu_e(n+u_ldown)=-gradivu_e(n+u_ldown)*cdivu
+        gradivu_e(n+u_right)=gradivu_e(n+u_right)*cdivu
+        gradivu_e(n+u_lup)=gradivu_e(n+u_lup)*cdivu
+        gradivu_e(n+u_ldown)=gradivu_e(n+u_ldown)*cdivu
       ENDDO
     ENDDO
@@ -308 +321 @@
   
   
-  SUBROUTINE gradrot(ue,gradrot_e)
+  SUBROUTINE gradrot(ue,gradrot_e,ind,sum1,sum2)
   USE domain_mod
   USE dimensions
@@ -316 +329 @@
     REAL(rstd),INTENT(IN)  :: ue(iim*3*jjm)
     REAL(rstd),INTENT(OUT) :: gradrot_e(iim*3*jjm)
+    INTEGER,INTENT(IN) :: ind
+    REAL,INTENT(OUT) ::sum1,sum2
     REAL(rstd)             :: rot_v(2*iim*jjm)
 
@@ -349 +364 @@
     ENDDO
 
+    
+    DO j=jj_begin,jj_end-1
+      DO i=ii_begin,ii_end-1
+        n=(j-1)*iim+i
+        sum1=sum1+rot_v(n+z_rup)**2*Av(n+z_rup)
+      ENDDO
+    ENDDO
+    
+    DO j=jj_begin+1,jj_end
+      DO i=ii_begin,ii_end-1
+        n=(j-1)*iim+i
+        sum1=sum1+rot_v(n+z_rdown)**2*Av(n+z_rdown)
+      ENDDO
+    ENDDO
+    
+    
+    DO j=jj_begin+1,jj_end-1
+      DO i=ii_begin,ii_end-1
+        n=(j-1)*iim+i
+        sum2=sum2+gradrot_e(n+u_right)*le(n+u_right)*de(n+u_right)*ue(n+u_right)
+      ENDDO
+    ENDDO
+    
+    DO j=jj_begin,jj_end-1
+      DO i=ii_begin+1,ii_end-1
+        n=(j-1)*iim+i
+        sum2=sum2+gradrot_e(n+u_rup)*le(n+u_rup)*de(n+u_rup)*ue(n+u_rup)
+      ENDDO
+    ENDDO
+
+    DO j=jj_begin,jj_end-1
+      DO i=ii_begin+1,ii_end
+        n=(j-1)*iim+i
+        sum2=sum2+gradrot_e(n+u_lup)*le(n+u_lup)*de(n+u_lup)*ue(n+u_lup)
+      ENDDO
+    ENDDO
+
+! frontière
+      j=jj_begin
+      DO i=ii_begin,ii_end-1
+        n=(j-1)*iim+i
+        if (domain(ind)%own(i,j)) sum2=sum2+gradrot_e(n+u_right)*le(n+u_right)*de(n+u_right)*ue(n+u_right)
+      ENDDO
+     
+      j=jj_end
+      DO i=ii_begin,ii_end-1
+        n=(j-1)*iim+i
+        if (domain(ind)%own(i,j)) sum2=sum2+gradrot_e(n+u_right)*le(n+u_right)*de(n+u_right)*ue(n+u_right)
+      ENDDO
+      
+     i=ii_begin
+     DO j=jj_begin,jj_end-1
+        n=(j-1)*iim+i
+        if (domain(ind)%own(i,j))sum2=sum2+gradrot_e(n+u_rup)*le(n+u_rup)*de(n+u_rup)*ue(n+u_rup)
+     ENDDO
+
+     i=ii_end
+     DO j=jj_begin,jj_end-1
+        n=(j-1)*iim+i
+        if (domain(ind)%own(i,j))sum2=sum2+gradrot_e(n+u_rup)*le(n+u_rup)*de(n+u_rup)*ue(n+u_rup)
+     ENDDO
+
+      
+      
+      
     DO j=jj_begin,jj_end
       DO i=ii_begin,ii_end
@@ -359 +439 @@
       ENDDO
     ENDDO  
+    
+    
   END SUBROUTINE 
 

codes/icosagcm/trunk/src/etat0_academic.f90

@@ -157 +157 @@
 
          CALL xyz2lonlat(xyz_e(ij+u_right,:),lon,lat)
-         IF (abs(sin(lat)<1.e-4)) lat=1e-4
+         IF (abs(sin(lat))<1.e-4) lat=1e-4
          x=cos(lat) ; x=x*x ; x=x*x ; x=x*x
          fact(ij+u_right)=(1.-x)/(2.*omega*sin(lat))
 
          CALL xyz2lonlat(xyz_e(ij+u_lup,:),lon,lat)
-         IF (abs(sin(lat)<1.e-4)) lat=1e-4
+         IF (abs(sin(lat))<1.e-4) lat=1e-4
          x=cos(lat) ; x=x*x ; x=x*x ; x=x*x
          fact(ij+u_lup)=(1.-x)/(2.*omega*sin(lat))
 
          CALL xyz2lonlat(xyz_e(ij+u_ldown,:),lon,lat)
-         IF (abs(sin(lat)<1.e-4)) lat=1e-4
+         IF (abs(sin(lat))<1.e-4) lat=1e-4
          x=cos(lat) ; x=x*x ; x=x*x ; x=x*x
          fact(ij+u_ldown)=(1.-x)/(2.*omega*sin(lat))

codes/icosagcm/trunk/src/icosa_sw.f90

@@ -53 +53 @@
 !  CALL WriteField("Ai",geom%Ai)
 !  CALL WriteField("sum_ne",sum_ne)
-
  CALL timeloop
   

codes/icosagcm/trunk/src/metric.f90

@@ -194 +194 @@
    ENDDO
 
-  CALL dist_cart(vertex_glo(1,1,1)%xyz,vertex_glo(2,1,1)%xyz,d1) 
-  CALL dist_cart(vertex_glo(2,1,1)%xyz,vertex_glo(3,1,1)%xyz,d2) 
-  CALL dist_cart(vertex_glo(3,1,1)%xyz,vertex_glo(4,1,1)%xyz,d3) 
-  CALL div_arc(vertex_glo(1,1,1)%xyz,vertex_glo(3,1,1)%xyz,0.5,p1)
+!  CALL dist_cart(vertex_glo(1,1,1)%xyz,vertex_glo(2,1,1)%xyz,d1) 
+!  CALL dist_cart(vertex_glo(2,1,1)%xyz,vertex_glo(3,1,1)%xyz,d2) 
+!  CALL dist_cart(vertex_glo(3,1,1)%xyz,vertex_glo(4,1,1)%xyz,d3) 
+!  CALL div_arc(vertex_glo(1,1,1)%xyz,vertex_glo(3,1,1)%xyz,0.5,p1)
 !  CALL div_arc(vertex_glo(2,1,1)%xyz,vertex_glo(1,3,1)%xyz,1./3,p1)
 !  CALL div_arc(vertex_glo(1,2,1)%xyz,vertex_glo(3,1,1)%xyz,1./3,p2)
@@ -206 +206 @@
 !  PRINT *,"dist",vertex_glo(2,1,1)%xyz
 !  PRINT *,"dist",p1/sqrt(sum(p1**2))
-  CALL Centroide(vertex_glo(1,1,1)%xyz,vertex_glo(2,1,1)%xyz,vertex_glo(1,2,1)%xyz,p1)
-  CALL Centroide(vertex_glo(1,2,1)%xyz,vertex_glo(2,1,1)%xyz,vertex_glo(1,1,1)%xyz,p1)
-  CALL Centroide(vertex_glo(1,2,1)%xyz,vertex_glo(1,1,1)%xyz,vertex_glo(2,1,1)%xyz,p1)
-  CALL dist_cart(vertex_glo(1,1,1)%xyz,p1,d1) 
-  CALL dist_cart(vertex_glo(2,1,1)%xyz,p1,d2) 
-  CALL dist_cart(vertex_glo(1,2,1)%xyz,p1,d3) 
+!  CALL Centroide(vertex_glo(1,1,1)%xyz,vertex_glo(2,1,1)%xyz,vertex_glo(1,2,1)%xyz,p1)
+!  CALL Centroide(vertex_glo(1,2,1)%xyz,vertex_glo(2,1,1)%xyz,vertex_glo(1,1,1)%xyz,p1)
+!  CALL Centroide(vertex_glo(1,2,1)%xyz,vertex_glo(1,1,1)%xyz,vertex_glo(2,1,1)%xyz,p1)
+!  CALL dist_cart(vertex_glo(1,1,1)%xyz,p1,d1) 
+!  CALL dist_cart(vertex_glo(2,1,1)%xyz,p1,d2) 
+!  CALL dist_cart(vertex_glo(1,2,1)%xyz,p1,d3) 
 !  PRINT *, "dist",d1
 !  PRINT *, "dist",d2

codes/icosagcm/trunk/src/spherical_geom.f90

@@ -125 +125 @@
   M(2,1)=xb ; M(2,2)=yb ; M(2,3)=zb
   M(3,1)=xc ; M(3,2)=yc ; M(3,3)=zc
-  
-  CALL DGESV(3,1,M,3,IPIV,C,3,info)
+  stop 'STOP'
+!  CALL DGESV(3,1,M,3,IPIV,C,3,info)
   
 END SUBROUTINE div_arc

codes/icosagcm/trunk/src/timeloop_sw.f90

@@ -70 +70 @@
 
   CALL allocate_caldyn
-  
+  CALL init_dissip(dt)
 
   CALL etat0_williamson(f_h,f_u)
@@ -78 +78 @@
 
     CALL caldyn(f_h, f_u, f_dh, f_du)
+    
 
     IF (scheme==Euler) THEN
@@ -89 +90 @@
     ENDIF
 
+    CALL dissip(f_u)   
+
   ENDDO