Index: /codes/icosagcm/trunk/src/caldyn_gcm.f90
===================================================================
--- /codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 127)
+++ /codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 128)
@@ -31,4 +31,5 @@
        STOP
     END SELECT
+    PRINT *, 'caldyn_conserv=',def
 
     def='direct'
@@ -108,9 +109,22 @@
     SELECT CASE(caldyn_conserv)
     CASE(1) ! energy-conserving
-    CASE(2) ! enstrophy-conserving
        DO ind=1,ndomain
           CALL swap_dimensions(ind)
           CALL swap_geometry(ind)
-          
+          ps=f_ps(ind)
+          rhodz=f_rhodz(ind)
+          p=f_p(ind)
+          qu=f_qu(ind)
+          u=f_u(ind)
+          !$OMP PARALLEL DEFAULT(SHARED)
+          CALL compute_pvort(ps, u, p,rhodz,qu)
+          !$OMP END PARALLEL
+       ENDDO
+
+       CALL transfert_request(f_qu,req_e1)
+
+       DO ind=1,ndomain
+          CALL swap_dimensions(ind)
+          CALL swap_geometry(ind)
           phis=f_phis(ind)
           ps=f_ps(ind)
@@ -125,5 +139,24 @@
           out_u=f_out_u(ind) 
           !$OMP PARALLEL DEFAULT(SHARED)
-!          CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, dps, dtheta_rhodz, du)
+          CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, dps, dtheta_rhodz, du)
+          !$OMP END PARALLEL
+       ENDDO       
+       
+    CASE(2) ! enstrophy-conserving
+       DO ind=1,ndomain
+          CALL swap_dimensions(ind)
+          CALL swap_geometry(ind)
+          phis=f_phis(ind)
+          ps=f_ps(ind)
+          dps=f_dps(ind)
+          theta_rhodz=f_theta_rhodz(ind)
+          dtheta_rhodz=f_dtheta_rhodz(ind)
+          rhodz=f_rhodz(ind)
+          p=f_p(ind)
+          qu=f_qu(ind)
+          u=f_u(ind)
+          du=f_du(ind)
+          out_u=f_out_u(ind) 
+          !$OMP PARALLEL DEFAULT(SHARED)
           CALL compute_pvort(ps, u, p,rhodz,qu)
           CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, dps, dtheta_rhodz, du)
@@ -131,19 +164,18 @@
        ENDDO
        
-       IF (mod(it,itau_out)==0 ) THEN
-          PRINT *,'CALL write_output_fields'
-          CALL write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, &
-               f_buf_i, f_buf_v, f_buf_u3d, f_buf_ulon, f_buf_ulat, f_buf_s, f_buf_p)
-       END IF
-    
     CASE DEFAULT
        STOP
     END SELECT
 
+    IF (mod(it,itau_out)==0 ) THEN
+       PRINT *,'CALL write_output_fields'
+       CALL write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, &
+            f_buf_i, f_buf_v, f_buf_u3d, f_buf_ulon, f_buf_ulat, f_buf_s, f_buf_p)
+    END IF
+    
     !    CALL check_mass_conservation(f_ps,f_dps)
     
 END SUBROUTINE caldyn
-  
-  
+    
 SUBROUTINE compute_pvort(ps, u, p,rhodz,qu)
   USE icosa
@@ -151,53 +183,51 @@
   USE exner_mod
   IMPLICIT NONE
