Index: codes/icosagcm/trunk/arch/arch-X64_YELLOWSTONE.fcm
===================================================================
--- codes/icosagcm/trunk/arch/arch-X64_YELLOWSTONE.fcm	(revision 404)
+++ codes/icosagcm/trunk/arch/arch-X64_YELLOWSTONE.fcm	(revision 405)
@@ -7,5 +7,5 @@
 %FPP_DEF             KEY_NONE CPP_NETCDF4 CPP_NF90_DEF_VAR_CHUNKING
 %BASE_FFLAGS         -i4 -r8 -auto -align all -I${MKLROOT}/include
-%PROD_FFLAGS         -g -traceback -O3 -vec-report2
+%PROD_FFLAGS         -g -traceback -O3
 %DEV_FFLAGS          -g -O1 -traceback
 %DEBUG_FFLAGS        -g -traceback -check bounds -fp-model strict
Index: codes/icosagcm/trunk/src/caldyn_kernels_base.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_kernels_base.f90	(revision 404)
+++ codes/icosagcm/trunk/src/caldyn_kernels_base.f90	(revision 405)
@@ -37,5 +37,5 @@
 
     INTEGER :: i,j,ij,l
-    REAL(rstd) :: Rd, p_ik, exner_ik, temp_ik
+    REAL(rstd) :: Rd, p_ik, exner_ik, temp_ik, qv, chi, Rmix
     INTEGER    :: ij_omp_begin_ext, ij_omp_end_ext
 
@@ -79,22 +79,45 @@
        ! uppermost layer
        
-       !DIR$ SIMD
-       DO ij=ij_omp_begin_ext,ij_omp_end_ext
-          pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm)
-       END DO
-       ! other layers
-       DO l = llm-1, 1, -1
-          !DIR$ SIMD
-          DO ij=ij_omp_begin_ext,ij_omp_end_ext         
-             pk(ij,l) = pk(ij,l+1) + (.5*g)*(rhodz(ij,l)+rhodz(ij,l+1))
-          END DO
-       END DO
-       ! surface pressure (for diagnostics)
-       IF(caldyn_eta==eta_lag) THEN
-          DO ij=ij_omp_begin_ext,ij_omp_end_ext         
-             ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1)
-          END DO
-       END IF
-       
+       SELECT CASE(caldyn_thermo)
+          CASE(thermo_theta, thermo_entropy)
+             !DIR$ SIMD
+             DO ij=ij_omp_begin_ext,ij_omp_end_ext
+                pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm)
+             END DO
+             ! other layers
+             DO l = llm-1, 1, -1
+                !DIR$ SIMD
+                DO ij=ij_omp_begin_ext,ij_omp_end_ext         
+                   pk(ij,l) = pk(ij,l+1) + (.5*g)*(rhodz(ij,l)+rhodz(ij,l+1))
+                END DO
+             END DO
+             ! surface pressure (for diagnostics)
+             IF(caldyn_eta==eta_lag) THEN
+                DO ij=ij_omp_begin_ext,ij_omp_end_ext         
+                   ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1)
+                END DO
+             END IF
+          CASE(thermo_moist) ! theta(ij,l,2) = qv = mv/md
+             !DIR$ SIMD
+             DO ij=ij_omp_begin_ext,ij_omp_end_ext
+                pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm)*(1.+theta(ij,l,2))
+             END DO
+             ! other layers
+             DO l = llm-1, 1, -1
+                !DIR$ SIMD
+                DO ij=ij_omp_begin_ext,ij_omp_end_ext         
+                   pk(ij,l) = pk(ij,l+1) + (.5*g)*(          &
+                        rhodz(ij,l)  *(1.+theta(ij,l,2)) +   &
+                        rhodz(ij,l+1)*(1.+theta(ij,l+1,2)) )
+                END DO
+             END DO
+             ! surface pressure (for diagnostics)
+             IF(caldyn_eta==eta_lag) THEN
+                DO ij=ij_omp_begin_ext,ij_omp_end_ext         
+                   ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1)*(1.+theta(ij,l,2))
+                END DO
+             END IF
+          END SELECT
+
        DO l = 1,llm
           SELECT CASE(caldyn_thermo)
@@ -117,4 +140,16 @@
                 geopot(ij,l+1) = geopot(ij,l) + (g*Rd)*rhodz(ij,l)*temp_ik/p_ik
              ENDDO
