Index: codes/icosagcm/trunk/src/caldyn_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 122)
+++ codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 123)
@@ -179,22 +179,4 @@
     REAL(rstd) :: etav,hv, du2
 
-!   REAL(rstd) :: theta(iim*jjm,llm)  ! potential temperature
-!   REAL(rstd) :: p(iim*jjm,llm+1)  ! pression
-!   REAL(rstd) :: pk(iim*jjm,llm)   ! Exner function
-!   REAL(rstd) :: pks(iim*jjm)
-!! Intermediate variable to compute exner function
-!   REAL(rstd) :: alpha(iim*jjm,llm)
-!   REAL(rstd) :: beta(iim*jjm,llm)
-!!   
-!   REAL(rstd) :: phi(iim*jjm,llm)    ! geopotential
-!   REAL(rstd) :: mass(iim*jjm,llm)   ! mass    
-!   REAL(rstd) :: rhodz(iim*jjm,llm)   ! mass density   
-!   REAL(rstd) :: Fe(3*iim*jjm,llm)   ! mass flux   
-!   REAL(rstd) :: Ftheta(3*iim*jjm,llm) ! theta flux   
-!   REAL(rstd) :: convm(iim*jjm,llm)    ! mass flux convergence
-!   REAL(rstd) :: w(iim*jjm,llm)       ! vertical velocity      
-!   REAL(rstd) :: qv(2*iim*jjm,llm)       ! potential velocity  
-!   REAL(rstd) :: berni(iim*jjm,llm)   ! bernouilli term
-    
     REAL(rstd),ALLOCATABLE,SAVE :: theta(:,:)  ! potential temperature
     REAL(rstd),ALLOCATABLE,SAVE :: p(:,:)  ! pression
@@ -206,5 +188,4 @@
  !   
     REAL(rstd),ALLOCATABLE,SAVE :: phi(:,:)    ! geopotential
-    REAL(rstd),ALLOCATABLE,SAVE :: mass(:,:)   ! mass    
     REAL(rstd),ALLOCATABLE,SAVE :: rhodz(:,:)   ! mass density   
     REAL(rstd),ALLOCATABLE,SAVE :: Fe(:,:)   ! mass flux   
@@ -228,5 +209,4 @@
       ALLOCATE(beta(iim*jjm,llm))
       ALLOCATE(phi(iim*jjm,llm))    ! geopotential
-      ALLOCATE(mass(iim*jjm,llm))   ! mass    
       ALLOCATE(rhodz(iim*jjm,llm))   ! mass density   
       ALLOCATE(Fe(3*iim*jjm,llm))   ! mass flux   
@@ -269,174 +249,17 @@
       ENDDO
     ENDDO
-      
- !!! Compute Exner function
-!    PRINT *, 'Computing Exner'
-    CALL compute_exner(ps,p,pks,pk,1) 
 
 !!! Compute mass
-   ! PRINT *, 'Computing mass'
+   ! PRINT *, 'Computing mass, theta'
    DO l = 1, llm
 !$OMP DO
-     DO j=jj_begin-1,jj_end+1
-       DO i=ii_begin-1,ii_end+1
-         ij=(j-1)*iim+i
-         mass(ij,l) = ( p(ij,l) - p(ij,l+1) ) * Ai(ij)/g
-         rhodz(ij,l) = mass(ij,l) / Ai(ij)
-       ENDDO
-     ENDDO
+      DO j=jj_begin-1,jj_end+1
+         DO i=ii_begin-1,ii_end+1
+            ij=(j-1)*iim+i
+            rhodz(ij,l) = ( p(ij,l) - p(ij,l+1) )/g
+            theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
+         ENDDO
+      ENDDO
    ENDDO
