Index: codes/icosagcm/trunk/src/caldyn_hevi.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_hevi.f90	(revision 376)
+++ codes/icosagcm/trunk/src/caldyn_hevi.f90	(revision 377)
@@ -153,4 +153,5 @@
        dtheta_rhodz=f_dtheta_rhodz(ind)
        du=f_du_slow(ind)
+
        IF(hydrostatic) THEN
           CALL compute_caldyn_slow_hydro(u,mass,hflux,du)
Index: codes/icosagcm/trunk/src/caldyn_kernels_base.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_kernels_base.f90	(revision 376)
+++ codes/icosagcm/trunk/src/caldyn_kernels_base.f90	(revision 377)
@@ -302,9 +302,15 @@
        !DIR$ SIMD
        DO ij=ij_begin,ij_end
-          dW(ij,l) = dW(ij,l) + W_etadot(ij,l-1) - W_etadot(ij,l)
           dPhi(ij,l) = dPhi(ij,l) - wflux(ij,l) &
                * (geopot(ij,l+1)-geopot(ij,l-1))/(mass(ij,l-1)+mass(ij,l))
        END DO
     END DO
+    DO l=ll_begin,ll_end
+       !DIR$ SIMD
+       DO ij=ij_begin,ij_end
+          dW(ij,l+1) = dW(ij,l+1) + W_etadot(ij,l) ! update inner+top interfaces
+          dW(ij,l)   = dW(ij,l)   - W_etadot(ij,l) ! update bottom+inner interfaces
+       END DO
+    END DO
     CALL trace_end("compute_caldyn_vert_nh")
 
Index: codes/icosagcm/trunk/src/caldyn_kernels_hevi.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_kernels_hevi.f90	(revision 376)
+++ codes/icosagcm/trunk/src/caldyn_kernels_hevi.f90	(revision 377)
@@ -155,5 +155,5 @@
 
     ! Newton-Raphson iteration : Phi_il contains current guess value
-    DO iter=1,2 ! 2 iterations should be enough
+    DO iter=1,5 ! 2 iterations should be enough
 
        ! Compute pressure, A_ik
@@ -673,7 +673,8 @@
           kup=l
        END IF
+       ! below : "checked" means "formula also valid when kup=kdown (top/bottom)"
        ! compute mil, wil=Wil/mil
        DO ij=ij_begin_ext, ij_end_ext
-          w_il(ij,l) = 2.*W(ij,l)/(rhodz(ij,kdown)+rhodz(ij,kup))
+          w_il(ij,l) = 2.*W(ij,l)/(rhodz(ij,kdown)+rhodz(ij,kup)) ! checked
        END DO
        ! compute DePhi, v_el, G_el, F_el
@@ -687,5 +688,5 @@
           W2_el = .5*le_de(ij+u_right) * &
                ( W(ij,l)*w_il(ij,l) + W(ij+t_right,l)*w_il(ij+t_right,l) )
-          v_el(ij+u_right) = .5*le_de(ij+u_right)*(u(ij+u_right,kup)+u(ij+u_right,kdown))
+          v_el(ij+u_right) = .5*le_de(ij+u_right)*(u(ij+u_right,kup)+u(ij+u_right,kdown)) ! checked
           G_el(ij+u_right) = v_el(ij+u_right)*W_el - DePhil(ij+u_right,l)*W2_el
           ! Compute on edge 'lup'
@@ -695,5 +696,5 @@
           W2_el = .5*le_de(ij+u_lup) * &
                ( W(ij,l)*w_il(ij,l) + W(ij+t_lup,l)*w_il(ij+t_lup,l) )
-          v_el(ij+u_lup) = .5*le_de(ij+u_lup)*( u(ij+u_lup,kup) + u(ij+u_lup,kdown) )
+          v_el(ij+u_lup) = .5*le_de(ij+u_lup)*( u(ij+u_lup,kup) + u(ij+u_lup,kdown)) ! checked
           G_el(ij+u_lup) = v_el(ij+u_lup)*W_el - DePhil(ij+u_lup,l)*W2_el
           ! Compute on edge 'ldown'
