Index: /codes/icosagcm/trunk/src/etat0.f90
===================================================================
--- /codes/icosagcm/trunk/src/etat0.f90	(revision 191)
+++ /codes/icosagcm/trunk/src/etat0.f90	(revision 192)
@@ -7,4 +7,5 @@
   SUBROUTINE etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_q)
     USE icosa
+    USE mpipara, ONLY : is_mpi_root
     USE disvert_mod
     USE etat0_williamson_mod, ONLY : etat0_williamson_new
@@ -61,4 +62,10 @@
        CALL etat0_dcmip3(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
      CASE ('dcmip4')
+        IF(nqtot<2) THEN
+           IF (is_mpi_root)  THEN
+              PRINT *, "nqtot must be at least 2 for test case DCMIP4"
+           END IF
+           STOP
+        END IF
        CALL etat0_dcmip4(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
      CASE ('dcmip5')
Index: /codes/icosagcm/trunk/src/etat0_dcmip5.f90
===================================================================
--- /codes/icosagcm/trunk/src/etat0_dcmip5.f90	(revision 191)
+++ /codes/icosagcm/trunk/src/etat0_dcmip5.f90	(revision 192)
@@ -96,5 +96,5 @@
     lonc=Pi 
 
-    Tv0=T0/(1+0.608*q0)
+    Tv0=T0*(1+0.608*q0)
     Tvt=Tv0-Gamma*zt
   
Index: /codes/icosagcm/trunk/src/physics_dcmip.f90
===================================================================
--- /codes/icosagcm/trunk/src/physics_dcmip.f90	(revision 191)
+++ /codes/icosagcm/trunk/src/physics_dcmip.f90	(revision 192)
@@ -16,5 +16,5 @@
   IMPLICIT NONE
 
-    testcase=1
+    testcase=1 ! OK for 4.2 (moist baroclinic instability)
     CALL getin("dcmip_physics",testcase)
     CALL allocate_field(f_out_i,field_t,type_real,llm)
@@ -106,5 +106,5 @@
     INTEGER :: i,j,l,ij
      
-    PRINT *,'Entering in DCMIP physics'    
+!    PRINT *,'Entering in DCMIP physics'    
     CALL compute_pression(ps,p,0)
     CALL compute_exner(ps,p,pks,pk,0)
@@ -112,7 +112,18 @@
     CALL compute_geopotential(phis,pks,pk,theta,phi,0)
     CALL compute_theta_rhodz2temperature(ps,theta_rhodz,T,0)
-   CALL compute_wind_centered(ue,uc)
-   CALL compute_wind_centered_lonlat_compound(uc, u, v)
-
+    ! Reconstruct wind vector at hexagons
+    CALL compute_wind_centered(ue,uc)
+    CALL compute_wind_centered_lonlat_compound(uc, u, v)
+
+    ! Convert from Tv to T
+    DO l=1,llm
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          T(ij,l)=T(ij,l)/(1+0.608*q(ij,l))
+        ENDDO
+      ENDDO
+    ENDDO       
+    
     DO j=jj_begin,jj_end
       DO i=ii_begin,ii_end
@@ -122,4 +133,7 @@
     ENDDO
    
+    ! bottom-up indexing (DYNAMICO) : u,utemp, v,vtemp                                                                     
+    ! top-down vertical indexing (DCMIP) : ufi, vfi, ...                                                             
+    ! => copy fields and mirror vertical index                                                       
     DO l=1,llm+1
       DO j=jj_begin,jj_end
@@ -131,4 +145,5 @@
     ENDDO
     
+    ! Pressure inside layers                                                                                             
     DO l=1,llm
       DO j=jj_begin,jj_end
@@ -140,4 +155,5 @@
     ENDDO
     
+    ! Pressure difference between two layers
     DO l=1,llm
       DO j=jj_begin,jj_end
@@ -149,17 +165,5 @@
     ENDDO       
        
-    
-!    ufi=u
-!    vfi=v    
-    
-    DO l=1,llm
-      DO j=jj_begin,jj_end
-        DO i=ii_begin,ii_end
-          ij=(j-1)*iim+i
-          T(ij,l)=T(ij,l)/(1+0.608*q(ij,l))
-        ENDDO
-      ENDDO
-    ENDDO       
-    
+    ! Copy T,u,v,q for input to physics
     DO l=1,llm
       DO j=jj_begin,jj_end
@@ -174,9 +178,7 @@
     ENDDO       
     
-!    q=0
-!    out_i=T
-    
     CALL simple_physics(iim*jjm, llm, dt, lat, tfi, qfi , ufi, vfi, pmid, pint, pdel, 1/pdel, ps, precl, testcase) 
     
+    ! Copy back T,u,v,q and mirror vertical index
     DO l=1,llm
       DO j=jj_begin,jj_end
@@ -192,4 +194,5 @@
 
 
+    ! Convert back T to Tv
     DO l=1,llm
       DO j=jj_begin,jj_end
@@ -201,9 +204,9 @@
     ENDDO       
 
-!    out_i=q
-        
+    ! Compute velocity update at hexagons                                                                                
     utemp=utemp-u
     vtemp=vtemp-v
 
+    ! lon-lat -> 3D
     DO l=1,llm
       DO j=jj_begin,jj_end
@@ -215,6 +218,5 @@
     ENDDO
 
-!    out_i=ufi
-    
+    ! Update velocity at velocity points
     DO l=1,llm
       DO j=jj_begin,jj_end
@@ -519,5 +521,4 @@
       end do
 
-      IF (test==0) return
 !===============================================================================
 ! Send variables to history file - THIS PROCESS WILL BE MODEL SPECIFIC