-    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
-    REAL(rstd),INTENT(IN)  :: ps(iim*jjm)
-    REAL(rstd),INTENT(OUT) :: p(iim*jjm,llm+1)
-    REAL(rstd),INTENT(OUT) :: rhodz(iim*jjm,llm)
-    REAL(rstd),INTENT(OUT) :: qu(iim*3*jjm,llm)
-    
-    INTEGER :: i,j,ij,l
-!    REAL(rstd) :: ww,uu 
-!    REAL(rstd) :: delta
-    REAL(rstd) :: etav,hv
-    REAL(rstd),ALLOCATABLE,SAVE :: qv(:,:)     ! potential velocity  
-    
-    LOGICAL,SAVE :: first=.TRUE.
-!$OMP THREADPRIVATE(first)
-        
-!$OMP BARRIER      
-!$OMP MASTER  
-!    IF (first) THEN
-    ALLOCATE(qv(2*iim*jjm,llm))     ! potential velocity  
-    
+  REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
+  REAL(rstd),INTENT(IN)  :: ps(iim*jjm)
+  REAL(rstd),INTENT(OUT) :: p(iim*jjm,llm+1)
+  REAL(rstd),INTENT(OUT) :: rhodz(iim*jjm,llm)
+  REAL(rstd),INTENT(OUT) :: qu(iim*3*jjm,llm)
+  
+  INTEGER :: i,j,ij,l
+  REAL(rstd) :: etav,hv
+  REAL(rstd),ALLOCATABLE,SAVE :: qv(:,:)     ! potential velocity  
+  
+  LOGICAL,SAVE :: first=.TRUE.
+  !$OMP THREADPRIVATE(first)
+  
+  !$OMP BARRIER      
+  !$OMP MASTER  
+  !    IF (first) THEN
+  ALLOCATE(qv(2*iim*jjm,llm))     ! potential velocity  
+  
 !!! Compute pressure
-    DO    l    = 1, llm+1
-       !$OMP DO
-       DO j=jj_begin-1,jj_end+1
-          DO i=ii_begin-1,ii_end+1
-             ij=(j-1)*iim+i
-             p(ij,l) = ap(l) + bp(l) * ps(ij)
-          ENDDO
-       ENDDO
-    ENDDO
-    
+  DO    l    = 1, llm+1
+     !$OMP DO
+     DO j=jj_begin-1,jj_end+1
+        DO i=ii_begin-1,ii_end+1
+           ij=(j-1)*iim+i
+           p(ij,l) = ap(l) + bp(l) * ps(ij)
+        ENDDO
+     ENDDO
+  ENDDO
+  
 !!! Compute mass
-   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
-            rhodz(ij,l) = ( p(ij,l) - p(ij,l+1) )/g
-         ENDDO
-      ENDDO
-   ENDDO
-
+  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
+           rhodz(ij,l) = ( p(ij,l) - p(ij,l+1) )/g
+        ENDDO
+     ENDDO
+  ENDDO
+  
 !!! Compute shallow-water potential vorticity
   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
-        
+       DO i=ii_begin-1,ii_end+1
+          ij=(j-1)*iim+i
+          
           etav= 1./Av(ij+z_up)*(  ne(ij,rup)        * u(ij+u_rup,l)        * de(ij+u_rup)         &
                                 + ne(ij+t_rup,left) * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left)  &
@@ -236,23 +266,9 @@
 !!$OMP BARRIER
 !!$OMP MASTER  
-    DEALLOCATE(theta)  ! potential temperature
-    DEALLOCATE(pk)   ! Exner function
-    DEALLOCATE(pks)
-    DEALLOCATE(alpha)
-    DEALLOCATE(beta)
-    DEALLOCATE(phi)    ! geopotential
-    DEALLOCATE(Fe)   ! mass flux   
-    DEALLOCATE(Ftheta) ! theta flux   
-    DEALLOCATE(convm)    ! mass flux convergence
-    DEALLOCATE(w)       ! vertical velocity      
     DEALLOCATE(qv)       ! potential velocity  
-    DEALLOCATE(berni)   ! bernouilli term
 !!$OMP END MASTER
 !!$OMP BARRIER                                                      
   END SUBROUTINE compute_pvort
   
-
-!----------- needed next : rhodz, p --------------!
-
   SUBROUTINE compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, dps, dtheta_rhodz, du)
   USE icosa
@@ -284,5 +300,4 @@
     REAL(rstd),ALLOCATABLE,SAVE :: convm(:,:)  ! mass flux convergence
     REAL(rstd),ALLOCATABLE,SAVE :: w(:,:)      ! vertical velocity      
