Index: codes/icosagcm/trunk/src/dissip/dissip_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/dissip/dissip_gcm.f90	(revision 899)
+++ codes/icosagcm/trunk/src/dissip/dissip_gcm.f90	(revision 900)
@@ -199,6 +199,6 @@
           u=f_u(ind)
           du=f_du(ind)
-          CALL compute_gradiv(u,du,1,1)
-          u=du
+          CALL compute_gradiv_inplace(u,1,1)
+          du=u
         ENDDO
       ENDDO
@@ -292,6 +292,6 @@
           u=f_u(ind)
           du=f_du(ind)
-          CALL compute_gradrot(u,du,1,1)
-          u=du
+          CALL compute_gradrot_inplace(u,1,1)
+          du=u
         ENDDO
       ENDDO
@@ -380,6 +380,6 @@
           theta=f_theta(ind)
           dtheta=f_dtheta(ind)
-          CALL compute_divgrad(theta,dtheta,1,1)
-          theta=dtheta
+          CALL compute_divgrad_inplace(theta,1,1)
+          dtheta=theta
         ENDDO
       ENDDO
@@ -673,5 +673,5 @@
         CALL swap_geometry(ind)
         due=f_due(ind) 
-        CALL compute_gradiv(due,due,ll_begin,ll_end)
+        CALL compute_gradiv_inplace(due,ll_begin,ll_end)
       ENDDO
     ENDDO
@@ -722,5 +722,5 @@
         CALL swap_geometry(ind)
         due=f_due(ind) 
-        CALL compute_gradrot(due,due,ll_begin,ll_end)
+        CALL compute_gradrot_inplace(due,ll_begin,ll_end)
       ENDDO
 
@@ -763,5 +763,5 @@
         CALL swap_geometry(ind)
         dtheta=f_dtheta(ind) 
-        CALL compute_divgrad(dtheta,dtheta,ll_begin,ll_end)
+        CALL compute_divgrad_inplace(dtheta,ll_begin,ll_end)
       ENDDO
 
@@ -814,5 +814,5 @@
         CALL swap_geometry(ind)
         dtheta_rhodz=f_dtheta_rhodz(ind) 
-        CALL compute_divgrad(dtheta_rhodz,dtheta_rhodz,ll_begin,ll_end)
+        CALL compute_divgrad_inplace(dtheta_rhodz,ll_begin,ll_end)
       ENDDO
 
@@ -846,144 +846,149 @@
     REAL(rstd),INTENT(IN)  :: ue(iim*3*jjm,llm)
     REAL(rstd),INTENT(OUT) :: gradivu_e(iim*3*jjm,llm)
+    
+    gradivu_e = ue
+    CALL compute_gradiv_inplace(gradivu_e, llb, lle)
+  END SUBROUTINE compute_gradiv
+  
+  SUBROUTINE compute_gradiv_inplace(ue_gradivu_e,llb,lle)
+  USE icosa
+  IMPLICIT NONE
+    INTEGER,INTENT(IN)     :: llb
+    INTEGER,INTENT(IN)     :: lle
+    REAL(rstd),INTENT(INOUT) :: ue_gradivu_e(iim*3*jjm,llm)
     REAL(rstd) :: divu_i(iim*jjm,llb:lle)
     
     INTEGER :: ij,l
 