-
-!! compute theta
-   ! PRINT *, 'Computing theta'
-    DO    l    = 1, llm
-!$OMP DO
-      DO j=jj_begin-1,jj_end+1
-        DO i=ii_begin-1,ii_end+1
-          ij=(j-1)*iim+i
-          theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
-        ENDDO
-      ENDDO
-    ENDDO
-
-
-!!! Compute geopotential
-
-  ! for first layer
-!$OMP DO
-   DO j=jj_begin-1,jj_end+1
-     DO i=ii_begin-1,ii_end+1
-       ij=(j-1)*iim+i
-       phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) )
-     ENDDO
-   ENDDO
-          
-  ! for other layers
-   DO l = 2, llm
-!$OMP DO
-     DO j=jj_begin-1,jj_end+1
-       DO i=ii_begin-1,ii_end+1
-         ij=(j-1)*iim+i
-         phi(ij,l) = phi(ij,l-1) + 0.5 * ( theta(ij,l)  + theta(ij,l-1) )  & 
-                                       * (  pk(ij,l-1) -  pk(ij,l)    )
-       ENDDO
-     ENDDO
-   ENDDO
-
-   
-!!!  Compute mass flux
-
-  DO l = 1, llm
-!$OMP DO
-    DO j=jj_begin-1,jj_end+1
-      DO i=ii_begin-1,ii_end+1
-        ij=(j-1)*iim+i
-        Fe(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right)
-        Fe(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup)
-        Fe(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown)
-      ENDDO
-    ENDDO
-  ENDDO
-
-!!! fisrt composante dtheta
-
- ! Flux on the edge 
-  DO l = 1, llm
-!$OMP DO
-    DO j=jj_begin-1,jj_end+1
-      DO i=ii_begin-1,ii_end+1
-        ij=(j-1)*iim+i  
-        Ftheta(ij+u_right,l)=0.5*(theta(ij,l)+theta(ij+t_right,l))*Fe(ij+u_right,l)
-        Ftheta(ij+u_lup,l)=0.5*(theta(ij,l)+theta(ij+t_lup,l))*Fe(ij+u_lup,l)
-        Ftheta(ij+u_ldown,l)=0.5*(theta(ij,l)+theta(ij+t_ldown,l))*Fe(ij+u_ldown,l)
-      ENDDO
-    ENDDO
-  ENDDO
- 
- 
- ! compute divergence 
-  DO l = 1, llm
-!$OMP DO
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i  
-! signe ? attention d (rho theta dz)
-        ! dtheta_rhodz =  -div(flux.theta)
-        dtheta_rhodz(ij,l)=-1./Ai(ij)*(ne(ij,right)*Ftheta(ij+u_right,l)    +  &
-                                 ne(ij,rup)*Ftheta(ij+u_rup,l)        +  &  
-                                 ne(ij,lup)*Ftheta(ij+u_lup,l)        +  &  
-                                 ne(ij,left)*Ftheta(ij+u_left,l)      +  &  
-                                 ne(ij,ldown)*Ftheta(ij+u_ldown,l)    +  &  
-                                 ne(ij,rdown)*Ftheta(ij+u_rdown,l))
-      ENDDO
-    ENDDO
-  ENDDO
-
-
-
-!!! mass flux convergence computation  
-
- ! horizontal convergence
-  DO l = 1, llm
-!$OMP DO
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
-        ! convm = +div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
-        convm(ij,l)= 1./Ai(ij)*(ne(ij,right)*Fe(ij+u_right,l)   +  &
-                                ne(ij,rup)*Fe(ij+u_rup,l)       +  &  
-                                ne(ij,lup)*Fe(ij+u_lup,l)       +  &  
-                                ne(ij,left)*Fe(ij+u_left,l)     +  &  
-                                ne(ij,ldown)*Fe(ij+u_ldown,l)   +  &  
-                                ne(ij,rdown)*Fe(ij+u_rdown,l))    
-       ENDDO
-     ENDDO
-   ENDDO
- 
-   
- ! vertical integration from up to down
-  DO  l = llm-1, 1, -1
-!$OMP DO
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
-        convm(ij,l) = convm(ij,l) + convm(ij,l+1)
-      ENDDO
-    ENDDO
-  ENDDO
-
-  
-!!! Compute dps
-!$OMP DO
-  DO j=jj_begin,jj_end
-    DO i=ii_begin,ii_end
-      ij=(j-1)*iim+i
-      ! dps/dt = -int(div flux)dz
-      dps(ij)=-convm(ij,1) * g 
-    ENDDO
-  ENDDO
-
-
-!!! Compute vertical velocity
-  DO l = 1,llm-1
-!$OMP DO
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
-        ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt
-        ! => w>0 for upward transport
-        w( ij, l+1 ) = convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 )
-      ENDDO
-    ENDDO
-  ENDDO
-
-!$OMP DO
-  ! vertical mass flux at the surface = 0
-  DO j=jj_begin,jj_end
-    DO i=ii_begin,ii_end
-      ij=(j-1)*iim+i
-      w(ij,1)  = 0.
-    ENDDO
-  ENDDO
-
 
 !!! Compute shallow-water potential vorticity
@@ -471,4 +294,81 @@
     ENDDO
 