+          CASE(thermo_moist) ! theta is moist pseudo-entropy per dry air mass
+             DO ij=ij_omp_begin_ext,ij_omp_end_ext
+                p_ik = pk(ij,l)
+                qv = theta(ij,l,2) ! water vaper mixing ratio = mv/md
+                Rmix = Rd+qv*Rv
+                chi = ( theta(ij,l,1) + Rmix*log(p_ik/preff) ) / (cpp + qv*cppv) ! log(T/Treff)
+                temp_ik = Treff*exp(chi)
+                pk(ij,l) = temp_ik
+                ! specific volume v = R*T/p = dphi/g/rhodz
+                ! R = (Rd + qv.Rv)/(1+qv)
+                geopot(ij,l+1) = geopot(ij,l) + g*Rmix*rhodz(ij,l)*temp_ik/(p_ik*(1+qv))
+             ENDDO
           CASE DEFAULT
              STOP
@@ -132,5 +167,5 @@
   SUBROUTINE compute_caldyn_vert(u,theta,rhodz,convm, wflux,wwuu, dps,dtheta_rhodz,du)
     REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
-    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)
+    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm,nqdyn)
     REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
     REAL(rstd),INTENT(INOUT)  :: convm(iim*jjm,llm)  ! mass flux convergence
@@ -139,9 +174,9 @@
     REAL(rstd),INTENT(INOUT) :: wwuu(iim*3*jjm,llm+1)
     REAL(rstd),INTENT(INOUT) :: du(iim*3*jjm,llm)
-    REAL(rstd),INTENT(INOUT) :: dtheta_rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(INOUT) :: dtheta_rhodz(iim*jjm,llm,nqdyn)
     REAL(rstd),INTENT(OUT) :: dps(iim*jjm)
 
     ! temporary variable    
-    INTEGER :: i,j,ij,l
+    INTEGER :: i,j,ij,l,iq
     REAL(rstd) :: p_ik, exner_ik
     INTEGER    :: ij_omp_begin, ij_omp_end
@@ -204,18 +239,20 @@
     !$OMP BARRIER
 
-    DO l=ll_begin,ll_endm1
-       !DIR$ SIMD
-       DO ij=ij_begin,ij_end         
-          dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  -   0.5 * (  wflux(ij,l+1) * (theta(ij,l) + theta(ij,l+1)))
-       ENDDO
-    ENDDO
-
-    DO l=ll_beginp1,ll_end
-       !DIR$ SIMD
-       DO ij=ij_begin,ij_end         
-          dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  +   0.5 * ( wflux(ij,l  ) * (theta(ij,l-1) + theta(ij,l) ) )
-       ENDDO
-    ENDDO
-
+    DO iq=1,nqdyn
+       DO l=ll_begin,ll_endm1
+       !DIR$ SIMD
+          DO ij=ij_begin,ij_end         
+             dtheta_rhodz(ij, l, iq) = dtheta_rhodz(ij, l, iq)  -   0.5 * &
+                  ( wflux(ij,l+1) * (theta(ij,l,iq) + theta(ij,l+1,iq)))
+          END DO
+       END DO
+       DO l=ll_beginp1,ll_end
+          !DIR$ SIMD
+          DO ij=ij_begin,ij_end         
+             dtheta_rhodz(ij, l, iq) = dtheta_rhodz(ij, l, iq)  +   0.5 * &
+                  ( wflux(ij,l) * (theta(ij,l-1,iq) + theta(ij,l,iq) ) )
+          END DO
+       END DO
+    END DO
 
     ! Compute vertical transport
Index: codes/icosagcm/trunk/src/caldyn_kernels_hevi.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_kernels_hevi.f90	(revision 404)
+++ codes/icosagcm/trunk/src/caldyn_kernels_hevi.f90	(revision 405)
@@ -382,18 +382,20 @@
     REAL(rstd),INTENT(INOUT) :: u(iim*3*jjm,llm)   ! OUT if tau>0
     REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm)
-    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm)
+    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm,nqdyn)
     REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm)
     REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1)
     REAL(rstd),INTENT(INOUT)   :: du(iim*3*jjm,llm)
     REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
+    REAL(rstd) :: berniv(iim*jjm,llm)  ! moist Bernoulli function
 
     INTEGER :: i,j,ij,l
-    REAL(rstd) :: due_right, due_lup, due_ldown
+    REAL(rstd) :: Rd, qv, temp, chi, nu, due_right, due_lup, due_ldown
 
     CALL trace_start("compute_caldyn_fast")
+
+    Rd=cpp*kappa
 
     ! Compute Bernoulli term
     IF(boussinesq) THEN
