Index: codes/icosagcm/trunk/src/etat0.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0.f90	(revision 168)
+++ codes/icosagcm/trunk/src/etat0.f90	(revision 170)
@@ -49,6 +49,6 @@
     CASE ('academic')
        CALL etat0_academic(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
-    CASE ('heldsz')
-	print*,"heldsz test case"
+    CASE ('held_suarez')
+       PRINT *,"Held & Suarez (1994) test case"
        CALL etat0_heldsz(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
     CASE ('dcmip1')
Index: codes/icosagcm/trunk/src/etat0_heldsz.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0_heldsz.f90	(revision 168)
+++ codes/icosagcm/trunk/src/etat0_heldsz.f90	(revision 170)
@@ -1,16 +1,26 @@
-    MODULE etat0_heldsz_mod
-         USE icosa
-      IMPLICIT NONE
-  REAL(rstd),ALLOCATABLE::knewt_t(:),kfrict(:)
-  REAL(rstd)::knewt_g
-  TYPE(t_field),POINTER :: f_tetarappel(:)
-  TYPE(t_field),POINTER :: f_clat(:)
-
+MODULE etat0_heldsz_mod
+  USE icosa
+  IMPLICIT NONE
+  PRIVATE
+
+  TYPE(t_field),POINTER :: f_theta_eq(:)
+  TYPE(t_field),POINTER :: f_theta(:)
+  TYPE(t_field),POINTER :: f_clat(:) ! FIXME, duplication
+
+  REAL(rstd),ALLOCATABLE :: knewt_t(:),kfrict(:)
+
+  LOGICAL, SAVE :: done=.FALSE.
+
+  REAL(rstd) :: teta0,ttp,delt_y,delt_z,eps
+  REAL(rstd) :: knewt_g, k_f,k_c_a,k_c_s
+
+  PUBLIC :: etat0, held_suarez
+  
 CONTAINS
-  
+
   SUBROUTINE test_etat0_heldsz
-  USE icosa
-  USE kinetic_mod
-  IMPLICIT NONE
+    USE icosa
+    USE kinetic_mod
+    IMPLICIT NONE
     TYPE(t_field),POINTER :: f_ps(:)
     TYPE(t_field),POINTER :: f_phis(:)
@@ -19,9 +29,8 @@
     TYPE(t_field),POINTER :: f_q(:)
     TYPE(t_field),POINTER :: f_Ki(:)
-  
+
     REAL(rstd),POINTER :: Ki(:,:)
     INTEGER :: ind
-    
-    
+
     CALL allocate_field(f_ps,field_t,type_real)
     CALL allocate_field(f_phis,field_t,type_real)
@@ -29,5 +38,5 @@
     CALL allocate_field(f_u,field_u,type_real,llm)
     CALL allocate_field(f_Ki,field_t,type_real,llm)
-    
+
     CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
     CALL kinetic(f_u,f_Ki)
@@ -36,9 +45,9 @@
     CALL writefield('theta',f_theta_rhodz)
   END SUBROUTINE test_etat0_heldsz
-   
-    
+
   SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
-  USE icosa
-  IMPLICIT NONE
+    USE icosa
+    USE theta2theta_rhodz_mod
+    IMPLICIT NONE
     TYPE(t_field),POINTER :: f_ps(:)
     TYPE(t_field),POINTER :: f_phis(:)
@@ -46,5 +55,5 @@
     TYPE(t_field),POINTER :: f_u(:)
     TYPE(t_field),POINTER :: f_q(:)
-  
+
     REAL(rstd),POINTER :: ps(:)
     REAL(rstd),POINTER :: phis(:)
@@ -52,148 +61,161 @@
     REAL(rstd),POINTER :: u(:,:)
     REAL(rstd),POINTER :: q(:,:,:)
+    REAL(rstd),POINTER :: clat(:) 
+    REAL(rstd),POINTER :: theta_eq(:,:) 
+    REAL(rstd),POINTER :: theta(:,:) 
+
     INTEGER :: ind
-    REAL(rstd),POINTER::clat(:) 
-    REAL(rstd),POINTER::tetarappel(:,:) 
-    
-    CALL allocate_field(f_tetarappel,field_t,type_real,llm)
-    CALL allocate_field(f_clat,field_t,type_real)
-    ALLOCATE(knewt_t(llm)); ALLOCATE( kfrict(llm)) 
-
+
+    CALL Init_Teq
     DO ind=1,ndomain
-      CALL swap_dimensions(ind)
-      CALL swap_geometry(ind)
-      ps=f_ps(ind)
-      phis=f_phis(ind)
-      theta_rhodz=f_theta_rhodz(ind)
-      tetarappel=f_tetarappel(ind) 
-      u=f_u(ind)
-      q=f_q(ind)
-      q=1e2
-      clat=f_clat(ind) 
-      CALL compute_etat0_heldsz(ps, phis, theta_rhodz, u,clat,tetarappel)
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+
+       theta_eq=f_theta_eq(ind) 
+       clat=f_clat(ind) 
+       CALL compute_Teq(clat,theta_eq) ! FIXME : already done by Init_Teq
+
+       ps=f_ps(ind)
+       phis=f_phis(ind)
+       u=f_u(ind)
+       ps(:)=1e5
+       phis(:)=0.
+       u(:,:)=0
+
+       theta_rhodz=f_theta_rhodz(ind)
+       theta=f_theta(ind)
+       CALL compute_etat0_heldsz(theta_eq,theta)
+       CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,1)
+       q=f_q(ind)
+       q(:,:,:)=1e2
     ENDDO
   END SUBROUTINE etat0
 
-  SUBROUTINE compute_etat0_heldsz(ps, phis, theta_rhodz, u,clat,tetarappel)
-  USE icosa
-  USE disvert_mod
-  USE pression_mod
-  USE exner_mod
-  USE geopotential_mod
-  USE theta2theta_rhodz_mod
-  IMPLICIT NONE  
-  REAL(rstd),INTENT(OUT) :: ps(iim*jjm)
-  REAL(rstd),INTENT(OUT) :: phis(iim*jjm)
-  REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
-  REAL(rstd),INTENT(OUT) :: u(3*iim*jjm,llm)
-  REAL(rstd),INTENT(OUT) :: clat(iim*jjm)
-  REAL(rstd),INTENT(OUT) :: tetarappel(iim*jjm,llm) 
-  
-  INTEGER :: i,j,l,ij
-  REAL(rstd) :: r
-  REAL(rstd) :: theta(iim*jjm,llm)
-  REAL(rstd) :: zsig
-  INTEGER    :: lsup
-  REAL(rstd) :: ddsin
-  REAL(rstd) :: lon,lat
-  REAL(rstd) :: p(iim*jjm,llm+1)
-  REAL(rstd) :: alpha(iim*jjm,llm),beta(iim*jjm,llm)
-  REAL(rstd) :: delta
-  REAL(rstd) :: pks(iim*jjm),pk(iim*jjm,llm)
-  REAL(rstd) :: phi(iim*jjm,llm)
-  REAL(rstd) :: x  
-  REAL(rstd) :: fact(3*iim*jjm)
-  REAL(rstd) :: ut(3*iim*jjm,llm)
-
-  REAL(rstd) :: teta0,ttp,delt_y,delt_z,eps
-  REAL(rstd) :: k_f,k_c_a,k_c_s
-  REAL(rstd) :: zz,ran1
-  REAL(rstd) :: tetastrat,tetajl(iim*jjm,llm) 
-  REAL(rstd) :: slat(iim*jjm) 
-!-------------choces of parametes and get it  
-     k_f=1.                !friction 
-     CALL getin('k_j',k_f)
-     k_f=1./(daysec*k_f)
-     k_c_s=4.  !cooling surface
-     CALL getin('k_c_s',k_c_s)
-     k_c_s=1./(daysec*k_c_s)
-     k_c_a=40. !cooling free atm
-     CALL getin('k_c_a',k_c_a)
-     k_c_a=1./(daysec*k_c_a)
-     ! Constants for Teta equilibrium profile
-     teta0=315.     ! mean Teta (S.H. 315K)
-     CALL getin('teta0',teta0)
-     ttp=200.       ! Tropopause temperature (S.H. 200K)
-     CALL getin('ttp',ttp)
-     eps=0.         ! Deviation to N-S symmetry(~0-20K)
-     CALL getin('eps',eps)
-     delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
-     CALL getin('delt_y',delt_y)
-     delt_z=10.     ! Vertical Gradient (S.H. 10K)
-     CALL getin('delt_z',delt_z)
-!-----------------------------------------------------------
-    knewt_g=k_c_a 
-    DO l=1,llm
-       zsig=ap(l)/preff+bp(l)
-       knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3) 
-       kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3) 
-    ENDDO
-         DO j=jj_begin-1,jj_end+1
-           DO i=ii_begin-1,ii_end+1
-             ij=(j-1)*iim+i
-             CALL xyz2lonlat(xyz_i(ij,:),lon,lat)
-	     clat(ij)=cos(lat) 
-	     slat(ij)=sin(lat) 
-            ENDDO
-         ENDDO
+  SUBROUTINE init_Teq
+    USE icosa
+    USE disvert_mod, ONLY : ap,bp
+    IMPLICIT NONE
+    REAL(rstd),POINTER :: clat(:) 
+    REAL(rstd),POINTER :: theta_eq(:,:) 
+    REAL(rstd) :: zsig
+    INTEGER :: ind, l
+
+    IF(.NOT.done) THEN
+       done = .TRUE.
+       
+       CALL allocate_field(f_theta,field_t,type_real,llm)
+       CALL allocate_field(f_theta_eq,field_t,type_real,llm)
+       CALL allocate_field(f_clat,field_t,type_real)
+       ALLOCATE(knewt_t(llm)); ALLOCATE( kfrict(llm)) 
+
+       k_f=1.                !friction 
+       CALL getin('k_j',k_f)
+       k_f=1./(daysec*k_f)
+       k_c_s=4.  !cooling surface
+       CALL getin('k_c_s',k_c_s)
+       k_c_s=1./(daysec*k_c_s)
+       k_c_a=40. !cooling free atm
+       CALL getin('k_c_a',k_c_a)
+       k_c_a=1./(daysec*k_c_a)
+       ! Constants for Teta equilibrium profile
+       teta0=315.     ! mean Teta (S.H. 315K)
+       CALL getin('teta0',teta0)
+       ttp=200.       ! Tropopause temperature (S.H. 200K)
+       CALL getin('ttp',ttp)
+       eps=0.         ! Deviation to N-S symmetry(~0-20K)
+       CALL getin('eps',eps)
+       delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
+       CALL getin('delt_y',delt_y)
+       delt_z=10.     ! Vertical Gradient (S.H. 10K)
+       CALL getin('delt_z',delt_z)
+
+       !-----------------------------------------------------------
+       knewt_g=k_c_a 
+       DO l=1,llm
+          zsig=ap(l)/preff+bp(l)
+          knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3) 
+          kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3) 
+       ENDDO
+
+       DO ind=1,ndomain
+          CALL swap_dimensions(ind)
+          CALL swap_geometry(ind)
+          clat=f_clat(ind)
+          theta_eq=f_theta_eq(ind)
+          CALL compute_Teq(clat,theta_eq)
+       ENDDO
+
+    ELSE
+       PRINT *, 'Init_Teq called twice'
+       CALL ABORT
+    END IF
+
+  END SUBROUTINE init_Teq
+
+  SUBROUTINE compute_Teq(clat,theta_eq)
+    USE icosa
+    USE disvert_mod
+    IMPLICIT NONE  
+    REAL(rstd),INTENT(OUT) :: clat(iim*jjm)
+    REAL(rstd),INTENT(OUT) :: theta_eq(iim*jjm,llm) 
+
+    REAL(rstd) :: lon, lat, r, zsig, ddsin, tetastrat, tetajl
+    REAL(rstd) :: slat(iim*jjm) 
+    INTEGER :: i,j,l,ij
+
+    DO j=jj_begin-1,jj_end+1
+       DO i=ii_begin-1,ii_end+1
+          ij=(j-1)*iim+i
+          CALL xyz2lonlat(xyz_i(ij,:),lon,lat)
+          clat(ij)=cos(lat) 
+          slat(ij)=sin(lat) 
+       ENDDO
+    ENDDO
 
     DO l=1,llm
        zsig=ap(l)/preff+bp(l)
        tetastrat=ttp*zsig**(-kappa)