+    DO l=llb,lle
+      !DIR$ SIMD
+      DO ij=ij_begin,ij_end
+          divu_i(ij,l)=1./Ai(ij)*(ne(ij,right)*ue_gradivu_e(ij+u_right,l)*le(ij+u_right)  +  &
+                             ne(ij,rup)*ue_gradivu_e(ij+u_rup,l)*le(ij+u_rup)        +  &  
+                             ne(ij,lup)*ue_gradivu_e(ij+u_lup,l)*le(ij+u_lup)        +  &  
+                             ne(ij,left)*ue_gradivu_e(ij+u_left,l)*le(ij+u_left)     +  &  
+                             ne(ij,ldown)*ue_gradivu_e(ij+u_ldown,l)*le(ij+u_ldown)  +  &  
+                             ne(ij,rdown)*ue_gradivu_e(ij+u_rdown,l)*le(ij+u_rdown))
+      ENDDO
+    ENDDO
+    
+    DO l=llb,lle
+      !DIR$ SIMD
+      DO ij=ij_begin,ij_end
+          ue_gradivu_e(ij+u_right,l)=-1/de(ij+u_right)*(ne(ij,right)*divu_i(ij,l)+ ne(ij+t_right,left)*divu_i(ij+t_right,l) ) 
+          ue_gradivu_e(ij+u_lup,l)=-1/de(ij+u_lup)*(ne(ij,lup)*divu_i(ij,l)+ ne(ij+t_lup,rdown)*divu_i(ij+t_lup,l))   
+          ue_gradivu_e(ij+u_ldown,l)=-1/de(ij+u_ldown)*(ne(ij,ldown)*divu_i(ij,l)+ne(ij+t_ldown,rup)*divu_i(ij+t_ldown,l) )
+      ENDDO
+    ENDDO
+
+    DO l=llb,lle
+      !DIR$ SIMD
+      DO ij=ij_begin,ij_end
+          ue_gradivu_e(ij+u_right,l)=-ue_gradivu_e(ij+u_right,l)*cgraddiv
+          ue_gradivu_e(ij+u_lup,l)=-ue_gradivu_e(ij+u_lup,l)*cgraddiv
+          ue_gradivu_e(ij+u_ldown,l)=-ue_gradivu_e(ij+u_ldown,l)*cgraddiv
+      ENDDO
+    ENDDO   
+  END SUBROUTINE compute_gradiv_inplace
+  
+  SUBROUTINE compute_divgrad(theta,divgrad_i,llb,lle)
+  USE icosa
+  IMPLICIT NONE
+    INTEGER,INTENT(IN)     :: llb
+    INTEGER,INTENT(IN)     :: lle
+    REAL(rstd),INTENT(IN) :: theta(iim*jjm,1:lle)
+    REAL(rstd),INTENT(OUT) :: divgrad_i(iim*jjm,1:lle)
+    
+    divgrad_i = theta
+    CALL compute_divgrad_inplace(divgrad_i, llb, lle)
+  END SUBROUTINE
+  
+  SUBROUTINE compute_divgrad_inplace(theta_divgrad_i,llb,lle)
+  USE icosa
+  IMPLICIT NONE
+    INTEGER,INTENT(IN)     :: llb
+    INTEGER,INTENT(IN)     :: lle
+    REAL(rstd),INTENT(INOUT) :: theta_divgrad_i(iim*jjm,1:lle)
+    REAL(rstd)  :: grad_e(3*iim*jjm,llb:lle)
+
+    INTEGER :: ij,l    
+    DO l=llb,lle
+      !DIR$ SIMD
+      DO ij=ij_begin_ext,ij_end_ext 
+          grad_e(ij+u_right,l)=-1/de(ij+u_right)*(ne(ij,right)*theta_divgrad_i(ij,l)+ ne(ij+t_right,left)*theta_divgrad_i(ij+t_right,l) )
+          grad_e(ij+u_lup,l)=-1/de(ij+u_lup)*(ne(ij,lup)*theta_divgrad_i(ij,l)+ ne(ij+t_lup,rdown)*theta_divgrad_i(ij+t_lup,l ))
+          grad_e(ij+u_ldown,l)=-1/de(ij+u_ldown)*(ne(ij,ldown)*theta_divgrad_i(ij,l)+ne(ij+t_ldown,rup)*theta_divgrad_i(ij+t_ldown,l) )
+      ENDDO
+    ENDDO    
+    
+    DO l=llb,lle
+      !DIR$ SIMD
+      DO ij=ij_begin,ij_end
+          theta_divgrad_i(ij,l)=1./Ai(ij)*(ne(ij,right)*grad_e(ij+u_right,l)*le(ij+u_right)  +  &
+                             ne(ij,rup)*grad_e(ij+u_rup,l)*le(ij+u_rup)              +  &
+                             ne(ij,lup)*grad_e(ij+u_lup,l)*le(ij+u_lup)              +  &
+                             ne(ij,left)*grad_e(ij+u_left,l)*le(ij+u_left)           +  &
+                             ne(ij,ldown)*grad_e(ij+u_ldown,l)*le(ij+u_ldown)        +  & 
+                             ne(ij,rdown)*grad_e(ij+u_rdown,l)*le(ij+u_rdown))
+      ENDDO
+    ENDDO
+   
+    DO l=llb,lle
+      DO ij=ij_begin,ij_end
+          theta_divgrad_i(ij,l)=-theta_divgrad_i(ij,l)*cdivgrad
+      ENDDO
+    ENDDO
+  END SUBROUTINE compute_divgrad_inplace
+
+  SUBROUTINE compute_gradrot(ue,gradrot_e,llb,lle)
+    INTEGER,INTENT(IN)     :: llb
+    INTEGER,INTENT(IN)     :: lle
+    REAL(rstd),INTENT(IN) :: ue(iim*3*jjm,lle)
+    REAL(rstd),INTENT(OUT) :: gradrot_e(iim*3*jjm,lle)
+    
+    gradrot_e = ue
+    CALL compute_gradrot_inplace(gradrot_e,llb,lle)
+  END SUBROUTINE
+    
+  SUBROUTINE compute_gradrot_inplace(ue_gradrot_e,llb,lle)
+  USE icosa
+  IMPLICIT NONE
+    INTEGER,INTENT(IN)     :: llb
+    INTEGER,INTENT(IN)     :: lle
+    REAL(rstd),INTENT(INOUT) :: ue_gradrot_e(iim*3*jjm,lle)
+    REAL(rstd) :: rot_v(2*iim*jjm,llb:lle)
+
+    INTEGER :: ij,l
+
+    DO l=llb,lle
+!DIR$ SIMD
+      DO ij=ij_begin_ext,ij_end_ext        
+          rot_v(ij+z_up,l)= 1./Av(ij+z_up)*(  ne(ij,rup)*ue_gradrot_e(ij+u_rup,l)*de(ij+u_rup) &
+                                + ne(ij+t_rup,left)*ue_gradrot_e(ij+t_rup+u_left,l)*de(ij+t_rup+u_left) &
+                                - ne(ij,lup)*ue_gradrot_e(ij+u_lup,l)*de(ij+u_lup) )                               
+          rot_v(ij+z_down,l) = 1./Av(ij+z_down)*( ne(ij,ldown)*ue_gradrot_e(ij+u_ldown,l)*de(ij+u_ldown) &
+                                     + ne(ij+t_ldown,right)*ue_gradrot_e(ij+t_ldown+u_right,l)*de(ij+t_ldown+u_right) &
+                                     - ne(ij,rdown)*ue_gradrot_e(ij+u_rdown,l)*de(ij+u_rdown) )
+      ENDDO
+    ENDDO                              
+  
     DO l=llb,lle
 !DIR$ SIMD
       DO ij=ij_begin,ij_end
