Index: src/transform_omega_etadot.f90
===================================================================
--- src/transform_omega_etadot.f90	(revision 5f2e8f6871e43df4ab42c948b56fffed9b5bcd15)
+++ src/transform_omega_etadot.f90	(revision 5f2e8f6871e43df4ab42c948b56fffed9b5bcd15)
@@ -0,0 +1,116 @@
+!**********************************************************************
+! Copyright 2016                                                      *
+! Andreas Stohl, Massimo Cassiani, Petra Seibert, A. Frank,           *
+! Gerhard Wotawa,  Caroline Forster, Sabine Eckhardt, John Burkhart,  *
+! Harald Sodemann, Ignacio Pisso                                      *
+!                                                                     *
+! This file is part of FLEXPART-NorESM                                *
+!                                                                     *
+! FLEXPART-NorESM is free software: you can redistribute it           *
+! and/or modify                                                       *
+! it under the terms of the GNU General Public License as published by*
+! the Free Software Foundation, either version 3 of the License, or   *
+! (at your option) any later version.                                 *
+!                                                                     *
+! FLEXPART-NorESM is distributed in the hope that it will be useful,  *
+! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
+! GNU General Public License for more details.                        *
+!                                                                     *
+! You should have received a copy of the GNU General Public License   *
+! along with FLEXPART-NorESM.                                         *
+!  If not, see <http://www.gnu.org/licenses/>.                        * 
+!**********************************************************************
+     
+
+    subroutine transform_omega_etadot(uuh,vvh,wwh,n)
+                                      !i   i  i/o i     
+
+!* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
+!*                                                                             *
+!*  calculate eta_dot*dp_eta = wwh-dp_dt-u*dp_dx-v*dp_dy on costant eta levels *
+!*  put eta_dot*dp_eta in wwh                                                  *
+!*  here we use ps to obatin p,                                                *
+!*  other variables used and in common block are:                              *
+!*  ps(ix,jy,1,n),akz(kz),bkz(kz),ps_tplus1_and_min1(ix,jy,indj-1 or indj+1)   *
+!*  if methodw is not  1 routne is not used                                    * 
+!*                                                                             *
+!*                                                                             *
+!*     Author:                                                                 *
+!*     M. Cassiani  2016                                                       *
+!*                                                                             *
+!*                                                                             *
+!*                                                                             *
+!* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
+    
+    use par_mod
+    use com_mod  
+    
+    implicit none
+    
+    integer :: ix,jy,kz,n,ixp,ixm,deltat
+    real(kind=4) :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
+    real(kind=4) :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
+    real(kind=4) :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
+    real(kind=4) :: pressure(0:nxmax-1,0:nymax-1,nwzmax)
+    real(kind=4) :: dpressuredt(0:nxmax-1,0:nymax-1,nwzmax)
+    real(kind=4) :: etadotdpdeta(0:nxmax-1,0:nymax-1,nwzmax) !not necessary but introduced for clarity eventually we coul rewrite wwh
+    real(kind=4) :: oneoverdxinm(0:nymax-1)
+    real(kind=4) :: oneoverdyinm,dpdx,dpdy
+    
+    
+    
+     do jy=1,ny-2
+       oneoverdxinm(jy)=dxconst/cos((real(jy)*dy+ylat0)*pi180)
+     end do
+        oneoverdyinm=dyconst 
+         
+     
+ 
+      do jy=0,nymin1
+        do ix=0, nxfield-1           
+          do kz=2,nuvz  !here nuvz is 27 (top of the domain) we skyp level one (extra kayer at the ground)
+            pressure(ix,jy,kz)=akz(kz)+bkz(kz)*ps(ix,jy,1,n) 
+            dpressuredt(ix,jy,kz)=bkz(kz)*(ps_tplus1_and_min1(ix,jy,2)-ps_tplus1_and_min1(ix,jy,1)) !ps_tplus1_and_min1 is already divided for the time interval
+          end do
+        end do
+      end do
+      
+      
+      do jy=1,nymin1-1 !skip the poles
+        do ix=0, nxfield-1  
+          ixp=ix+1
+          ixm=ix-1
+          if (ixm.eq.-1)  ixm=nxfield-1 !special treatment for 0 and nxfield-1,
+          if (ixp.eq.nxfield) ixp=0    !special treatment for 0 and nxfield-1,
+          do kz=2,nuvz  !here nuvz is 27 (top of the domain) we skip level one (extra layer at the ground)
+            dpdx=(pressure(ixp,jy,kz)-pressure(ixm,jy,kz))*oneoverdxinm(jy)*0.5 !ix increase eastward  i.e. in the direction positive u (positive towards the east)
+            dpdy=(pressure(ix,jy+1,kz)-pressure(ix,jy-1,kz))*oneoverdyinm*0.5  !jy incrrease northward i.e. in the direction positive v (positivwe towards the north)
+            etadotdpdeta(ix,jy,kz)=wwh(ix,jy,kz)-dpressuredt(ix,jy,kz)-uuh(ix,jy,kz)*dpdx-vvh(ix,jy,kz)*dpdy
+            ixp=ixp
+          end do
+        end do
+      end do
+      
+      do jy=1,nymin1-1 !skip the poles
+        do ix=0, nxfield-1  !special treatment for 0 and nxfield-1,          
+          do kz=2,nuvz  !here nuvz is 27 (top of the domain) we skip level one (extra layer at the ground)
+            wwh(ix,jy,kz)=etadotdpdeta(ix,jy,kz)          
+          end do
+        end do
+      end do
+       
+      wwh(:,0,:)=0.
+      wwh(:,nymin1,:)=0.
+       
+      do ix=0, nxfield-1  !         
+        do kz=2,nuvz  !here nuvz is 27 (top of the domain) we skip level one (extra layer at the ground)
+          wwh(ix,0,kz)=wwh(ix,0,kz)+etadotdpdeta(ix,1,kz)/nxfield  !set to averaged value of next circle 
+          wwh(ix,nymin1,kz)= wwh(ix,nymin1,kz)+etadotdpdeta(ix,nymin1-1,kz)/nxfield  !set to averaged value of previous circle 
+        end do
+      end do
+      
+       
+       
+      return
+      end