-         DO j=jj_begin-1,jj_end+1
-           DO i=ii_begin-1,ii_end+1
+       DO j=jj_begin-1,jj_end+1
+          DO i=ii_begin-1,ii_end+1
              ij=(j-1)*iim+i
 	     ddsin=slat(ij) 
-	     tetajl(ij,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin &
-                -delt_z*(1.-ddsin*ddsin)*log(zsig)
-	     tetajl(ij,l)=MAX(tetajl(ij,l),tetastrat) 
-	     tetarappel(ij,l)=tetajl(ij,l) 
-	    ENDDO
-	  ENDDO
-	ENDDO 
-     
-     DO j=jj_begin-1,jj_end+1
-       DO i=ii_begin-1,ii_end+1
-         ij=(j-1)*iim+i
-         ps(ij)=100000.0
-         phis(ij)=0.0
-       ENDDO
-     ENDDO
-     
-
-    CALL compute_pression(ps,p,1)     
-    CALL compute_exner(ps,p,pks,pk,1)  
-	theta(:,:)=tetarappel(:,:) 
-    CALL compute_geopotential(phis,pks,pk,theta,phi,1)
-
-	u=0.0 !!wind 0 
-!============================================================
-     DO l=1,llm
-      DO j=jj_begin-1,jj_end+1
-        DO i=ii_begin-1,ii_end+1
-          ij=(j-1)*iim+i
-	  CALL RANDOM_NUMBER(r); r = 0.0 
-          theta(ij,l)=theta(ij,l)*(1.+0.0005*r)
-       ENDDO
-      ENDDO
-     ENDDO
-    CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,1)
+	     tetajl=teta0-delt_y*ddsin*ddsin+eps*ddsin &
+                  -delt_z*(1.-ddsin*ddsin)*log(zsig)
+	     theta_eq(ij,l)=MAX(tetajl,tetastrat) 
+          ENDDO
+       ENDDO
+    ENDDO
+  END SUBROUTINE compute_Teq
+
+  SUBROUTINE compute_etat0_heldsz(theta_eq, theta)
+    USE icosa
+    USE disvert_mod
+    USE pression_mod
+    USE exner_mod
+    USE geopotential_mod
+    USE theta2theta_rhodz_mod
+    IMPLICIT NONE  
+    REAL(rstd),INTENT(IN) :: theta_eq(iim*jjm,llm) 
+    REAL(rstd),INTENT(OUT) :: theta(iim*jjm,llm) 
+
+    REAL(rstd) :: r  ! random number
+    INTEGER :: i,j,l,ij
+
+    DO l=1,llm
+       DO j=jj_begin-1,jj_end+1
+          DO i=ii_begin-1,ii_end+1
+             ij=(j-1)*iim+i
+             CALL RANDOM_NUMBER(r); r = 0.0 
+             theta(ij,l)=theta_eq(ij,l)*(1.+0.0005*r)
+          ENDDO
+       ENDDO
+    ENDDO
 
   END SUBROUTINE compute_etat0_heldsz
 
 
