Index: codes/icosagcm/trunk/src/caldyn_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 49)
+++ codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 50)
@@ -2,30 +2,26 @@
   USE icosa
 
-  PRIVATE 
-  TYPE(t_field),POINTER :: f_out(:)
-  REAL(rstd),POINTER :: out(:,:)
+  INTEGER :: itau_out
   TYPE(t_field),POINTER :: f_out_u(:)
   REAL(rstd),POINTER :: out_u(:,:)
-  TYPE(t_field),POINTER :: f_out_z(:)
-  REAL(rstd),POINTER :: out_z(:,:)
- 
-  INTEGER :: itau_out
-
-  PUBLIC init_caldyn, caldyn
+
+  TYPE(t_field),POINTER :: f_buf_i(:), f_buf_ulon(:), f_buf_ulat(:), f_buf_u3d(:)
+  TYPE(t_field),POINTER :: f_buf_v(:), f_buf_s(:), f_buf_p(:)
+
+  PUBLIC init_caldyn, caldyn, write_output_fields
 
 CONTAINS
   
   SUBROUTINE init_caldyn(dt)
-  USE icosa
-  IMPLICIT NONE
+    USE icosa
+    IMPLICIT NONE
     REAL(rstd),INTENT(IN) :: dt
     INTEGER :: write_period
-
-    CALL allocate_caldyn
+    
+    write_period=0
     CALL getin('write_period',write_period)
     write_period=write_period/scale_factor
-    
     itau_out=INT(write_period/dt)
-    
+  
     CALL allocate_caldyn
   
@@ -36,8 +32,14 @@
   IMPLICIT NONE
 
-    CALL allocate_field(f_out,field_t,type_real,llm) 
-    CALL allocate_field(f_out_u,field_u,type_real,llm) 
-    CALL allocate_field(f_out_z,field_z,type_real,llm) 
-    
+  CALL allocate_field(f_out_u,field_u,type_real,llm) 
+  
+  CALL allocate_field(f_buf_i,field_t,type_real,llm)
+  CALL allocate_field(f_buf_p,field_t,type_real,llm+1)
+  CALL allocate_field(f_buf_u3d,field_t,type_real,3,llm)  ! 3D vel at cell centers
+  CALL allocate_field(f_buf_ulon,field_t,type_real,llm)
+  CALL allocate_field(f_buf_ulat,field_t,type_real,llm)
+  CALL allocate_field(f_buf_v,field_z,type_real,llm)
+  CALL allocate_field(f_buf_s,field_t,type_real)
+
   END SUBROUTINE allocate_caldyn
   
@@ -45,7 +47,7 @@
   IMPLICIT NONE
     INTEGER,INTENT(IN) :: ind
-    out=f_out(ind) 
+!    out=f_out(ind) 
     out_u=f_out_u(ind) 
-    out_z=f_out_z(ind) 
+!    out_z=f_out_z(ind) 
       
   END SUBROUTINE swap_caldyn
@@ -88,6 +90,4 @@
 
   END SUBROUTINE check_mass_conservation
-  
-  
   
   SUBROUTINE caldyn(it,f_phis, f_ps, f_theta_rhodz, f_u, f_dps, f_dtheta_rhodz, f_du)
@@ -159,5 +159,5 @@
     ENDDO
 
-    CALL transfert_request(f_out_u,req_e1)
+!    CALL transfert_request(f_out_u,req_e1)
 !    CALL transfert_request(f_out_u,req_e1) 
 
@@ -165,24 +165,10 @@
 
     IF (mod(it,itau_out)==0 ) THEN
-      CALL writefield("ps",f_ps)
-      CALL writefield("dps",f_dps)
-!      CALL writefield("theta_rhodz",f_theta_rhodz)
-!    CALL kinetic(f_u,f_out)
-!      CALL writefield("Ki",f_out)
-!      CALL writefield("dtheta_rhodz",f_dtheta_rhodz)
-      CALL vorticity(f_u,f_out_z)
-      CALL writefield("vort",f_out_z)
-!      CALL writefield("theta",f_out)
-      CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_out) ;
-      CALL writefield("T",f_out)
-    
-!      CALL writefield("out",f_out)
-!      DO ind=1,ndomain
-!        CALL writefield("Ki",f_out,ind)
-!        CALL writefield("vort",f_out_z,ind)
-!        CALL writefield("dps",f_dps,ind)
-!      ENDDO
-
-    ENDIF
+       PRINT *,'CALL write_output_fields'
+       CALL write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, &
+            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)
 
