Index: codes/icosagcm/trunk/src/caldyn_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 95)
+++ codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 96)
@@ -641,4 +641,5 @@
     USE theta2theta_rhodz_mod
     USE pression_mod
+    USE omega_mod
     USE write_field
     TYPE(t_field),POINTER :: f_ps(:), f_phis(:), f_u(:), f_theta_rhodz(:), f_dps(:), &
@@ -650,4 +651,7 @@
     CALL vorticity(f_u,f_buf_v)
     CALL writefield("vort",f_buf_v)
+
+    CALL w_omega(f_ps, f_u, f_buf_i)
+    CALL writefield('omega', f_buf_i)
     
     ! Temperature
Index: codes/icosagcm/trunk/src/omega.f90
===================================================================
--- codes/icosagcm/trunk/src/omega.f90	(revision 96)
+++ codes/icosagcm/trunk/src/omega.f90	(revision 96)
@@ -0,0 +1,146 @@
+MODULE omega_mod
+
+  USE icosa
+  PRIVATE
+
+  PUBLIC :: w_omega, compute_omega
+
+CONTAINS
+
+  SUBROUTINE w_omega(f_ps, f_u, f_omega) ! Compute omega = Dp/Dt
+    TYPE(t_field),POINTER :: f_ps(:), f_u(:), f_omega(:)
+    INTEGER :: ind
+    REAL(rstd),POINTER :: ps(:), u(:,:), om(:,:)
+    DO ind=1,ndomain
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       ps=f_ps(ind)
+       u=f_u(ind)
+       om=f_omega(ind)
+       CALL compute_omega(ps,u,om)
+    END DO
+  END SUBROUTINE W_omega
+
+  SUBROUTINE compute_omega(ps,u, w)
+    USE disvert_mod, ONLY : ap,bp
+    REAL(rstd),INTENT(IN) :: u(iim*3*jjm,llm), ps(iim*jjm)
+    REAL(rstd),INTENT(OUT):: w(iim*jjm,llm)
+    REAL(rstd):: convm(iim*jjm,llm+1)
+    REAL(rstd):: p(iim*jjm,llm+1), rhodz(iim*jjm,llm), Fe(iim*3*jjm,llm)
+    REAL(rstd):: ugradps
+    DO    l    = 1, llm+1
+       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
+       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 mass flux
+    DO l = 1, llm
+       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
+    
+!!! mass flux convergence computation  
+    
+    ! horizontal convergence
+    DO l = 1, llm
+       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
+       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
+    convm(:,llm+1)=0.
+
+!!! Compute dps
+    !  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 
+    !      convm(ij,llm+1)=0.
+    !    ENDDO
+    !  ENDDO
+    
+    ! Compute Omega = Dp/Dt
+    ! with p = A(eta)+B(eta)ps
+    !      Dp/Dt = dp/deta.Deta/Dt + B(eta)Dps/Dt
+    !            = -mg.Deta/Dt + B.Dps/Dt
+    ! By definition the mass flux through model levels is W=m.Deta/Dt with m=-1/g dp/deta
+    ! therefore 
+    !  Dp/Dt = -g.W + B.dps/dt + Bu.grad ps
+    !        = B.u.grad ps - g*convm
+    
+!!! Compute vertical flux through model layers
+    !  DO l = 1,llm-1
+    !    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 ) ! g.W = g.convm + B dps/dt
+    !      ENDDO
+    !    ENDDO
+    !  ENDDO
+    
+    
+!!! Compute omega
+    ! -grad ps : ( ne(ij,ldown)*ps(ij,l) + ne(ij+t_ldown,rup)*ps(ij+t_ldown,l) ) ) / de(ij+u_ldown)
+
+    DO l = 1,llm
+       DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+             toto = 1
+             ij=(j-1)*iim+i
+             ugradps = &
+                  le(ij+u_right)*u(ij+u_right,l)*( ne(ij,right)*ps(ij) + ne(ij+t_right,left)*ps(ij+t_right) )  &
+                  + le(ij+u_rup)*u(ij+u_rup,l)*( ne(ij,rup)*ps(ij) + ne(ij+t_rup,ldown)*ps(ij+t_rup) )  &
+                  + le(ij+u_lup)*u(ij+u_lup,l)*( ne(ij,lup)*ps(ij) + ne(ij+t_lup,rdown)*ps(ij+t_lup) )         &
+                  + le(ij+u_left)*u(ij+u_left,l)*( ne(ij,left)*ps(ij) + ne(ij+t_left,right)*ps(ij+t_left) )    &
+                  + le(ij+u_ldown)*u(ij+u_ldown,l)*( ne(ij,ldown)*ps(ij) + ne(ij+t_ldown,rup)*ps(ij+t_ldown) ) &
+                  + le(ij+u_rdown)*u(ij+u_rdown,l)*( ne(ij,rdown)*ps(ij) + ne(ij+t_rdown,lup)*ps(ij+t_rdown) )
+             ugradps = .5*(bp(l)+bp(l+1)) *ugradps/(-4*Ai(ij))  ! sign convention as in Ringler et al. 2010, Eq. 22 p.3072
+             w( ij, l) =  ugradps - .5*(convm( ij,l+1)+convm(ij,l)) 
+          ENDDO
+       ENDDO
+    ENDDO
+    
+  END SUBROUTINE compute_omega
+  
+END MODULE omega_mod