-          divu_i(ij,l)=1./Ai(ij)*(ne(ij,right)*ue(ij+u_right,l)*le(ij+u_right)  +  &
-                             ne(ij,rup)*ue(ij+u_rup,l)*le(ij+u_rup)        +  &  
-                             ne(ij,lup)*ue(ij+u_lup,l)*le(ij+u_lup)        +  &  
-                             ne(ij,left)*ue(ij+u_left,l)*le(ij+u_left)     +  &  
-                             ne(ij,ldown)*ue(ij+u_ldown,l)*le(ij+u_ldown)  +  &  
-                             ne(ij,rdown)*ue(ij+u_rdown,l)*le(ij+u_rdown))
-      ENDDO
-    ENDDO
-    
+          ue_gradrot_e(ij+u_right,l)=1/le(ij+u_right)*ne(ij,right)*(rot_v(ij+z_rdown,l)-rot_v(ij+z_rup,l))       
+          ue_gradrot_e(ij+u_lup,l)=1/le(ij+u_lup)*ne(ij,lup)*(rot_v(ij+z_up,l)-rot_v(ij+z_lup,l))        
+          ue_gradrot_e(ij+u_ldown,l)=1/le(ij+u_ldown)*ne(ij,ldown)*(rot_v(ij+z_ldown,l)-rot_v(ij+z_down,l))
+      ENDDO
+    ENDDO
+
     DO l=llb,lle
 !DIR$ SIMD