@@ -703,5 +704,5 @@
           W2_el = .5*le_de(ij+u_ldown) * &
                ( W(ij,l)*w_il(ij,l) + W(ij+t_ldown,l)*w_il(ij+t_ldown,l) )
-          v_el(ij+u_ldown) = .5*le_de(ij+u_ldown)*( u(ij+u_ldown,kup) + u(ij+u_ldown,kdown) )
+          v_el(ij+u_ldown) = .5*le_de(ij+u_ldown)*( u(ij+u_ldown,kup) + u(ij+u_ldown,kdown)) ! checked
           G_el(ij+u_ldown) = v_el(ij+u_ldown)*W_el - DePhil(ij+u_ldown,l)*W2_el
        END DO
@@ -715,5 +716,7 @@
                le_de(ij+u_ldown)*DePhil(ij+u_ldown,l)**2 +  &
                le_de(ij+u_rdown)*DePhil(ij+u_rdown,l)**2 )
-          dPhi(ij,l) = gradPhi2(ij,l)*w_il(ij,l)-1/(2*Ai(ij))* & 
+!          gradPhi2(ij,l) = 0. ! FIXME !!
+
+          dPhi(ij,l) = gradPhi2(ij,l)*w_il(ij,l) -1/(2*Ai(ij))* & 
                ( DePhil(ij+u_right,l)*v_el(ij+u_right) + & ! -v.gradPhi, 
                  DePhil(ij+u_rup,l)*v_el(ij+u_rup) +     & ! v_el already has le_de
@@ -722,4 +725,5 @@
                  DePhil(ij+u_ldown,l)*v_el(ij+u_ldown) + &
                  DePhil(ij+u_rdown,l)*v_el(ij+u_rdown) )
