Changeset 134 for codes/icosagcm/trunk/src/caldyn_gcm.f90

Timestamp:

02/14/13 18:26:54 (11 years ago)

Author:

dubos

Message:

Time-accumulation of mass fluxes active

File:

• codes/icosagcm/trunk/src/caldyn_gcm.f90 (modified) (20 diffs)

codes/icosagcm/trunk/src/caldyn_gcm.f90

@@ -35 +35 @@
        STOP
     END SELECT
-     IF (is_mpi_root) PRINT *, 'caldyn_conserv=',def
+    IF (is_mpi_root) PRINT *, 'caldyn_conserv=',def
 
     def='direct'
@@ -84 +84 @@
   END SUBROUTINE allocate_caldyn
    
-  SUBROUTINE caldyn(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du)
+  SUBROUTINE caldyn(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, &
+       f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du)
     USE icosa
     USE vorticity_mod
@@ -97 +98 @@
     TYPE(t_field),POINTER :: f_u(:)
     TYPE(t_field),POINTER :: f_q(:)
+    TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:)
     TYPE(t_field),POINTER :: f_dps(:)
     TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
@@ -103 +105 @@
     REAL(rstd),POINTER :: phis(:), ps(:), dps(:)
     REAL(rstd),POINTER :: theta_rhodz(:,:), dtheta_rhodz(:,:)
-    REAL(rstd),POINTER :: u(:,:), du(:,:)
+    REAL(rstd),POINTER :: u(:,:), du(:,:), hflux(:,:), wflux(:,:)
     REAL(rstd),POINTER :: p(:,:), rhodz(:,:), qu(:,:)
     INTEGER :: ind,ij
@@ -133 +135 @@
           CALL swap_geometry(ind)
           phis=f_phis(ind)
+          hflux=f_hflux(ind)
+          wflux=f_wflux(ind)
           ps=f_ps(ind)
           dps=f_dps(ind)
@@ -144 +148 @@
           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, &
+               hflux, wflux, dps, dtheta_rhodz, du)
           !$OMP END PARALLEL
        ENDDO       
@@ -155 +160 @@
           ps=f_ps(ind)
           dps=f_dps(ind)
+          hflux=f_hflux(ind)
+          wflux=f_wflux(ind)
           theta_rhodz=f_theta_rhodz(ind)
           dtheta_rhodz=f_dtheta_rhodz(ind)
@@ -165 +172 @@
           !$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)
+          CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, &
+               hflux, wflux, dps, dtheta_rhodz, du)
           !$OMP END PARALLEL
        ENDDO
@@ -276 +284 @@
   END SUBROUTINE compute_pvort
   
-  SUBROUTINE compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, dps, dtheta_rhodz, du)
+  SUBROUTINE compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, hflux, wflux, dps, dtheta_rhodz, du)
   USE icosa
   USE disvert_mod
@@ -289 +297 @@
     REAL(rstd),INTENT(IN)  :: qu(iim*3*jjm,llm)
 
-    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm), hflux(iim*3*jjm,llm)
     REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm)
     REAL(rstd),INTENT(OUT) :: dps(iim*jjm)
+    REAL(rstd),INTENT(OUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux
     
     INTEGER :: i,j,ij,l
@@ -301 +310 @@
     REAL(rstd),ALLOCATABLE,SAVE :: pk(:,:), pks(:) ! Exner function
     REAL(rstd),ALLOCATABLE,SAVE :: alpha(:,:), beta(:,:)
-    REAL(rstd),ALLOCATABLE,SAVE :: phi(:,:)              ! geopotential
-    REAL(rstd),ALLOCATABLE,SAVE :: Fe(:,:), Ftheta(:,:)  ! mass flux, theta flux
+    REAL(rstd),ALLOCATABLE,SAVE :: phi(:,:)    ! geopotential
+    REAL(rstd),ALLOCATABLE,SAVE :: Ftheta(:,:) ! theta flux
     REAL(rstd),ALLOCATABLE,SAVE :: convm(:,:)  ! mass flux convergence
-    REAL(rstd),ALLOCATABLE,SAVE :: w(:,:)      ! vertical mass flux
     REAL(rstd),ALLOCATABLE,SAVE :: berni(:,:)  ! Bernouilli function
     
@@ -319 +327 @@
     ALLOCATE(beta(iim*jjm,llm))
     ALLOCATE(phi(iim*jjm,llm))      ! geopotential
-    ALLOCATE(Fe(3*iim*jjm,llm))     ! mass flux   
     ALLOCATE(Ftheta(3*iim*jjm,llm)) ! theta flux   
     ALLOCATE(convm(iim*jjm,llm))    ! mass flux convergence
-    ALLOCATE(w(iim*jjm,llm))        ! vertical velocity      
     ALLOCATE(berni(iim*jjm,llm))    ! bernouilli term
 
@@ -342 +348 @@
       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)