-      DO ij=ij_begin,ij_end
- 
-          gradivu_e(ij+u_right,l)=-1/de(ij+u_right)*(ne(ij,right)*divu_i(ij,l)+ ne(ij+t_right,left)*divu_i(ij+t_right,l) )       
-
-          gradivu_e(ij+u_lup,l)=-1/de(ij+u_lup)*(ne(ij,lup)*divu_i(ij,l)+ ne(ij+t_lup,rdown)*divu_i(ij+t_lup,l))       
-    
-          gradivu_e(ij+u_ldown,l)=-1/de(ij+u_ldown)*(ne(ij,ldown)*divu_i(ij,l)+ne(ij+t_ldown,rup)*divu_i(ij+t_ldown,l) )
-
-      ENDDO
-    ENDDO
-
-    DO l=llb,lle
-!DIR$ SIMD
-      DO ij=ij_begin,ij_end
-          gradivu_e(ij+u_right,l)=-gradivu_e(ij+u_right,l)*cgraddiv
-          gradivu_e(ij+u_lup,l)=-gradivu_e(ij+u_lup,l)*cgraddiv
-          gradivu_e(ij+u_ldown,l)=-gradivu_e(ij+u_ldown,l)*cgraddiv
-      ENDDO
-    ENDDO
-
-   
-  END SUBROUTINE compute_gradiv
-  
-  SUBROUTINE compute_divgrad(theta,divgrad_i,llb,lle)
-  USE icosa
-  IMPLICIT NONE
-    INTEGER,INTENT(IN)     :: llb
-    INTEGER,INTENT(IN)     :: lle
-    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)
-    REAL(rstd),INTENT(OUT) :: divgrad_i(iim*jjm,llm)
-    REAL(rstd)  :: grad_e(3*iim*jjm,llb:lle)
-
-    INTEGER :: ij,l
-
-       
-    DO l=llb,lle
-!DIR$ SIMD
-      DO ij=ij_begin_ext,ij_end_ext
- 
-          grad_e(ij+u_right,l)=-1/de(ij+u_right)*(ne(ij,right)*theta(ij,l)+ ne(ij+t_right,left)*theta(ij+t_right,l) )        
-
-          grad_e(ij+u_lup,l)=-1/de(ij+u_lup)*(ne(ij,lup)*theta(ij,l)+ ne(ij+t_lup,rdown)*theta(ij+t_lup,l ))       
-    
-          grad_e(ij+u_ldown,l)=-1/de(ij+u_ldown)*(ne(ij,ldown)*theta(ij,l)+ne(ij+t_ldown,rup)*theta(ij+t_ldown,l) )
-
-      ENDDO
-    ENDDO
-    
-    
-    DO l=llb,lle
-!DIR$ SIMD
-      DO ij=ij_begin,ij_end
-
-          divgrad_i(ij,l)=1./Ai(ij)*(ne(ij,right)*grad_e(ij+u_right,l)*le(ij+u_right)  +  &
-                             ne(ij,rup)*grad_e(ij+u_rup,l)*le(ij+u_rup)              +  &  
-                             ne(ij,lup)*grad_e(ij+u_lup,l)*le(ij+u_lup)              +  &  
-                             ne(ij,left)*grad_e(ij+u_left,l)*le(ij+u_left)           +  &  
-                             ne(ij,ldown)*grad_e(ij+u_ldown,l)*le(ij+u_ldown)        +  &  
-                             ne(ij,rdown)*grad_e(ij+u_rdown,l)*le(ij+u_rdown))
-      ENDDO
-    ENDDO
-   
-    DO l=llb,lle
-      DO ij=ij_begin,ij_end
-          divgrad_i(ij,l)=-divgrad_i(ij,l)*cdivgrad
-      ENDDO
-    ENDDO
-
-  END SUBROUTINE compute_divgrad
-
-    
-  SUBROUTINE compute_gradrot(ue,gradrot_e,llb,lle)
-  USE icosa
-  IMPLICIT NONE
-    INTEGER,INTENT(IN)     :: llb
-    INTEGER,INTENT(IN)     :: lle
-    REAL(rstd),INTENT(IN)  :: ue(iim*3*jjm,llm)
-    REAL(rstd),INTENT(OUT) :: gradrot_e(iim*3*jjm,llm)
-    REAL(rstd) :: rot_v(2*iim*jjm,llb:lle)
-
-    INTEGER :: ij,l
-     
-    DO l=llb,lle
-!DIR$ SIMD
-      DO ij=ij_begin_ext,ij_end_ext
-        
-          rot_v(ij+z_up,l)= 1./Av(ij+z_up)*(  ne(ij,rup)*ue(ij+u_rup,l)*de(ij+u_rup)                     &
-                                + ne(ij+t_rup,left)*ue(ij+t_rup+u_left,l)*de(ij+t_rup+u_left)          &
-                                - ne(ij,lup)*ue(ij+u_lup,l)*de(ij+u_lup) ) 
-                              
-          rot_v(ij+z_down,l) = 1./Av(ij+z_down)*( ne(ij,ldown)*ue(ij+u_ldown,l)*de(ij+u_ldown)                 &
-                                     + ne(ij+t_ldown,right)*ue(ij+t_ldown+u_right,l)*de(ij+t_ldown+u_right)  &
-                                     - ne(ij,rdown)*ue(ij+u_rdown,l)*de(ij+u_rdown) )
-  
-      ENDDO
-    ENDDO                              
-  
-    DO l=llb,lle
-!DIR$ SIMD
-      DO ij=ij_begin,ij_end
-  
-          gradrot_e(ij+u_right,l)=1/le(ij+u_right)*ne(ij,right)*(rot_v(ij+z_rdown,l)-rot_v(ij+z_rup,l)) 
-         
-          gradrot_e(ij+u_lup,l)=1/le(ij+u_lup)*ne(ij,lup)*(rot_v(ij+z_up,l)-rot_v(ij+z_lup,l)) 
-        
-          gradrot_e(ij+u_ldown,l)=1/le(ij+u_ldown)*ne(ij,ldown)*(rot_v(ij+z_ldown,l)-rot_v(ij+z_down,l))
-        
-      ENDDO
-    ENDDO
-
-    DO l=llb,lle
-!DIR$ SIMD
-      DO ij=ij_begin,ij_end
-    
-          gradrot_e(ij+u_right,l)=-gradrot_e(ij+u_right,l)*cgradrot       
-          gradrot_e(ij+u_lup,l)=-gradrot_e(ij+u_lup,l)*cgradrot
-          gradrot_e(ij+u_ldown,l)=-gradrot_e(ij+u_ldown,l)*cgradrot
-        
+      DO ij=ij_begin,ij_end    
+          ue_gradrot_e(ij+u_right,l)=-ue_gradrot_e(ij+u_right,l)*cgradrot       
+          ue_gradrot_e(ij+u_lup,l)=-ue_gradrot_e(ij+u_lup,l)*cgradrot
+          ue_gradrot_e(ij+u_ldown,l)=-ue_gradrot_e(ij+u_ldown,l)*cgradrot        
       ENDDO
     ENDDO  
-   
-  END SUBROUTINE compute_gradrot
+  END SUBROUTINE compute_gradrot_inplace
 
 