-	SUBROUTINE held_saurez(f_ps,f_theta_rhodz,f_u) 
-	USE icosa
-	    IMPLICIT NONE
+  SUBROUTINE held_suarez(f_ps,f_theta_rhodz,f_u) 
+    USE icosa
+    IMPLICIT NONE
     TYPE(t_field),POINTER :: f_theta_rhodz(:)
     TYPE(t_field),POINTER :: f_u(:)
@@ -202,49 +224,52 @@
     REAL(rstd),POINTER :: u(:,:)
     REAL(rstd),POINTER :: ps(:)
-    REAL(rstd),POINTER :: tetarappel(:,:)
+    REAL(rstd),POINTER :: theta_eq(:,:)
+    REAL(rstd),POINTER :: theta(:,:)
     REAL(rstd),POINTER :: clat(:)
     INTEGER::ind
 
     DO ind=1,ndomain
-      CALL swap_dimensions(ind)
-      CALL swap_geometry(ind)
-      theta_rhodz=f_theta_rhodz(ind)
-      u=f_u(ind)
-      ps=f_ps(ind) 
-      tetarappel=f_tetarappel(ind) 
-      clat=f_clat(ind) 
-      CALL compute_heldsz(ps,theta_rhodz,u,clat,tetarappel) 
-    ENDDO
-        END SUBROUTINE held_saurez 
-
-	SUBROUTINE compute_heldsz(ps,theta_rhodz,u,clat,tetarappel) 
-	USE icosa
-	USE theta2theta_rhodz_mod
-           IMPLICIT NONE
-   REAL(rstd),INTENT(IN)::ps(iim*jjm) 
-   REAL(rstd),INTENT(INOUT) :: theta_rhodz(iim*jjm,llm)
-   REAL(rstd),INTENT(INOUT) :: u(3*iim*jjm,llm)
-   REAL(rstd)::theta(iim*jjm,llm) 
-   REAL(rstd),INTENT(IN)::tetarappel(iim*jjm,llm) 
-   REAL(rstd),INTENT(IN):: clat(iim*jjm) 
-   INTEGER :: i,j,l,ij
-
-      CALL compute_theta_rhodz2theta(ps,theta_rhodz,theta,1)
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       theta_rhodz=f_theta_rhodz(ind)
+       u=f_u(ind)
+       ps=f_ps(ind) 
+       theta_eq=f_theta_eq(ind) 
+       theta=f_theta(ind) 
+       clat=f_clat(ind) 
+       CALL compute_heldsz(ps,theta_eq,clat, theta_rhodz,u, theta) 
+    ENDDO
+  END SUBROUTINE held_suarez
+
+  SUBROUTINE compute_heldsz(ps,theta_eq,clat, theta_rhodz,u, theta) 
+    USE icosa
+    USE theta2theta_rhodz_mod
+    IMPLICIT NONE
+    REAL(rstd),INTENT(IN)    :: ps(iim*jjm) 
+    REAL(rstd),INTENT(IN)    :: theta_eq(iim*jjm,llm) 
+    REAL(rstd),INTENT(IN)    :: clat(iim*jjm) 
+    REAL(rstd),INTENT(INOUT) :: theta_rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(INOUT) :: u(3*iim*jjm,llm)
+    REAL(rstd),INTENT(OUT)   :: theta(iim*jjm,llm) 
+
+    INTEGER :: i,j,l,ij
+
+    CALL compute_theta_rhodz2theta(ps,theta_rhodz,theta,1)
     DO l=1,llm