-
-        Ftheta(ij+u_right,l)=0.5*(theta(ij,l)+theta(ij+t_right,l))*Fe(ij+u_right,l)
-        Ftheta(ij+u_lup,l)=0.5*(theta(ij,l)+theta(ij+t_lup,l))*Fe(ij+u_lup,l)
-        Ftheta(ij+u_ldown,l)=0.5*(theta(ij,l)+theta(ij+t_ldown,l))*Fe(ij+u_ldown,l)
+        hflux(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right)
+        hflux(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup)
+        hflux(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown)
+
+        Ftheta(ij+u_right,l)=0.5*(theta(ij,l)+theta(ij+t_right,l))*hflux(ij+u_right,l)
+        Ftheta(ij+u_lup,l)=0.5*(theta(ij,l)+theta(ij+t_lup,l))*hflux(ij+u_lup,l)
+        Ftheta(ij+u_ldown,l)=0.5*(theta(ij,l)+theta(ij+t_ldown,l))*hflux(ij+u_ldown,l)
       ENDDO
     ENDDO
@@ -358 +364 @@
         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))    
+        convm(ij,l)= 1./Ai(ij)*(ne(ij,right)*hflux(ij+u_right,l)   +  &
+                                ne(ij,rup)*hflux(ij+u_rup,l)       +  &  
+                                ne(ij,lup)*hflux(ij+u_lup,l)       +  &  
+                                ne(ij,left)*hflux(ij+u_left,l)     +  &  
+                                ne(ij,ldown)*hflux(ij+u_ldown,l)   +  &  
+                                ne(ij,rdown)*hflux(ij+u_rdown,l))    
 
         ! signe ? attention d (rho theta dz)
@@ -396 +402 @@
         ! 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 )
+        wflux( ij, l+1 ) = convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 )
       ENDDO
     ENDDO
@@ -406 +412 @@
     DO i=ii_begin,ii_end
       ij=(j-1)*iim+i
-      w(ij,1)  = 0.
+      wflux(ij,1)  = 0.
       ! dps/dt = -int(div flux)dz
       dps(ij)=-convm(ij,1) * g 
@@ -423 +429 @@
                 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))
+                uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))+                             &
+                     wee(ij+u_right,2,1)*hflux(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))+                             &
+                     wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+                           &
+                     wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+u_ldown,l))+                         &
+                     wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+u_rdown,l))+                         & 
+                     wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_ldown,l))+         &
+                     wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rdown,l))+         &
+                     wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_right,l))+         &
+                     wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rup,l))+             &
+                     wee(ij+u_right,5,2)*hflux(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))
+                uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) +                        &
+                     wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) +                      &
+                     wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) +                      &
+                     wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*(qu(ij+u_lup,l)+qu(ij+u_right,l)) +                      &
+                     wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+u_rup,l)) +                          & 
+                     wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_right,l)) +          &
+                     wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_rup,l)) +              &
+                     wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_lup,l)) +              &
+                     wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_left,l)) +            &
+                     wee(ij+u_lup,5,2)*hflux(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))
+                uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) +                      &
+                     wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+u_right,l)) +                      &
+                     wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) +                          &
+                     wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) +                          &
+                     wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+u_left,l)) +                        & 
+                     wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_lup,l)) +          &
+                     wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_left,l)) +        &
+                     wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_ldown,l)) +      &
+                     wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_rdown,l)) +      &
+                     wee(ij+u_ldown,5,2)*hflux(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)
                 