-
        DO l=ll_begin,ll_end
           !DIR$ SIMD
@@ -417,6 +419,22 @@
              DO ij=ij_begin,ij_end         
                 berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
-                     + pk(ij,l)*(cpp-theta(ij,l)) ! pk=temperature, theta=entropy
+                     + pk(ij,l)*(cpp-theta(ij,l,1)) ! pk=temperature, theta=entropy
              END DO
+          CASE(thermo_moist) 
+             !DIR$ SIMD
+             DO ij=ij_begin,ij_end         
+                ! du/dt = grad(Bd)+rv.grad(Bv)+s.grad(T)
+                ! Bd = Phi + gibbs_d
+                ! Bv = Phi + gibbs_v
+                ! pk=temperature, theta=entropy
+                qv = theta(ij,l,2)
+                temp = pk(ij,l)
+                chi = log(temp/Treff)
+                nu = (chi*(cpp+qv*cppv)-theta(ij,l,1))/(Rd+qv*Rv) ! log(p/preff)
+                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
+                     + temp*(cpp*(1.-chi)+Rd*nu)
+                berniv(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
+                     + temp*(cppv*(1.-chi)+Rv*nu)
+            END DO
           END SELECT
        END DO
@@ -426,24 +444,54 @@
 !!! u:=u+tau*du, du = -grad(B)-theta.grad(pi)
     DO l=ll_begin,ll_end
-       !DIR$ SIMD
-       DO ij=ij_begin,ij_end         
-          due_right =  0.5*(theta(ij,l)+theta(ij+t_right,l))  &
-                          *(pk(ij+t_right,l)-pk(ij,l))        &
-                         +  berni(ij+t_right,l)-berni(ij,l)
-          due_lup = 0.5*(theta(ij,l)+theta(ij+t_lup,l))    &
-                       *(pk(ij+t_lup,l)-pk(ij,l))          &
-                         +  berni(ij+t_lup,l)-berni(ij,l)
-          due_ldown = 0.5*(theta(ij,l)+theta(ij+t_ldown,l)) &
-                         *(pk(ij+t_ldown,l)-pk(ij,l))      &
-                        +   berni(ij+t_ldown,l)-berni(ij,l)
-          du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right
-          du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup
-          du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown
-          u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l)
-          u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l)
-          u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l)
-       ENDDO
-    ENDDO
-
+       IF(caldyn_thermo == thermo_moist) THEN
+          !DIR$ SIMD
+          DO ij=ij_begin,ij_end         
+             due_right =  berni(ij+t_right,l)-berni(ij,l)       &
+                  + 0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))   &
+                       *(pk(ij+t_right,l)-pk(ij,l))             &
+                  + 0.5*(theta(ij,l,2)+theta(ij+t_right,l,2))   &
+                       *(berniv(ij+t_right,l)-berniv(ij,l))
+             
+             due_lup = berni(ij+t_lup,l)-berni(ij,l)            &
+                  + 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))     &
+                       *(pk(ij+t_lup,l)-pk(ij,l))               &
+                  + 0.5*(theta(ij,l,2)+theta(ij+t_lup,l,2))     &
+                       *(berniv(ij+t_lup,l)-berniv(ij,l))
+             
+             due_ldown = berni(ij+t_ldown,l)-berni(ij,l)        &
+                  + 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1))   &
+                       *(pk(ij+t_ldown,l)-pk(ij,l))             &
+                  + 0.5*(theta(ij,l,2)+theta(ij+t_ldown,l,2))   &
+                       *(berniv(ij+t_ldown,l)-berniv(ij,l))
+             
+             du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right
+             du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup
+             du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown
+             u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l)
+             u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l)
+             u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l)
+          END DO
+       ELSE
+          !DIR$ SIMD
+          DO ij=ij_begin,ij_end         
+             due_right =  0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))  &
+                  *(pk(ij+t_right,l)-pk(ij,l))        &
+                  +  berni(ij+t_right,l)-berni(ij,l)
+             due_lup = 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))    &
+                  *(pk(ij+t_lup,l)-pk(ij,l))          &
+                  +  berni(ij+t_lup,l)-berni(ij,l)
+             due_ldown = 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1)) &
+                  *(pk(ij+t_ldown,l)-pk(ij,l))      &
+                  +   berni(ij+t_ldown,l)-berni(ij,l)
+             du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right
+             du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup
+             du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown
+             u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l)
+             u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l)
+             u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l)
+          END DO
+       END IF
+    END DO
+       
     CALL trace_end("compute_caldyn_fast")
 