-    REAL(rstd),ALLOCATABLE,SAVE :: qv(:,:)     ! potential velocity  
     REAL(rstd),ALLOCATABLE,SAVE :: berni(:,:)  ! Bernouilli function
     
@@ -303,5 +318,4 @@
     ALLOCATE(convm(iim*jjm,llm))    ! mass flux convergence
     ALLOCATE(w(iim*jjm,llm))        ! vertical velocity      
-    ALLOCATE(qv(2*iim*jjm,llm))     ! potential velocity  
     ALLOCATE(berni(iim*jjm,llm))    ! bernouilli term
 
@@ -394,52 +408,107 @@
 
 !!! Compute potential vorticity (Coriolis) contribution to du
-  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) = qu(ij+u_right,l)/de(ij+u_right) *      &  
-                      ( wee(ij+u_right,1,1)*Fe(ij+u_rup,l)+                           &
-                        wee(ij+u_right,2,1)*Fe(ij+u_lup,l)+                           &
-                        wee(ij+u_right,3,1)*Fe(ij+u_left,l)+                          &
-                        wee(ij+u_right,4,1)*Fe(ij+u_ldown,l)+                         &
-                        wee(ij+u_right,5,1)*Fe(ij+u_rdown,l)+                         & 
-                        wee(ij+u_right,1,2)*Fe(ij+t_right+u_ldown,l)+                 &
-                        wee(ij+u_right,2,2)*Fe(ij+t_right+u_rdown,l)+                 &
-                        wee(ij+u_right,3,2)*Fe(ij+t_right+u_right,l)+                 &
-                        wee(ij+u_right,4,2)*Fe(ij+t_right+u_rup,l)+                   &
-                        wee(ij+u_right,5,2)*Fe(ij+t_right+u_lup,l) )
-
-
-       du(ij+u_lup,l) = qu(ij+u_lup,l)/de(ij+u_lup) *         &  
-                      ( wee(ij+u_lup,1,1)*Fe(ij+u_left,l)+                        &
-                        wee(ij+u_lup,2,1)*Fe(ij+u_ldown,l)+                       &
-                        wee(ij+u_lup,3,1)*Fe(ij+u_rdown,l)+                       &
-                        wee(ij+u_lup,4,1)*Fe(ij+u_right,l)+                       &
-                        wee(ij+u_lup,5,1)*Fe(ij+u_rup,l)+                         & 
-                        wee(ij+u_lup,1,2)*Fe(ij+t_lup+u_right,l)+                 &
-                        wee(ij+u_lup,2,2)*Fe(ij+t_lup+u_rup,l)+                   &
-                        wee(ij+u_lup,3,2)*Fe(ij+t_lup+u_lup,l)+                   &
-                        wee(ij+u_lup,4,2)*Fe(ij+t_lup+u_left,l)+                  &
-                        wee(ij+u_lup,5,2)*Fe(ij+t_lup+u_ldown,l) )
-
-
-       du(ij+u_ldown,l) = qu(ij+u_ldown,l)/de(ij+u_ldown) *   &  
-                      ( wee(ij+u_ldown,1,1)*Fe(ij+u_rdown,l)+                      &
-                        wee(ij+u_ldown,2,1)*Fe(ij+u_right,l)+                      &
-                        wee(ij+u_ldown,3,1)*Fe(ij+u_rup,l)+                        &
-                        wee(ij+u_ldown,4,1)*Fe(ij+u_lup,l)+                        &
-                        wee(ij+u_ldown,5,1)*Fe(ij+u_left,l)+                       & 
-                        wee(ij+u_ldown,1,2)*Fe(ij+t_ldown+u_lup,l)+                &
-                        wee(ij+u_ldown,2,2)*Fe(ij+t_ldown+u_left,l)+               &
-                        wee(ij+u_ldown,3,2)*Fe(ij+t_ldown+u_ldown,l)+              &
-                        wee(ij+u_ldown,4,2)*Fe(ij+t_ldown+u_rdown,l)+              &
-                        wee(ij+u_ldown,5,2)*Fe(ij+t_ldown+u_right,l) )
-
+
+    SELECT CASE(caldyn_conserv)
+    CASE(1) ! energy-conserving TRiSK
+
+       DO l=1,llm
+          !$OMP DO
+          DO j=jj_begin,jj_end
+             DO i=ii_begin,ii_end
+                ij=(j-1)*iim+i
+
+                uu = wee(ij+u_right,1,1)*Fe(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))+                             &
+                     wee(ij+u_right,2,1)*Fe(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))+                             &
+                     wee(ij+u_right,3,1)*Fe(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+                           &
+                     wee(ij+u_right,4,1)*Fe(ij+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+u_ldown,l))+                         &
+                     wee(ij+u_right,5,1)*Fe(ij+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+u_rdown,l))+                         & 
+                     wee(ij+u_right,1,2)*Fe(ij+t_right+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_ldown,l))+         &
+                     wee(ij+u_right,2,2)*Fe(ij+t_right+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rdown,l))+         &
+                     wee(ij+u_right,3,2)*Fe(ij+t_right+u_right,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_right,l))+         &
+                     wee(ij+u_right,4,2)*Fe(ij+t_right+u_rup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rup,l))+             &
+                     wee(ij+u_right,5,2)*Fe(ij+t_right+u_lup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_lup,l))
+                du(ij+u_right,l) = .5*uu/de(ij+u_right)
+                
+                uu = wee(ij+u_lup,1,1)*Fe(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) +                        &
+                     wee(ij+u_lup,2,1)*Fe(ij+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) +                      &
+                     wee(ij+u_lup,3,1)*Fe(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) +                      &
+                     wee(ij+u_lup,4,1)*Fe(ij+u_right,l)*(qu(ij+u_lup,l)+qu(ij+u_right,l)) +                      &
+                     wee(ij+u_lup,5,1)*Fe(ij+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+u_rup,l)) +                          & 
+                     wee(ij+u_lup,1,2)*Fe(ij+t_lup+u_right,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_right,l)) +          &
+                     wee(ij+u_lup,2,2)*Fe(ij+t_lup+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_rup,l)) +              &
+                     wee(ij+u_lup,3,2)*Fe(ij+t_lup+u_lup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_lup,l)) +              &
+                     wee(ij+u_lup,4,2)*Fe(ij+t_lup+u_left,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_left,l)) +            &
+                     wee(ij+u_lup,5,2)*Fe(ij+t_lup+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_ldown,l))
+                du(ij+u_lup,l) = .5*uu/de(ij+u_lup)
+
+                
+                uu = wee(ij+u_ldown,1,1)*Fe(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) +                      &
+                     wee(ij+u_ldown,2,1)*Fe(ij+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+u_right,l)) +                      &
+                     wee(ij+u_ldown,3,1)*Fe(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) +                          &
+                     wee(ij+u_ldown,4,1)*Fe(ij+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) +                          &
+                     wee(ij+u_ldown,5,1)*Fe(ij+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+u_left,l)) +                        & 
+                     wee(ij+u_ldown,1,2)*Fe(ij+t_ldown+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_lup,l)) +          &
+                     wee(ij+u_ldown,2,2)*Fe(ij+t_ldown+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_left,l)) +        &
+                     wee(ij+u_ldown,3,2)*Fe(ij+t_ldown+u_ldown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_ldown,l)) +      &
+                     wee(ij+u_ldown,4,2)*Fe(ij+t_ldown+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_rdown,l)) +      &
+                     wee(ij+u_ldown,5,2)*Fe(ij+t_ldown+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_right,l))
+                du(ij+u_ldown,l) = .5*uu/de(ij+u_ldown)
+                
+             ENDDO
+          ENDDO
+       ENDDO
+
+    CASE(2) ! enstrophy-conserving TRiSK
+  
+       DO l=1,llm
+          !$OMP DO
+          DO j=jj_begin,jj_end
+             DO i=ii_begin,ii_end
+                ij=(j-1)*iim+i
+
+                uu = wee(ij+u_right,1,1)*Fe(ij+u_rup,l)+                           &
+                     wee(ij+u_right,2,1)*Fe(ij+u_lup,l)+                           &
+                     wee(ij+u_right,3,1)*Fe(ij+u_left,l)+                          &
+                     wee(ij+u_right,4,1)*Fe(ij+u_ldown,l)+                         &
+                     wee(ij+u_right,5,1)*Fe(ij+u_rdown,l)+                         & 
+                     wee(ij+u_right,1,2)*Fe(ij+t_right+u_ldown,l)+                 &
+                     wee(ij+u_right,2,2)*Fe(ij+t_right+u_rdown,l)+                 &
+                     wee(ij+u_right,3,2)*Fe(ij+t_right+u_right,l)+                 &
+                     wee(ij+u_right,4,2)*Fe(ij+t_right+u_rup,l)+                   &
+                     wee(ij+u_right,5,2)*Fe(ij+t_right+u_lup,l)
+                du(ij+u_right,l) = qu(ij+u_right,l)*uu/de(ij+u_right)
+                
+                
+                uu = wee(ij+u_lup,1,1)*Fe(ij+u_left,l)+                        &
+                     wee(ij+u_lup,2,1)*Fe(ij+u_ldown,l)+                       &
+                     wee(ij+u_lup,3,1)*Fe(ij+u_rdown,l)+                       &
+                     wee(ij+u_lup,4,1)*Fe(ij+u_right,l)+                       &
+                     wee(ij+u_lup,5,1)*Fe(ij+u_rup,l)+                         & 
+                     wee(ij+u_lup,1,2)*Fe(ij+t_lup+u_right,l)+                 &
+                     wee(ij+u_lup,2,2)*Fe(ij+t_lup+u_rup,l)+                   &
+                     wee(ij+u_lup,3,2)*Fe(ij+t_lup+u_lup,l)+                   &
+                     wee(ij+u_lup,4,2)*Fe(ij+t_lup+u_left,l)+                  &
+                     wee(ij+u_lup,5,2)*Fe(ij+t_lup+u_ldown,l)
+                du(ij+u_lup,l) = qu(ij+u_lup,l)*uu/de(ij+u_lup)
+
+                uu = wee(ij+u_ldown,1,1)*Fe(ij+u_rdown,l)+                      &
+                     wee(ij+u_ldown,2,1)*Fe(ij+u_right,l)+                      &
+                     wee(ij+u_ldown,3,1)*Fe(ij+u_rup,l)+                        &
+                     wee(ij+u_ldown,4,1)*Fe(ij+u_lup,l)+                        &
+                     wee(ij+u_ldown,5,1)*Fe(ij+u_left,l)+                       & 
+                     wee(ij+u_ldown,1,2)*Fe(ij+t_ldown+u_lup,l)+                &
+                     wee(ij+u_ldown,2,2)*Fe(ij+t_ldown+u_left,l)+               &
+                     wee(ij+u_ldown,3,2)*Fe(ij+t_ldown+u_ldown,l)+              &
+                     wee(ij+u_ldown,4,2)*Fe(ij+t_ldown+u_rdown,l)+              &
+                     wee(ij+u_ldown,5,2)*Fe(ij+t_ldown+u_right,l)
+                du(ij+u_ldown,l) = qu(ij+u_ldown,l)*uu/de(ij+u_ldown)
      
-      ENDDO
-    ENDDO
-  ENDDO
+             ENDDO
+          ENDDO
+       ENDDO
+       
+    CASE DEFAULT
+       STOP
+    END SELECT
   
 !!! Compute Exner function
@@ -567,5 +636,4 @@
     DEALLOCATE(convm)    ! mass flux convergence
     DEALLOCATE(w)       ! vertical velocity      
-    DEALLOCATE(qv)       ! potential velocity  
     DEALLOCATE(berni)   ! bernouilli term
 !!$OMP END MASTER