@@ -193,4 +179,5 @@
   USE icosa
   USE disvert_mod
+  USE exner_mod
   IMPLICIT NONE
     REAL(rstd),INTENT(IN) :: phis(iim*jjm)
@@ -297,60 +284,7 @@
     ENDDO
       
-
-
- 
- !!! Compute Exnher function
-  
- !! Compute Alpha and Beta
-  
- ! for llm layer
-!$OMP DO
-   DO j=jj_begin-1,jj_end+1
-     DO i=ii_begin-1,ii_end+1
-       ij=(j-1)*iim+i
-       alpha(ij,llm) = 0.
-       beta (ij,llm) = 1./ (1+ 2*kappa)
-     ENDDO
-   ENDDO
-   
- ! for other layer   
-   DO l = llm-1 , 2 , -1
-!$OMP DO
-     DO j=jj_begin-1,jj_end+1
-       DO i=ii_begin-1,ii_end+1
-          ij=(j-1)*iim+i
-          delta = p(ij,l)* (1+2*kappa) + p(ij,l+1)* ( beta(ij,l+1)- (1+2*kappa) )
-          alpha(ij,l)  = - p(ij,l+1) / delta * alpha(ij,l+1)
-          beta (ij,l)  =   p(ij,l  ) / delta   
-       ENDDO
-     ENDDO
-   ENDDO
-   
- !! Compute pk
- 
- ! 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
-        pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
-        pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )     /               &
-                   (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-(1.+2*kappa) )* p(ij,2) )
-     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
-         pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
-       ENDDO
-     ENDDO
-   ENDDO
-   
-  
+ !!! Compute Exner function
+    CALL compute_exner(ps,p,pks,pk,1) 
+
 !!! Compute mass
    DO l = 1, llm
@@ -484,13 +418,13 @@
   
 
-  DO l = 1, llm
-!$OMP DO
-   DO j=jj_begin,jj_end
-    DO i=ii_begin,ii_end
-      ij=(j-1)*iim+i
-      out(ij,l)=theta(ij,l)-288
-    ENDDO
-  ENDDO
- ENDDO
+!  DO l = 1, llm
+!!$OMP DO
+!   DO j=jj_begin,jj_end
+!    DO i=ii_begin,ii_end
+!      ij=(j-1)*iim+i
+!      out(ij,l)=theta(ij,l)-288
+!    ENDDO
+!  ENDDO
+! ENDDO
 
 
@@ -680,4 +614,5 @@
     ENDDO
   ENDDO
+
 !!! contribution due to vertical advection
  
@@ -743,5 +678,85 @@
   END SUBROUTINE compute_caldyn
   