+
+  DO l = 1, llm
+!!!  Compute mass and theta fluxes
+!$OMP DO
+    DO j=jj_begin-1,jj_end+1
+      DO i=ii_begin-1,ii_end+1
+        ij=(j-1)*iim+i
+        Fe(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right)
+        Fe(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup)
+        Fe(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown)
+
+        Ftheta(ij+u_right,l)=0.5*(theta(ij,l)+theta(ij+t_right,l))*Fe(ij+u_right,l)
+        Ftheta(ij+u_lup,l)=0.5*(theta(ij,l)+theta(ij+t_lup,l))*Fe(ij+u_lup,l)
+        Ftheta(ij+u_ldown,l)=0.5*(theta(ij,l)+theta(ij+t_ldown,l))*Fe(ij+u_ldown,l)
+      ENDDO
+    ENDDO
+
+!!! compute horizontal divergence of fluxes
+!$OMP DO
+    DO j=jj_begin,jj_end
+      DO i=ii_begin,ii_end
+        ij=(j-1)*iim+i  
+        ! convm = +div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
+        convm(ij,l)= 1./Ai(ij)*(ne(ij,right)*Fe(ij+u_right,l)   +  &
+                                ne(ij,rup)*Fe(ij+u_rup,l)       +  &  
+                                ne(ij,lup)*Fe(ij+u_lup,l)       +  &  
+                                ne(ij,left)*Fe(ij+u_left,l)     +  &  
+                                ne(ij,ldown)*Fe(ij+u_ldown,l)   +  &  
+                                ne(ij,rdown)*Fe(ij+u_rdown,l))    
+
+        ! signe ? attention d (rho theta dz)
+        ! dtheta_rhodz =  -div(flux.theta)
+        dtheta_rhodz(ij,l)=-1./Ai(ij)*(ne(ij,right)*Ftheta(ij+u_right,l)    +  &
+                                 ne(ij,rup)*Ftheta(ij+u_rup,l)        +  &  
+                                 ne(ij,lup)*Ftheta(ij+u_lup,l)        +  &  
+                                 ne(ij,left)*Ftheta(ij+u_left,l)      +  &  
+                                 ne(ij,ldown)*Ftheta(ij+u_ldown,l)    +  &  
+                                 ne(ij,rdown)*Ftheta(ij+u_rdown,l))
+      ENDDO
+    ENDDO
+  ENDDO
+   
+!!! cumulate mass flux convergence from top to bottom
+  DO  l = llm-1, 1, -1
+!$OMP DO
+    DO j=jj_begin,jj_end
+      DO i=ii_begin,ii_end
+        ij=(j-1)*iim+i
+        convm(ij,l) = convm(ij,l) + convm(ij,l+1)
+      ENDDO
+    ENDDO
+  ENDDO
+  
+!!! Compute vertical mass flux
+  DO l = 1,llm-1
+!$OMP DO
+    DO j=jj_begin,jj_end
+      DO i=ii_begin,ii_end
+        ij=(j-1)*iim+i
+        ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt
+        ! => w>0 for upward transport
+        w( ij, l+1 ) = convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 )
+      ENDDO
+    ENDDO
+  ENDDO
+
+! compute dps, vertical mass flux at the surface = 0
+!$OMP DO
+  DO j=jj_begin,jj_end
+    DO i=ii_begin,ii_end
+      ij=(j-1)*iim+i
+      w(ij,1)  = 0.
+      ! dps/dt = -int(div flux)dz
+      dps(ij)=-convm(ij,1) * g 
+    ENDDO
+  ENDDO
+
 !!! Compute potential vorticity (Coriolis) contribution to du
   DO l=1,llm
@@ -521,5 +421,32 @@
   ENDDO
   
-
+!!! Compute Exner function
+!    PRINT *, 'Computing Exner'
+    CALL compute_exner(ps,p,pks,pk,1) 
+
+!!! Compute geopotential
+
+  ! for first layer
+!$OMP DO
+   DO j=jj_begin-1,jj_end+1
+     DO i=ii_begin-1,ii_end+1
+       ij=(j-1)*iim+i
+       phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) )
+     ENDDO
+   ENDDO
+          
+  ! for other layers
+   DO l = 2, llm
+!$OMP DO
+     DO j=jj_begin-1,jj_end+1
+       DO i=ii_begin-1,ii_end+1
+         ij=(j-1)*iim+i
+         phi(ij,l) = phi(ij,l-1) + 0.5 * ( theta(ij,l)  + theta(ij,l-1) )  & 
+                                       * (  pk(ij,l-1) -  pk(ij,l)    )
+       ENDDO
+     ENDDO
+   ENDDO
+
+   
 !!! Compute bernouilli term = Kinetic Energy + geopotential
   DO l=1,llm