+
           dW(ij,l) = -1./Ai(ij)*(           & ! -div(G_el), 
                ne_right*G_el(ij+u_right) +  & ! G_el already has le_de
@@ -731,5 +735,9 @@
        END DO
     END DO
-    
+    ! FIXME !!
+ !   F_el(:,:)=0. 
+ !   dPhi(:,:)=0.
+ !   dW(:,:)=0.
+
     DO l=ll_begin, ll_end ! compute on k levels (layers)
        ! Compute berni at scalar points
@@ -765,4 +773,7 @@
        END DO
     END DO
+    ! FIXME !!
+    ! du(:,:)=0.
+    ! hflux(:,:)=0.
 
     CALL trace_end("compute_caldyn_slow_NH")
Index: codes/icosagcm/trunk/src/dcmip_initial_conditions_test_1_2_3_v5.f90
===================================================================
--- codes/icosagcm/trunk/src/dcmip_initial_conditions_test_1_2_3_v5.f90	(revision 376)
+++ codes/icosagcm/trunk/src/dcmip_initial_conditions_test_1_2_3_v5.f90	(revision 377)
@@ -1196,5 +1196,5 @@
 ! Test 3-1
 !==========
-SUBROUTINE test3_gravity_wave (lon,lat,p,z,zcoords,u,v,w,t,phis,ps,rho,q)
+SUBROUTINE test3_gravity_wave (X,lon,lat,p,z,zcoords,u,v,w,t,phis,ps,rho,q)
 
 IMPLICIT NONE
@@ -1204,5 +1204,6 @@
 
 	real(rstd), intent(in)  :: lon, &	! Longitude (radians)
-                                   lat		! Latitude (radians)
+                                   lat, &	! Latitude (radians)
+                                   X            ! Reduced Earth reduction factor (DCMIP value = 125)
 
 	real(rstd), intent(inout) :: p, &	! Pressure  (Pa)
@@ -1227,7 +1228,6 @@
 !     test case parameters
 !----------------------------------------------------------------------- 
-	real(rstd), parameter :: 	X       = 125.d0,		&	! Reduced Earth reduction factor
+	real(rstd), parameter :: 	&
 			    	Om      = 0.d0,			&	! Rotation Rate of Earth
-                            	as      = a/X,			&	! New Radius of small Earth     
 			    	u0      = 20.d0,		&	! Reference Velocity 
 !			    	u0      = 0.d0,		        &	! FIXME : no zonal wind for NH tests
@@ -1243,5 +1243,7 @@
 				N2      = N*N, 		 	&	! Brunt-Vaisala frequency Squared
 				bigG    = (g*g)/(N2*cp)			! Constant
-                            
+
+      real(rstd) :: as                                                  ! New Radius of small Earth     
+
       real(rstd) :: height							! Model level height
       real(rstd) :: sin_tmp, cos_tmp					! Calculation of great circle distance
@@ -1251,4 +1253,6 @@
       real(rstd) :: theta_pert						! Pot-temp perturbation
 
+      as = a/X
+
 !-----------------------------------------------------------------------
 !    THE VELOCITIES
Index: codes/icosagcm/trunk/src/disvert_std.f90
===================================================================
--- codes/icosagcm/trunk/src/disvert_std.f90	(revision 376)
+++ codes/icosagcm/trunk/src/disvert_std.f90	(revision 377)
@@ -1,4 +1,7 @@
 MODULE disvert_std_mod
   USE icosa
+  IMPLICIT NONE
+  PRIVATE
+
   REAL(rstd), SAVE, ALLOCATABLE,TARGET :: ap(:)
 !$OMP THREADPRIVATE(ap)
@@ -7,4 +10,6 @@
   REAL(rstd), SAVE, ALLOCATABLE,TARGET :: presnivs(:)
 !$OMP THREADPRIVATE(presnivs)
+
+  PUBLIC :: init_disvert, ap, bp, presnivs
 
 CONTAINS
Index: codes/icosagcm/trunk/src/earth_const.f90
===================================================================
--- codes/icosagcm/trunk/src/earth_const.f90	(revision 376)
+++ codes/icosagcm/trunk/src/earth_const.f90	(revision 377)
@@ -17,4 +17,5 @@
 
   LOGICAL, SAVE :: boussinesq, hydrostatic
+  !$OMP THREADPRIVATE(boussinesq, hydrostatic)
 
 CONTAINS
Index: codes/icosagcm/trunk/src/etat0.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0.f90	(revision 376)
+++ codes/icosagcm/trunk/src/etat0.f90	(revision 377)
@@ -21,4 +21,5 @@
     USE etat0_dcmip4_mod, ONLY : getin_etat0_dcmip4=>getin_etat0
     USE etat0_dcmip5_mod, ONLY : getin_etat0_dcmip5=>getin_etat0
+    USE etat0_bubble_mod, ONLY : getin_etat0_bubble=>getin_etat0
     USE etat0_williamson_mod, ONLY : getin_etat0_williamson=>getin_etat0
     USE etat0_temperature_mod, ONLY: getin_etat0_temperature=>getin_etat0
@@ -71,4 +72,6 @@
     CASE ('dcmip5')
         CALL getin_etat0_dcmip5
+    CASE ('bubble')
+        CALL getin_etat0_bubble
     CASE ('williamson91.6')
        autoinit_mass=.FALSE.
@@ -96,5 +99,6 @@
       ELSE
          PRINT*, 'Bad selector for variable etat0 <',etat0_type, &
-              '> options are <jablonowsky06>, <academic>, <dcmip[1-4]> '
+              '> options are <isothermal>, <temperature_profile>, <held_suarez>, &
+              &<bubble>, <venus>, <jablonowsky06>, <academic>, <dcmip[1-4]> '
          STOP
       END IF
@@ -176,4 +180,5 @@
     USE etat0_dcmip4_mod, ONLY : compute_dcmip4 => compute_etat0
     USE etat0_dcmip5_mod, ONLY : compute_dcmip5 => compute_etat0
+    USE etat0_bubble_mod, ONLY : compute_bubble => compute_etat0  
     USE etat0_williamson_mod, ONLY : compute_w91_6 => compute_etat0
     USE etat0_temperature_mod, ONLY: compute_etat0_temperature => compute_etat0
@@ -205,5 +210,5 @@
 
     ! For NH geopotential and vertical momentum must be initialized.
-    ! Unless autoinit_mass is set to .FALSE. , they will be initialized
+    ! Unless autoinit_NH is set to .FALSE. , they will be initialized
     ! to hydrostatic geopotential and zero
     autoinit_mass = .TRUE.
@@ -239,4 +244,8 @@
        CALL compute_dcmip5(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q)
        CALL compute_dcmip5(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
+    CASE('bubble')
+       CALL compute_bubble(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, geopot, q)
+       CALL compute_bubble(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, geopot_e, q_e)
+!       autoinit_NH = .FALSE. ! compute_bubble initializes geopot
     CASE('williamson91.6')
        CALL compute_w91_6(iim*jjm,lon_i,lat_i, phis, mass(:,1), temp_i(:,1), ulon_i(:,1), ulat_i(:,1))
Index: codes/icosagcm/trunk/src/etat0_bubble.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0_bubble.f90	(revision 377)
+++ codes/icosagcm/trunk/src/etat0_bubble.f90	(revision 377)
@@ -0,0 +1,90 @@
+MODULE etat0_bubble_mod
+  USE icosa
+  IMPLICIT NONE
+  PRIVATE
+
+  REAL(rstd), PARAMETER :: latc=0., lonc=0., &
+       width        = 0.3, & ! width of mountain (radian)
+       width_bubble = 0.5, & ! width of bubble (radian), curvature effects O(width^2)
+       zc           = 350., & ! height at bubble center
+       zw           = 250.    ! bubble vertical scale => radius=4000m for circular bubble 
+  REAL(rstd), SAVE :: T0 ! adiabatic atmosphere at theta=T0 (default 300K)
+  !$OMP THREADPRIVATE(T0)
+  REAL(rstd), SAVE :: ztop ! mountain top (default 0.)
+  !$OMP THREADPRIVATE(ztop)
+  REAL(rstd), SAVE :: delta_T ! bubble extra temperature (default 5K)
+  !$OMP THREADPRIVATE(delta_T)
+  
+  PUBLIC getin_etat0, compute_etat0
+
+CONTAINS
+
+  SUBROUTINE getin_etat0
+    T0=300.
+    CALL getin("bubble_T0",T0)
+    ztop=0.
+    CALL getin("bubble_mountain_top",ztop)
+    delta_T=5.
+    CALL getin("bubble_delta_T",delta_T)
+  END SUBROUTINE getin_etat0
+
+  SUBROUTINE compute_etat0(ngrid,lon,lat, phis,ps,temp,ulon,ulat,geopot,q)
+    USE icosa
+    USE disvert_mod
+    USE omp_para
+    INTEGER, INTENT(IN) :: ngrid
+    REAL(rstd),INTENT(IN) :: lon(ngrid)
+    REAL(rstd),INTENT(IN) :: lat(ngrid)
+    REAL(rstd),INTENT(OUT) :: phis(ngrid)
+    REAL(rstd),INTENT(OUT) :: ps(ngrid)
+    REAL(rstd),INTENT(OUT) :: temp(ngrid,llm)
+    REAL(rstd),INTENT(OUT) :: ulon(ngrid,llm)
+    REAL(rstd),INTENT(OUT) :: ulat(ngrid,llm)
+    REAL(rstd),INTENT(OUT) :: geopot(ngrid,llm+1)
+    REAL(rstd),INTENT(OUT) :: q(ngrid,llm,nqtot)
+    
+    INTEGER :: l,ij
+    REAL(rstd) :: sinlat, coslat, K, Ts, p_ij, T_ij, z_ij
+    ! adiabatic profile
+    ! p=preff, T=T0 at Phi=0
+    ! => theta = T0, cp.T + Phi = cst = cp.T0
+    DO ij=1,ngrid
+       sinlat=SIN(lat(ij))
+       coslat=COS(lat(ij))
+       K=sin(latc)*sinlat+cos(latc)*coslat*cos(lon(ij)-lonc)
+       K=acos(K)/width
+       phis(ij) = g*ztop*exp(-(K**2))
+       Ts = T0-phis(ij)/cpp ! adiabatic 
+       ps(ij) = preff * (Ts/T0)**(1/kappa) ! theta=T0=T.(p/preff)^-kappa
+    ENDDO
+    
+    DO l=1,llm+1
+       DO ij=1,ngrid
+          p_ij = ap(l)+bp(l)*ps(ij)
+          T_ij = T0*(p_ij/preff)**kappa ! adiabatic 
+          geopot(ij,l) = cpp*(T0-T_ij)
+       END DO
+    END DO
+    
+    DO l=1, llm
+       DO ij=1,ngrid
+          sinlat=SIN(lat(ij))
+          coslat=COS(lat(ij))
+          K=sin(latc)*sinlat+cos(latc)*coslat
+          K=acos(K)/width_bubble
+          p_ij = 0.5 *( ap(l)+bp(l)*ps(ij) + ap(l+1)+bp(l+1)*ps(ij) )
+          T_ij = T0*(p_ij/preff)**kappa
+          z_ij = (T0-T_ij)*cpp/g
+          z_ij = (z_ij - zc)/zw
+          K=MIN(1.,sqrt(K**2+z_ij**2))
+          K=.5*(1.+COS(pi*K))
+          Temp(ij,l) = T_ij + delta_T*K
+          ulat(ij,l)=0.
+          ulon(ij,l)=0.
+          q(ij,l,:)=0.
+       END DO
+    END DO
+    
+  END SUBROUTINE compute_etat0
+  
+END MODULE etat0_bubble_mod
Index: codes/icosagcm/trunk/src/etat0_dcmip3.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0_dcmip3.f90	(revision 376)
+++ codes/icosagcm/trunk/src/etat0_dcmip3.f90	(revision 377)
@@ -8,5 +8,4 @@
   
   SUBROUTINE compute_etat0(ngrid,lon,lat, phis,ps,temp,ulon,ulat,geopot,q)
-    USE genmod
     USE dcmip_initial_conditions_test_1_2_3
     USE disvert_mod
@@ -27,5 +26,5 @@
     pp=peq
     DO ij=1,ngrid
-       CALL test3_gravity_wave(lon(ij),lat(ij),pp,dummy,0, &
+       CALL test3_gravity_wave(scale_factor, lon(ij),lat(ij),pp,dummy,0, &
             dummy,dummy,dummy,dummy,phis(ij),ps(ij),dummy,dummy)
     END DO
@@ -33,5 +32,5 @@
        DO ij=1,ngrid
           pp = ap(l) + bp(l)*ps(ij) ! half-layer pressure
-          CALL test3_gravity_wave(lon(ij),lat(ij),pp,zz,0, &
+          CALL test3_gravity_wave(scale_factor, lon(ij),lat(ij),pp,zz,0, &
                dummy,dummy,dummy,dummy,dummy,dummy,dummy,dummy)
           geopot(ij,l) = g*zz ! initialize geopotential for NH
@@ -41,5 +40,5 @@
        DO ij=1,ngrid
           pp = .5*(ap(l)+ap(l+1)) + .5*(bp(l)+bp(l+1))*ps(ij) ! full-layer pressure
-          CALL test3_gravity_wave(lon(ij),lat(ij),pp,dummy,0, &
+          CALL test3_gravity_wave(scale_factor, lon(ij),lat(ij),pp,dummy,0, &
                ulon(ij,l),ulat(ij,l),dummy,Temp(ij,l),dummy,dummy,dummy,dummy)
        END DO
Index: codes/icosagcm/trunk/src/etat0_temperature.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0_temperature.f90	(revision 376)
+++ codes/icosagcm/trunk/src/etat0_temperature.f90	(revision 377)
@@ -3,6 +3,11 @@
   USE icosa, ONLY: llm,nqtot
   IMPLICIT NONE
+  PRIVATE
+
   REAL(rstd), SAVE, ALLOCATABLE :: t_profile(:)
 !$OMP THREADPRIVATE(t_profile)
+
+  PUBLIC :: getin_etat0, compute_etat0
+
 CONTAINS
   
@@ -13,5 +18,4 @@
     USE transfert_omp_mod, ONLY: bcast_omp
     USE free_unit_mod, ONLY : free_unit
-    IMPLICIT NONE
     INTEGER :: unit,ok
     INTEGER :: l
@@ -59,5 +63,4 @@
   SUBROUTINE compute_etat0(ngrid, phis, ps, temp, ulon, ulat, q)
     USE earth_const, ONLY: preff
-    IMPLICIT NONE  
     INTEGER, INTENT(IN)    :: ngrid
     REAL(rstd),INTENT(OUT) :: phis(ngrid)
Index: codes/icosagcm/trunk/src/observable.f90
===================================================================
--- codes/icosagcm/trunk/src/observable.f90	(revision 376)
+++ codes/icosagcm/trunk/src/observable.f90	(revision 377)
@@ -169,4 +169,6 @@
     REAL(rstd) :: F_el(3*iim*jjm,llm+1)
     REAL(rstd) :: uu_right, uu_lup, uu_ldown, W_el, DePhil
+    ! NB : u and uh are not in DEC form, they are normal components    
+    ! => we must divide by de
     IF(hydrostatic) THEN
        uh(:,:)=u(:,:)
@@ -178,13 +180,13 @@
              W_el   = .5*( W(ij,l)+W(ij+t_right,l) )
              DePhil = ne_right*(Phi(ij+t_right,l)-Phi(ij,l))
-             F_el(ij+u_right,l)   = DePhil*W_el
+             F_el(ij+u_right,l)   = DePhil*W_el/de(ij+u_right)
              ! Compute on edge 'lup'
              W_el   = .5*( W(ij,l)+W(ij+t_lup,l) )
              DePhil = ne_lup*(Phi(ij+t_lup,l)-Phi(ij,l))
-             F_el(ij+u_lup,l)   = DePhil*W_el
+             F_el(ij+u_lup,l)   = DePhil*W_el/de(ij+u_lup)
              ! Compute on edge 'ldown'
              W_el   = .5*( W(ij,l)+W(ij+t_ldown,l) )
              DePhil = ne_ldown*(Phi(ij+t_ldown,l)-Phi(ij,l))
-             F_el(ij+u_ldown,l) = DePhil*W_el
+             F_el(ij+u_ldown,l) = DePhil*W_el/de(ij+u_ldown)
           END DO
        END DO
@@ -192,5 +194,5 @@
        DO l=ll_begin, ll_end ! compute on k levels (full levels)
           DO ij=ij_begin_ext, ij_end_ext
-             ! w = vertical momentum = g^-2*dPhi/dt = g*uz
+             ! w = vertical momentum = g^-2*dPhi/dt = uz/g
              uz(ij,l) = (.5*g)*(W(ij,l)+W(ij,l+1))/rhodz(ij,l)
              ! uh = u-w.grad(Phi) = u - uz.grad(z)
Index: codes/icosagcm/trunk/src/timeloop_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 376)
+++ codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 377)
@@ -11,5 +11,6 @@
   
   INTEGER, PARAMETER :: itau_sync=10
-  TYPE(t_message),SAVE :: req_ps0, req_mass0, req_theta_rhodz0, req_u0, req_q0
+  TYPE(t_message),SAVE :: req_ps0, req_mass0, req_theta_rhodz0, req_u0, req_q0, req_W0, req_geopot0
+  LOGICAL, SAVE :: positive_theta
 
   PUBLIC :: init_timeloop, timeloop
@@ -33,5 +34,12 @@
     IF (xios_output) itau_out=1
     IF (.NOT. enable_io) itau_out=HUGE(itau_out)
-    
+
+    positive_theta=.FALSE.
+    CALL getin('positive_theta',positive_theta)
+    IF(positive_theta .AND. nqtot<1) THEN
+       PRINT *, 'nqtot must be >0 if positive_theta is .TRUE.'
+       STOP
+    END IF
+
     def='ARK2.3'
     CALL getin('time_scheme',def)
@@ -147,4 +155,6 @@
     CALL init_message(f_u,req_e0_vect,req_u0)
     CALL init_message(f_q,req_i0,req_q0)
+    CALL init_message(f_geopot,req_i0,req_geopot0)
+    CALL init_message(f_W,req_i0,req_W0)
 
   END SUBROUTINE init_timeloop
@@ -198,4 +208,6 @@
     fluxt_zero=.TRUE.
 
+    IF(positive_theta) CALL copy_theta_to_q1(f_theta_rhodz,f_rhodz,f_q)
+
     !$OMP MASTER
     CALL SYSTEM_CLOCK(start_clock)
@@ -226,5 +238,11 @@
           CALL wait_message(req_u0)
           CALL send_message(f_q,req_q0) 
-          CALL wait_message(req_q0) 
+          CALL wait_message(req_q0)
+          IF(.NOT. hydrostatic) THEN
+!             CALL send_message(f_geopot,req_geopot0)
+!             CALL wait_message(req_geopot0)
+!             CALL send_message(f_W,req_W0)
+!             CALL wait_message(req_W0)
+          END IF
        ENDIF
 
@@ -285,4 +303,5 @@
              END DO
           ENDIF
+          IF(positive_theta) CALL copy_q1_to_theta(f_theta_rhodz,f_rhodz,f_q)
        END IF
        
@@ -341,3 +360,33 @@
   END SUBROUTINE print_iteration
 
+  SUBROUTINE copy_theta_to_q1(f_theta_rhodz,f_rhodz,f_q)
+    TYPE(t_field),POINTER :: f_theta_rhodz(:),f_rhodz(:), f_q(:)
+    REAL(rstd), POINTER :: theta_rhodz(:,:), rhodz(:,:), q(:,:,:)
+    INTEGER :: ind
+    DO ind=1, ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       theta_rhodz=f_theta_rhodz(ind)
+       rhodz=f_rhodz(ind)
+       q=f_q(ind)
+       q(:,:,1)  = theta_rhodz(:,:)/rhodz(:,:)
+    END DO
+  END SUBROUTINE copy_theta_to_q1
+
+  SUBROUTINE copy_q1_to_theta(f_theta_rhodz,f_rhodz,f_q)
+    TYPE(t_field),POINTER :: f_theta_rhodz(:),f_rhodz(:), f_q(:)
+    REAL(rstd), POINTER :: theta_rhodz(:,:), rhodz(:,:), q(:,:,:)
+    INTEGER :: ind
+    DO ind=1, ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       theta_rhodz=f_theta_rhodz(ind)
+       rhodz=f_rhodz(ind)
+       q=f_q(ind)
+       theta_rhodz(:,:) = rhodz(:,:)*q(:,:,1)
+    END DO
+  END SUBROUTINE copy_q1_to_theta
+
 END MODULE timeloop_gcm_mod