-  
+  SUBROUTINE write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, &
+       f_buf_i, f_buf_v, f_buf_i3, f_buf1_i, f_buf2_i, f_buf_s, f_buf_p)
+    USE icosa
+    USE vorticity_mod
+    USE theta2theta_rhodz_mod
+    USE pression_mod
+    USE write_field
+    TYPE(t_field),POINTER :: f_ps(:), f_phis(:), f_u(:), f_theta_rhodz(:), f_dps(:), &
+         f_buf_i(:), f_buf_v(:), f_buf_i3(:), f_buf1_i(:), f_buf2_i(:), f_buf_s(:), f_buf_p(:)
+    
+    CALL writefield("dps",f_dps)
+    CALL vorticity(f_u,f_buf_v)
+    CALL writefield("vort",f_buf_v)
+    
+    ! Temperature
+    CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i) ;
+    CALL writefield("T",f_buf_i)
+    
+    ! velocity components
+    CALL un2ulonlat(f_u, f_buf_i3, f_buf1_i, f_buf2_i)
+    CALL writefield("ulon",f_buf1_i)
+    CALL writefield("ulat",f_buf2_i)
+    
+    ! geopotential
+    CALL thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_buf_s,f_buf_p,f_buf1_i,f_buf2_i,f_buf_i)
+    CALL writefield("p",f_buf_p)
+    CALL writefield("phi",f_buf_i)
+    CALL writefield("theta",f_buf1_i) ! potential temperature
+    CALL writefield("pk",f_buf2_i) ! Exner pressure
+  
+  END SUBROUTINE write_output_fields
+  
+  SUBROUTINE thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_pks,f_p,f_theta,f_pk,f_phi) 
+    USE field_mod
+    USE pression_mod
+    USE exner_mod
+    USE geopotential_mod
+    USE theta2theta_rhodz_mod
+    TYPE(t_field), POINTER :: f_ps(:), f_phis(:), f_theta_rhodz(:), &  ! IN
+         f_pks(:), f_p(:), f_theta(:), f_pk(:), f_phi(:)               ! OUT
+    REAL(rstd),POINTER :: pk(:,:), p(:,:), theta(:,:), theta_rhodz(:,:), &
+         phi(:,:), phis(:), ps(:), pks(:)
+    INTEGER :: ind
+
+    DO ind=1,ndomain
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       ps = f_ps(ind)
+       p  = f_p(ind)
+       CALL compute_pression(ps,p,0)
+       pk = f_pk(ind)
+       pks = f_pks(ind)
+       CALL compute_exner(ps,p,pks,pk,0)
+       theta_rhodz = f_theta_rhodz(ind)
+       theta = f_theta(ind)
+       CALL compute_theta_rhodz2theta(ps, theta_rhodz,theta,0)
+       phis = f_phis(ind)
+       phi = f_phi(ind)
+       CALL compute_geopotential(phis,pks,pk,theta,phi,0)
+    END DO
+
+  END SUBROUTINE thetarhodz2geopot
+  
+  SUBROUTINE un2ulonlat(f_u, f_u3d, f_ulon, f_ulat)
+    USE field_mod
+    USE wind_mod
+    TYPE(t_field), POINTER :: f_u(:), & ! IN  : normal velocity components on edges
+         f_u3d(:), f_ulon(:), f_ulat(:) ! OUT : velocity reconstructed at hexagons
+    REAL(rstd),POINTER :: u(:,:), u3d(:,:,:), ulon(:,:), ulat(:,:)
+    INTEGER :: ind
+    DO ind=1,ndomain
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       u=f_u(ind)
+       u3d=f_u3d(ind)
+       CALL compute_wind_centered(u,u3d)
+       ulon=f_ulon(ind)
+       ulat=f_ulat(ind)
+       CALL compute_wind_centered_lonlat_compound(u3d, ulon, ulat)
+    END DO
+  END SUBROUTINE un2ulonlat
  
 END MODULE caldyn_gcm_mod
Index: codes/icosagcm/trunk/src/exner.f90
===================================================================
--- codes/icosagcm/trunk/src/exner.f90	(revision 49)
+++ codes/icosagcm/trunk/src/exner.f90	(revision 50)
@@ -43,52 +43,78 @@
     REAL(rstd) :: delta 
     
+    IF(.FALSE.) THEN          ! LMD-Z style calculation of Exner pressure
+       !! Compute Alpha and Beta
+       
+       ! for llm layer
+       !$OMP DO
+       DO j=jj_begin-offset,jj_end+offset
+          DO i=ii_begin-offset,ii_end+offset
+             ij=(j-1)*iim+i
+             alpha(ij,llm) = 0.
+             beta (ij,llm) = 1./ (1+ 2*kappa)
+          ENDDO
+       ENDDO
+       
+       ! for other layer   
+       DO l = llm-1 , 2 , -1
+          !$OMP DO
+          DO j=jj_begin-offset,jj_end+offset
+             DO i=ii_begin-offset,ii_end+offset
+                ij=(j-1)*iim+i
+                delta = p(ij,l)* (1+2*kappa) + p(ij,l+1)* ( beta(ij,l+1)- (1+2*kappa) )
+                alpha(ij,l)  = - p(ij,l+1) / delta * alpha(ij,l+1)
+                beta (ij,l)  =   p(ij,l  ) / delta   
+             ENDDO
+          ENDDO
+       ENDDO
+       
+       !! Compute pk
+       
+       ! for first layer
+       !$OMP DO
+       DO j=jj_begin-offset,jj_end+offset
+          DO i=ii_begin-offset,ii_end+offset
+             ij=(j-1)*iim+i
+             pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
+             pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )     /    &
+                  (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-(1.+2*kappa) )* p(ij,2) )
+          ENDDO
+       ENDDO
+       
+       ! for other layers
+       
+       DO l = 2, llm
+          !$OMP DO
+          DO j=jj_begin-offset,jj_end+offset
+             DO i=ii_begin-offset,ii_end+offset
+                ij=(j-1)*iim+i
+                pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
+             ENDDO
+          ENDDO
+       ENDDO
+       
+    ELSE ! Simple calculation of Exner pressure based on centered average
+       ! surface : pks
+       !$OMP DO
+       DO j=jj_begin-offset,jj_end+offset
+          DO i=ii_begin-offset,ii_end+offset
+             ij=(j-1)*iim+i
+             pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
+          ENDDO
+       ENDDO
 