@@ -542,30 +469,5 @@
   
  
-!!! second contribution to du (gradients of Bernoulli and Exner functions)
-  DO l=1,llm
-!$OMP DO
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
-        
-         du(ij+u_right,l)= du(ij+u_right,l)+ 1/de(ij+u_right) * (                              &
-                           0.5*(theta(ij,l)+theta(ij+t_right,l))                             &
-                              *( ne(ij,right)*pk(ij,l)+ne(ij+t_right,left)*pk(ij+t_right,l)) &
-                         + ne(ij,right)*berni(ij,l)+ne(ij+t_right,left)*berni(ij+t_right,l) )
-         
-         du(ij+u_lup,l)= du(ij+u_lup,l)+ 1/de(ij+u_lup) * (                                  &
-                           0.5*(theta(ij,l)+theta(ij+t_lup,l))                             &
-                              *( ne(ij,lup)*pk(ij,l)+ne(ij+t_lup,rdown)*pk(ij+t_lup,l))    &
-                         + ne(ij,lup)*berni(ij,l)+ne(ij+t_lup,rdown)*berni(ij+t_lup,l) )
-                         
-         du(ij+u_ldown,l)= du(ij+u_ldown,l)+ 1/de(ij+u_ldown) * (                                &
-                           0.5*(theta(ij,l)+theta(ij+t_ldown,l))                               &
-                              *( ne(ij,ldown)*pk(ij,l)+ne(ij+t_ldown,rup)*pk(ij+t_ldown,l))    &
-                         + ne(ij,ldown)*berni(ij,l)+ne(ij+t_ldown,rup)*berni(ij+t_ldown,l) )
-      ENDDO
-    ENDDO
-  ENDDO
- 
-!!! save second contribution to du for debugging output
+!!! gradients of Bernoulli and Exner functions
   DO l=1,llm
 !$OMP DO
@@ -579,40 +481,37 @@
                         + ne(ij,right)*berni(ij,l)+ne(ij+t_right,left)*berni(ij+t_right,l) )
         
+        du(ij+u_right,l) = du(ij+u_right,l) + out_u(ij+u_right,l)
+        
         out_u(ij+u_lup,l)=  1/de(ij+u_lup) * (                                  &
                           0.5*(theta(ij,l)+theta(ij+t_lup,l))                             &
                              *( ne(ij,lup)*pk(ij,l)+ne(ij+t_lup,rdown)*pk(ij+t_lup,l))    &
                         + ne(ij,lup)*berni(ij,l)+ne(ij+t_lup,rdown)*berni(ij+t_lup,l) )
-                        
+        
+        du(ij+u_lup,l) = du(ij+u_lup,l) + out_u(ij+u_lup,l)
+                
         out_u(ij+u_ldown,l)= 1/de(ij+u_ldown) * (                                &
                           0.5*(theta(ij,l)+theta(ij+t_ldown,l))                               &
                              *( ne(ij,ldown)*pk(ij,l)+ne(ij+t_ldown,rup)*pk(ij+t_ldown,l))    &
                         + ne(ij,ldown)*berni(ij,l)+ne(ij+t_ldown,rup)*berni(ij+t_ldown,l) )
+        
+        du(ij+u_ldown,l) = du(ij+u_ldown,l) + out_u(ij+u_ldown,l)
+
       ENDDO
     ENDDO
   ENDDO
-
-!!! contributions due to vertical advection
  
- ! Contribution to dtheta
+!!! contributions of vertical advection to du, dtheta
+ 
   DO l=1,llm-1
 !$OMP DO
     DO j=jj_begin,jj_end
       DO i=ii_begin,ii_end
+        ij=(j-1)*iim+i
          ! ww>0 <=> upward transport
-        ij=(j-1)*iim+i
+
         ww            =   0.5 * w(ij,l+1) * (theta(ij,l) + theta(ij,l+1) )
         dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  -  ww
         dtheta_rhodz(ij,l+1) = dtheta_rhodz(ij,l+1)  +  ww
-      ENDDO
-    ENDDO
-  ENDDO       
-
-
- ! Contribution to du
-  DO l=1,llm-1
-!$OMP DO
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
+
         ww  = 0.5 * ( w(ij,l+1) + w(ij+t_right,l+1))
         uu  = u(ij+u_right,l+1) - u(ij+u_right,l)  
@@ -643,5 +542,5 @@
     DEALLOCATE(beta)
     DEALLOCATE(phi)    ! geopotential
-    DEALLOCATE(mass)   ! mass    
+!    DEALLOCATE(mass)   ! mass    
     DEALLOCATE(rhodz)   ! mass density   
     DEALLOCATE(Fe)   ! mass flux   
