Index: codes/icosagcm/trunk/src/etat0.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0.f90	(revision 381)
+++ codes/icosagcm/trunk/src/etat0.f90	(revision 382)
@@ -24,4 +24,5 @@
     USE etat0_williamson_mod, ONLY : getin_etat0_williamson=>getin_etat0
     USE etat0_temperature_mod, ONLY: getin_etat0_temperature=>getin_etat0
+    USE etat0_dcmip2016_baroclinic_wave_mod, ONLY : getin_etat0_dcmip2016_baroclinic_wave=>getin_etat0
     ! Ad hoc interfaces
     USE etat0_academic_mod, ONLY : etat0_academic=>etat0  
@@ -77,4 +78,6 @@
        autoinit_mass=.FALSE.
        CALL getin_etat0_williamson
+    CASE ('dcmip2016_baroclinic_wave')
+        CALL getin_etat0_dcmip2016_baroclinic_wave
     CASE DEFAULT
        collocated=.FALSE.
@@ -183,4 +186,5 @@
     USE etat0_williamson_mod, ONLY : compute_w91_6 => compute_etat0
     USE etat0_temperature_mod, ONLY: compute_etat0_temperature => compute_etat0
+    USE etat0_dcmip2016_baroclinic_wave_mod, ONLY : compute_dcmip2016_baroclinic_wave => compute_etat0
     IMPLICIT NONE
     REAL(rstd),INTENT(INOUT) :: ps(iim*jjm)
@@ -252,4 +256,7 @@
        CALL compute_w91_6(3*iim*jjm,lon_e,lat_e, phis_e, mass_e(:,1), temp_e(:,1), ulon_e(:,1), ulat_e(:,1))
        autoinit_mass = .FALSE. ! do not overwrite mass