-      DO j=jj_begin-1,jj_end+1
-        DO i=ii_begin-1,ii_end+1
-          ij=(j-1)*iim+i
-        theta(ij,l)=theta(ij,l) - dt*(theta(ij,l)-tetarappel(ij,l))* &
-        (knewt_g+knewt_t(l)*clat(ij)**4 )
-	ENDDO
-      ENDDO
-    ENDDO 
-       CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,1)
-
-	Do l=1,llm
-	u(:,l)=u(:,l)*(1.-dt*kfrict(l))
-	END DO
-
-	  END SUBROUTINE compute_heldsz
+       DO j=jj_begin-1,jj_end+1
+          DO i=ii_begin-1,ii_end+1
+             ij=(j-1)*iim+i
+             theta(ij,l)=theta(ij,l) - dt*(theta(ij,l)-theta_eq(ij,l))* &
+                  (knewt_g+knewt_t(l)*clat(ij)**4 )
+          ENDDO
+       ENDDO
+    ENDDO
+    CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,1)
+
+    Do l=1,llm
+       u(:,l)=u(:,l)*(1.-dt*kfrict(l))
+    END DO
+
+  END SUBROUTINE compute_heldsz
 
 END MODULE etat0_heldsz_mod
Index: codes/icosagcm/trunk/src/physics.f90
===================================================================
--- codes/icosagcm/trunk/src/physics.f90	(revision 168)
+++ codes/icosagcm/trunk/src/physics.f90	(revision 170)
@@ -1,5 +1,5 @@
 MODULE physics_mod
 