-   !! Compute Alpha and Beta
-  
-   ! for llm layer
-     DO j=jj_begin-offset,jj_end+offset
-       DO i=ii_begin-offset,ii_end+offset
-         ij=(j-1)*iim+i
-         alpha(ij,llm) = 0.
-         beta (ij,llm) = 1./ (1+ 2*kappa)
+       ! 3D : pk
+       DO l = 1, llm
+          !$OMP DO
+          DO j=jj_begin-offset,jj_end+offset
+             DO i=ii_begin-offset,ii_end+offset
+                ij=(j-1)*iim+i
+                pk(ij,l) = cpp * ((.5/preff)*(p(ij,l)+p(ij,l+1))) ** kappa
+             ENDDO
+          ENDDO
        ENDDO
-     ENDDO
-   
-   ! for other layer   
-     DO l = llm-1 , 2 , -1
-       DO j=jj_begin-offset,jj_end+offset
-         DO i=ii_begin-offset,ii_end+offset
-            ij=(j-1)*iim+i
-            delta = p(ij,l)* (1+2*kappa) + p(ij,l+1)* ( beta(ij,l+1)- (1+2*kappa) )
-            alpha(ij,l)  = - p(ij,l+1) / delta * alpha(ij,l+1)
-            beta (ij,l)  =   p(ij,l  ) / delta   
-         ENDDO
-       ENDDO
-     ENDDO     
-      
-   !! Compute pk
- 
-   ! for first layer
-     DO j=jj_begin-offset,jj_end+offset
-       DO i=ii_begin-offset,ii_end+offset
-          ij=(j-1)*iim+i
-          pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
-          pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )     /               &
-                     (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-(1.+2*kappa) )* p(ij,2) )
-       ENDDO
-     ENDDO
-   
-   ! for other layers
-     
-     DO l = 2, llm
-       DO j=jj_begin-offset,jj_end+offset
-         DO i=ii_begin-offset,ii_end+offset
-           ij=(j-1)*iim+i
-           pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
-         ENDDO
-       ENDDO
-     ENDDO
+    END IF
     
   END SUBROUTINE compute_exner
-
+  
 END MODULE exner_mod
Index: codes/icosagcm/trunk/src/timeloop_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 49)
+++ codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 50)
@@ -1,4 +1,3 @@
 MODULE timeloop_gcm_mod
-
   
 CONTAINS
@@ -8,5 +7,4 @@
   USE dissip_gcm_mod
   USE caldyn_mod
-  USE theta2theta_rhodz_mod
   USE etat0_mod
   USE guided_mod
@@ -33,4 +31,5 @@
   REAL(rstd),POINTER :: du(:,:), dum1(:,:), dum2(:,:)
   REAL(rstd),POINTER :: dtheta_rhodz(:,:),dtheta_rhodzm1(:,:),dtheta_rhodzm2(:,:)
+
   REAL(rstd) :: dt, run_length
   INTEGER :: ind
@@ -39,6 +38,4 @@
   INTEGER :: matsuno_period
   INTEGER :: itaumax
-  INTEGER :: write_period  
-  INTEGER :: itau_out
 
   dt=90.
@@ -54,10 +51,4 @@
   dt=dt/scale_factor
 
-  write_period=0
-  CALL getin('write_period',write_period)
-  write_period=write_period/scale_factor
-  
-  itau_out=INT(write_period/dt)
-  
   scheme='adam_bashforth'
   CALL getin('scheme',scheme)
@@ -101,5 +92,5 @@
   
   CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
-  
+
   DO it=0,itaumax
     PRINT *,"It No :",It,"   t :",dt*It
@@ -242,12 +233,9 @@
         dpsm1(:)=dps(:) ; dum1(:,:)=du(:,:) ; dtheta_rhodzm1(:,:)=dtheta_rhodz(:,:)
 
-      ENDDO
-     
+      ENDDO     
       
     END SUBROUTINE adam_bashforth_scheme
 
-  END SUBROUTINE timeloop
-    
-    
+  END SUBROUTINE timeloop    
   
 END MODULE timeloop_gcm_mod