+    CASE('dcmip2016_baroclinic_wave')
+       CALL compute_dcmip2016_baroclinic_wave(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q)
+       CALL compute_dcmip2016_baroclinic_wave(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
     END SELECT
 
Index: codes/icosagcm/trunk/src/etat0_dcmip2016_baroclinic_wave.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0_dcmip2016_baroclinic_wave.f90	(revision 382)
+++ codes/icosagcm/trunk/src/etat0_dcmip2016_baroclinic_wave.f90	(revision 382)
@@ -0,0 +1,107 @@
+MODULE etat0_dcmip2016_baroclinic_wave_mod
+  USE icosa
+  IMPLICIT NONE
+  PRIVATE
+  
+  REAL(rstd),PARAMETER :: eta0=0.252
+  REAL(rstd),PARAMETER :: etat=0.2
+  REAL(rstd),PARAMETER :: ps0=1e5
+  REAL(rstd),PARAMETER :: u0=35
+  REAL(rstd),PARAMETER :: T0=288
+  REAL(rstd),PARAMETER :: DeltaT=4.8e5
+  REAL(rstd),PARAMETER :: Rd=287
+  REAL(rstd),PARAMETER :: Gamma=0.005
+  REAL(rstd),PARAMETER :: up0=1
+  REAL(rstd),PARAMETER :: lonc=Pi/9, latc=2*Pi/9, latw=2*Pi/9
+  REAL(rstd),PARAMETER :: pw=34000
+  REAL(rstd),PARAMETER :: q0=0.021
+  
+  INTEGER,SAVE :: testcase
+  !$OMP THREADPRIVATE(testcase)  
+  
+  PUBLIC getin_etat0, compute_etat0
+
+CONTAINS
+
+  SUBROUTINE getin_etat0
+    USE mpipara, ONLY : is_mpi_root
+    USE tracer_mod
+    IF(nqtot<5) THEN
+       IF (is_mpi_root)  THEN
+          PRINT *, "nqtot must be at least 5 for test case dcmip2016_baroclinic_wave"
+       END IF
+       STOP
+    END IF
+    CALL set_advection_scheme(3,advect_none)
+    
+  END SUBROUTINE getin_etat0
+
+  SUBROUTINE compute_etat0(ngrid,lon,lat, phis,ps,temp,ulon,ulat,q)
+    USE icosa
+    USE disvert_mod
+    USE omp_para
+    USE terminator
+    USE tracer_mod
+    INTEGER, INTENT(IN) :: ngrid
+    REAL(rstd),INTENT(IN) :: lon(ngrid)
+    REAL(rstd),INTENT(IN) :: lat(ngrid)
+    REAL(rstd),INTENT(OUT) :: phis(ngrid)
+    REAL(rstd),INTENT(OUT) :: ps(ngrid)
+    REAL(rstd),INTENT(OUT) :: temp(ngrid,llm)
+    REAL(rstd),INTENT(OUT) :: ulon(ngrid,llm)
+    REAL(rstd),INTENT(OUT) :: ulat(ngrid,llm)
+    REAL(rstd),INTENT(OUT) :: q(ngrid,llm,nqtot)
+    
+    INTEGER :: l,ij
+    REAL(rstd) :: etal, etavl, etas, etavs, sinlat, coslat, &
+         Y, Tave, T, phis_ave, vort, r2, utot, &
+         dthetaodeta_ave, dthetaodeta, dthetaodlat, duodeta, K, r
+    
+    etas=ap(1)/preff+bp(1)
+    etavs=(etas-eta0)*Pi/2
+    phis_ave=T0*g/Gamma*(1-etas**(Rd*Gamma/g))
+    
+    DO ij=1,ngrid
+       sinlat=SIN(lat(ij))
+       coslat=COS(lat(ij))
+       phis(ij)=phis_ave+u0*cos(etavs)**1.5*( (-2*sinlat**6 * (coslat**2+1./3) + 10./63 )*u0*cos(etavs)**1.5  &
+            +(8./5*coslat**3 * (sinlat**2 + 2./3) - Pi/4)*radius*Omega )
+       ps(ij)=ps0
+    ENDDO
+    
+    DO l=ll_begin,ll_end
+       etal = 0.5 *( ap(l)/preff+bp(l) + ap(l+1)/preff+bp(l+1) )
+       etavl=(etal-eta0)*Pi/2
+       
+       Tave=T0*etal**(Rd*Gamma/g)
+       dthetaodeta_ave = T0 *( Rd*Gamma/g - kappa)* etal**(Rd*Gamma/g-kappa-1)
+       IF (etat>etal) THEN
+          Tave=Tave+DeltaT*(etat-etal)**5
+          dthetaodeta_ave = dthetaodeta_ave - DeltaT * ( 5*(etat-etal)**4 * etal**(-kappa)  &
+               + kappa * (etat-etal)**5 * etal**(-kappa-1))
+       END IF
+       
+       DO ij=1,ngrid
+          sinlat=SIN(lat(ij))
+          coslat=COS(lat(ij))
+          
+          K=sin(latc)*sinlat+cos(latc)*coslat*cos(lon(ij)-lonc)
+          r=radius*acos(K)
+          utot=u0*cos(etavl)**1.5*sin(2*lat(ij))**2 + up0*exp(-(r/(0.1*radius))**2)
+          ulon(ij,l) = utot
+          ulat(ij,l) = 0.
+          Y = ((-2*sinlat**6*(coslat**2+1./3)+10./63)*2*u0*cos(etavl)**1.5     &
+               + (8./5*coslat**3*(sinlat**2+2./3)-Pi/4)*radius*Omega)
+          T = Tave + 0.75*(etal*Pi*u0/Rd)*sin(etavl)*cos(etavl)**0.5 * Y
+          temp(ij,l)=T
+          
+          q(ij,l,1)=q0*exp(-(lat(ij)/latw)**4)*exp(-((etal-1)*preff/pw)**2)
+          q(ij,l,2)=0.
+          q(ij,l,3)=0.
+          CALL initial_value_Terminator(lat(ij),lon(ij),q(ij,l,4),q(ij,l,5))
+       END DO
+    END DO
+    
+  END SUBROUTINE compute_etat0
+  
+END MODULE etat0_dcmip2016_baroclinic_wave_mod