@@ -472 +478 @@
                 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)
+                uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)+                           &
+                     wee(ij+u_right,2,1)*hflux(ij+u_lup,l)+                           &
+                     wee(ij+u_right,3,1)*hflux(ij+u_left,l)+                          &
+                     wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)+                         &
+                     wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)+                         & 
+                     wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)+                 &
+                     wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)+                 &
+                     wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)+                 &
+                     wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)+                   &
+                     wee(ij+u_right,5,2)*hflux(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)
+                uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)+                        &
+                     wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)+                       &
+                     wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)+                       &
+                     wee(ij+u_lup,4,1)*hflux(ij+u_right,l)+                       &
+                     wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)+                         & 
+                     wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)+                 &
+                     wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)+                   &
+                     wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)+                   &
+                     wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)+                  &
+                     wee(ij+u_lup,5,2)*hflux(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)
+                uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)+                      &
+                     wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)+                      &
+                     wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)+                        &
+                     wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)+                        &
+                     wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)+                       & 
+                     wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)+                &
+                     wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)+               &
+                     wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)+              &
+                     wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)+              &
+                     wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)
                 du(ij+u_ldown,l) = qu(ij+u_ldown,l)*uu/de(ij+u_ldown)
      
@@ -606 +612 @@
          ! ww>0 <=> upward transport
 
-        ww            =   0.5 * w(ij,l+1) * (theta(ij,l) + theta(ij,l+1) )
+        ww            =   0.5 * wflux(ij,l+1) * (theta(ij,l) + theta(ij,l+1) )
         dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  -  ww
         dtheta_rhodz(ij,l+1) = dtheta_rhodz(ij,l+1)  +  ww
 
-        ww  = 0.5 * ( w(ij,l+1) + w(ij+t_right,l+1))
+        ww  = 0.5 * ( wflux(ij,l+1) + wflux(ij+t_right,l+1))
         uu  = u(ij+u_right,l+1) - u(ij+u_right,l)  
         du(ij+u_right, l )   = du(ij+u_right,l)   - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_right,l)))
         du(ij+u_right, l+1 ) = du(ij+u_right,l+1) - 0.5 * ww * uu / (0.5*(rhodz(ij,l+1)+rhodz(ij+t_right,l+1)))
 
-        ww  = 0.5 * ( w(ij,l+1) + w(ij+t_lup,l+1))
+        ww  = 0.5 * ( wflux(ij,l+1) + wflux(ij+t_lup,l+1))
         uu  = u(ij+u_lup,l+1) - u(ij+u_lup,l)  
         du(ij+u_lup, l )   = du(ij+u_lup,l)   - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l)))
         du(ij+u_lup, l+1 ) = du(ij+u_lup,l+1) - 0.5 * ww * uu / (0.5*(rhodz(ij,l+1)+rhodz(ij+t_lup,l+1)))
 
-        ww  = 0.5 * ( w(ij,l+1) + w(ij+t_ldown,l+1))
+        ww  = 0.5 * ( wflux(ij,l+1) + wflux(ij+t_ldown,l+1))
         uu  = u(ij+u_ldown,l+1) - u(ij+u_ldown,l)  
         du(ij+u_ldown, l )   = du(ij+u_ldown,l)   - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l)))
@@ -631 +637 @@
 !!$OMP BARRIER
 !!$OMP MASTER  
-    DEALLOCATE(theta)  ! potential temperature
-    DEALLOCATE(pk)   ! Exner function
+    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(phi)     ! geopotential
+    DEALLOCATE(Ftheta)  ! theta flux   
+    DEALLOCATE(convm)   ! mass flux convergence
     DEALLOCATE(berni)   ! bernouilli term
 !!$OMP END MASTER