-  CHARACTER(LEN=255) :: physics_type="none"
+  CHARACTER(LEN=255) :: physics_type="automatic"
 
 
@@ -7,35 +7,35 @@
 
   SUBROUTINE init_physics
-  USE icosa
-  USE physics_dcmip_mod,init_physics_dcmip=>init_physics
-  USE physics_dry_mod
-  IMPLICIT NONE
-    
+    USE icosa
+    USE physics_dcmip_mod,init_physics_dcmip=>init_physics
+    USE physics_dry_mod
+    IMPLICIT NONE
+
     CALL getin("physics",physics_type)
-    
+
     SELECT CASE(TRIM(physics_type))
-      CASE ('none')
-    
-      CASE ('dcmip')
-        CALL init_physics_dcmip
+    CASE ('automatic')
 
-      CASE ('lmd')
-	CALL init_physics_dry
-      
-      CASE DEFAULT
-         PRINT*, 'Bad selector for variable physics init <',physics_type, &
-              '> options are <none>, <dcmip>,'
+    CASE ('dcmip')
+       CALL init_physics_dcmip
 
+    CASE ('dry')
+       CALL init_physics_dry
+
+    CASE DEFAULT
+       PRINT*, 'init_physics : Bad selector for variable physics <',TRIM(physics_type), &
+            '> options are <automatic>, <dcmip>, <dry>'
+       STOP
     END SELECT
-    
+
   END SUBROUTINE init_physics
 
   SUBROUTINE physics(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
-  USE icosa
-  USE physics_dry_mod
-  USE physics_dcmip_mod, physics_dcmip=>physics
-  USE etat0_mod
-  USE etat0_heldsz_mod
-  IMPLICIT NONE
+    USE icosa
+    USE physics_dry_mod
+    USE physics_dcmip_mod, physics_dcmip=>physics
+    USE etat0_mod
+    USE etat0_heldsz_mod
+    IMPLICIT NONE
     INTEGER, INTENT(IN)   :: it
     REAL(rstd),INTENT(IN)::jD_cur,jH_cur
@@ -46,30 +46,30 @@
     TYPE(t_field),POINTER :: f_q(:)
     LOGICAL:: firstcall,lastcall
-    
+
     SELECT CASE(TRIM(physics_type))
-      CASE ('none')
+    CASE ('automatic')
 
-	SELECT CASE(TRIM(etat0_type))
-	CASE('heldsz') 
-!	CALL transfert_request(f_ps,req_i1)
-!	CALL transfert_request(f_theta_rhodz,req_i1)
-!	CALL transfert_request(f_ue,req_e1_vect)
-!	CALL held_saurez(f_ps,f_theta_rhodz,f_ue) 
-	CASE DEFAULT
-	PRINT*,"NO PHYSICAL PACAKAGE USED" 
-	END SELECT
-    
-      CASE ('dcmip')
-        CALL physics_dcmip(it, f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
+       SELECT CASE(TRIM(etat0_type))
+       CASE('held_suarez')
+          !	CALL transfert_request(f_ps,req_i1)
+          !	CALL transfert_request(f_theta_rhodz,req_i1)
+          !	CALL transfert_request(f_ue,req_e1_vect)
+          CALL held_suarez(f_ps,f_theta_rhodz,f_ue) 
+       CASE DEFAULT
+          PRINT*,"NO PHYSICAL PACAKAGE USED" ! FIXME MPI
+       END SELECT
 
-      CASE ('dry')
-        CALL physics_dry(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
-      
-      CASE DEFAULT
-         PRINT*, 'Bad selector for variable physics <',physics_type, &
-              '> options are <none>, <dcmip>,'
-  STOP
+    CASE ('dcmip')
+       CALL physics_dcmip(it, f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
+
+    CASE ('dry')
+       CALL physics_dry(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
+
+    CASE DEFAULT
+       PRINT*, 'Bad selector for variable physics <',TRIM(physics_type), &
+            '> options are <automatic>, <dcmip>, <dry>'
+       STOP
     END SELECT
-    
+
   END SUBROUTINE physics
 
Index: codes/icosagcm/trunk/src/timeloop_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 168)
+++ codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 170)
@@ -303,11 +303,9 @@
 
 !----------------------------------------------------
-!    jD_cur = jD_ref + day_ini - day_ref + it/day_step
-!    jH_cur = jH_ref + start_time + mod(it,day_step)/float(day_step)
-!    jD_cur = jD_cur + int(jH_cur)
-!    jH_cur = jH_cur - int(jH_cur)
-!    CALL physics(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_u, f_q)
-
-!    CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q)
+    jD_cur = jD_ref + day_ini - day_ref + it/day_step
+    jH_cur = jH_ref + start_time + mod(it,day_step)/float(day_step)
+    jD_cur = jD_cur + int(jH_cur)
+    jH_cur = jH_cur - int(jH_cur)
+    CALL physics(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_u, f_q)
     ENDDO
 
