Index: codes/icosagcm/trunk/src/etat0.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0.f90	(revision 202)
+++ codes/icosagcm/trunk/src/etat0.f90	(revision 203)
@@ -15,6 +15,8 @@
     USE mpipara, ONLY : is_mpi_root
     USE disvert_mod
+    ! New interface
+    USE etat0_dcmip5_mod, ONLY : getin_etat0_dcmip5=>getin_etat0
+    ! Old interface
     USE etat0_williamson_mod, ONLY : etat0_williamson_new
-    USE etat0_jablonowsky06_mod, ONLY : etat0_jablonowsky06=>etat0
     USE etat0_academic_mod, ONLY : etat0_academic=>etat0  
     USE etat0_dcmip1_mod, ONLY : etat0_dcmip1=>etat0
@@ -22,5 +24,4 @@
     USE etat0_dcmip3_mod, ONLY : etat0_dcmip3=>etat0  
     USE etat0_dcmip4_mod, ONLY : etat0_dcmip4=>etat0  
-    USE etat0_dcmip5_mod, ONLY : etat0_dcmip5=>etat0  
     USE etat0_heldsz_mod, ONLY : etat0_heldsz=>etat0  
     USE dynetat0_gcm_mod, ONLY : dynetat0_start=>etat0 
@@ -51,13 +52,18 @@
     
     SELECT CASE (TRIM(etat0_type))
+       !------------------- New interface ---------------------
     CASE ('isothermal')
        CALL getin_etat0_isothermal
        CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q)
+    CASE ('jablonowsky06')
+       CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q)
+     CASE ('dcmip5')
+        CALL getin_etat0_dcmip5
+        CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q)
+        
+        !------------------- Old interface --------------------
     CASE ('williamson91.6')
        init_mass=.FALSE.
        CALL etat0_williamson_new(f_phis,f_mass,f_theta_rhodz,f_u, f_q)
-    CASE ('jablonowsky06')
-!       CALL etat0_jablonowsky06(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
-       CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q)
     CASE ('academic')
        CALL etat0_academic(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
@@ -79,6 +85,4 @@
         END IF
        CALL etat0_dcmip4(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
-     CASE ('dcmip5')
-       CALL etat0_dcmip5(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
      CASE ('readnf_start') 
 	  print*,"readnf_start used"	
@@ -135,5 +139,4 @@
       q=f_q(ind)
       CALL compute_etat0_collocated(ps,mass, phis, theta_rhodz, u, q)
-
     ENDDO
   END SUBROUTINE etat0_collocated
@@ -143,5 +146,6 @@
     USE theta2theta_rhodz_mod
     USE wind_mod
-    USE etat0_jablonowsky06_mod, ONLY : compute_jablonowsky06 => compute_etat0_new
+    USE etat0_jablonowsky06_mod, ONLY : compute_jablonowsky06 => compute_etat0
+    USE etat0_dcmip5_mod, ONLY : compute_dcmip5 => compute_etat0
     IMPLICIT NONE
     REAL(rstd),INTENT(INOUT) :: ps(iim*jjm)
@@ -189,4 +193,7 @@
        CALL compute_jablonowsky06(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i)
        CALL compute_jablonowsky06(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e)
+    CASE('dcmip5')
+       CALL compute_dcmip5(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q)
+       CALL compute_dcmip5(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_dcmip5.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0_dcmip5.f90	(revision 202)
+++ codes/icosagcm/trunk/src/etat0_dcmip5.f90	(revision 203)
@@ -16,235 +16,72 @@
   REAL(rstd),PARAMETER :: epsilon=1e-25
   REAL(rstd),PARAMETER :: Rd=287
-  REAL(rstd),SAVE :: lonc
-!$OMP THREADPRIVATE(lonc)  
-  REAL(rstd),SAVE :: latc
-!$OMP THREADPRIVATE(latc)  
-  TYPE(t_field),POINTER :: f_out_i(:)
-  REAL(rstd),POINTER,SAVE :: out_i(:,:)
-!$OMP THREADPRIVATE(out_i)  
+  
+  REAL(rstd), PARAMETER :: lonc=Pi, latc=Pi/18, &
+       Tv0=T0*(1+0.608*q0), Tvt=Tv0-Gamma*zt
 
-  PUBLIC  etat0
+  INTEGER, SAVE :: dcmip5_testcase
+
+  PUBLIC getin_etat0, compute_etat0
+
 CONTAINS
   
-   
-    
-  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
-  USE icosa
-  IMPLICIT NONE
-    TYPE(t_field),POINTER :: f_ps(:)
-    TYPE(t_field),POINTER :: f_phis(:)
-    TYPE(t_field),POINTER :: f_theta_rhodz(:)
-    TYPE(t_field),POINTER :: f_u(:)
-    TYPE(t_field),POINTER :: f_q(:)
-  
-    REAL(rstd),POINTER :: ps(:)
-    REAL(rstd),POINTER :: phis(:)
-    REAL(rstd),POINTER :: theta_rhodz(:,:)
-    REAL(rstd),POINTER :: u(:,:)
-    REAL(rstd),POINTER :: q(:,:,:)
-    INTEGER :: ind
- 
-    CALL allocate_field(f_out_i,field_t,type_real,llm)
- 
-    DO ind=1,ndomain
-      IF (.NOT. assigned_domain(ind)) CYCLE
-      CALL swap_dimensions(ind)
-      CALL swap_geometry(ind)
-      ps=f_ps(ind)
-      phis=f_phis(ind)
-      theta_rhodz=f_theta_rhodz(ind)
-      u=f_u(ind)
-      q=f_q(ind)
-      out_i=f_out_i(ind)
-      CALL compute_etat0_dcmip5(ps, phis, theta_rhodz, u, q)
-    ENDDO
-    CALL writefield("out_i",f_out_i)
-    CALL deallocate_field(f_out_i)
-  END SUBROUTINE etat0
-  
-  SUBROUTINE compute_etat0_dcmip5(ps, phis, theta_rhodz, ue, q)
-  USE icosa
-  USE pression_mod
-  USE wind_mod
-  USE theta2theta_rhodz_mod
-  USE exner_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) :: ue(3*iim*jjm,llm)
-  REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm,nqtot)
-  
-  INTEGER :: i,j,l,ij
-  INTEGER :: testcase
-  REAL(rstd) :: Tv0, Tvt 
-  REAL(rstd) :: r 
-  REAL(rstd) :: z(iim*jjm,llm),zz 
-  REAL(rstd) :: p(iim*jjm,llm+1) 
-  REAL(rstd) :: Tave,Tp 
-  REAL(rstd) :: T(iim*jjm,llm) 
-  REAL(rstd) :: vt 
-  REAL(rstd) :: d1,d2,d 
-  REAL(rstd) :: ulon(3*iim*jjm,llm) 
-  REAL(rstd) :: ulat(3*iim*jjm,llm)
-  REAL(rstd) :: lon,lat
-  REAL(rstd) :: pks(iim*jjm) 
-  REAL(rstd) :: pk(iim*jjm,llm) 
+  SUBROUTINE getin_etat0
+    dcmip5_testcase=1
+    CALL getin("dcmip5_testcase",dcmip5_testcase)
+  END SUBROUTINE getin_etat0
 
+  SUBROUTINE compute_etat0(ngrid,lon,lat, phis, ps, Temp, ulon, ulat, q)
+    USE disvert_mod
+    IMPLICIT NONE
+    INTEGER, INTENT(IN)    :: ngrid
+    REAL(rstd),INTENT(IN)  :: lon(ngrid)
+    REAL(rstd),INTENT(IN)  :: lat(ngrid)
+    REAL(rstd),INTENT(OUT) :: ps(ngrid)
+    REAL(rstd),INTENT(OUT) :: phis(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)
 
-    latc=Pi/18
-    lonc=Pi 
+    INTEGER :: l,ij
+    REAL(rstd) :: r, d, d1,d2, zz, Tave, Tp, vt, aa, bb
 
-    Tv0=T0*(1+0.608*q0)
-    Tvt=Tv0-Gamma*zt
-  
-    testcase=1
-    CALL getin("dcmip5_testcase",testcase)
-    
     phis(:)=0.
-   
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
-        CALL xyz2lonlat(xyz_i(ij,:),lon,lat)
-        r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc))
-        ps(ij)=pb-DeltaP*exp(-(r/rp)**1.5)
-      ENDDO
-    ENDDO
-    
-    CALL compute_pression(ps,p,0)
+
+    DO ij=1,ngrid
+       r=radius*arc(lonc,latc, lon(ij),lat(ij))
+       ps(ij)=pb-DeltaP*exp(-(r/rp)**1.5)
+    END DO
 
     DO l=1,llm
-      DO j=jj_begin,jj_end
-        DO i=ii_begin,ii_end
-          ij=(j-1)*iim+i
-          CALL xyz2lonlat(xyz_i(ij,:),lon,lat)
-          r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc))
-          CALL compute_z(0.5*(p(ij,l)+p(ij,l+1)),ps(ij),r,z(ij,l))
-        ENDDO
-      ENDDO
-    ENDDO
-    out_i=z
-    
-    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
-          IF (z(ij,l)<=zt) THEN
-            q(ij,l,1)=q0*exp(-z(ij,l)/zq1)*exp(-(z(ij,l)/zq2)**2)
+       aa=.5*(ap(l)+ap(l+1))
+       bb=.5*(bp(l)+bp(l+1))
+       DO ij=1,ngrid
+          r=radius*arc(lonc,latc, lon(ij),lat(ij))
+          CALL compute_z(aa+bb*ps(ij),ps(ij),r,zz)
+          IF (zz<=zt) THEN
+             q(ij,l,1)=q0*exp(-zz/zq1)*exp(-(zz/zq2)**2)
+             Tave=Tv0-Gamma*zz
+             Tp=(Tv0-Gamma*zz) * ( ( 1 + 2*Rd*(Tv0-Gamma*zz)*zz /                          &
+                  (g*zp**2*(1-pb/Deltap*exp((r/rp)**1.5+(zz/zp)**2))))**(-1) - 1) 
+             vt=-omega*sin(latc)*r+sqrt( (omega*sin(latc)*r)**2 -  1.5*(r/rp)**1.5 * (Tv0-Gamma*zz)*Rd    &
+                  / (1 + 2*Rd*(Tv0-Gamma*zz)*zz/(g*zp**2)        &
+                  - pb/Deltap * exp((r/rp)**1.5 + (zz/zp)**2 )))
           ELSE
-            q(ij,l,1)=qt
+             q(ij,l,1)=qt
+             Tave=Tvt
+             Tp=0
+             vt=0
           ENDIF
-        ENDDO
-      ENDDO
-    ENDDO
+          Temp(ij,l)=Tave+Tp 
 
-    
-     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 xyz2lonlat(xyz_i(ij,:),lon,lat)
-          r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc))
-          IF (z(ij,l)<=zt) THEN
-            Tave=Tv0-Gamma*z(ij,l)
-            Tp=(Tv0-Gamma*z(ij,l)) * ( ( 1 + 2*Rd*(Tv0-Gamma*z(ij,l))*z(ij,l) /                           &
-                                            (g*zp**2*(1-pb/Deltap*exp((r/rp)**1.5+(z(ij,l)/zp)**2))))**(-1) - 1) 
-          ELSE
-            Tave=Tvt
-            Tp=0
-          ENDIF
-          T(ij,l)=Tave+Tp 
-          out_i(ij,l)=Tave+Tp
-        ENDDO
-      ENDDO
-    ENDDO
-
-    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 xyz2lonlat(xyz_e(ij+u_right,:),lon,lat)
-          zz=0.5*(z(ij,l)+z(ij+t_right,l))
-          r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc))
-          IF (zz<=zt) THEN
-            vt=-omega*sin(latc)*r+sqrt( (omega*sin(latc)*r)**2 -  1.5*(r/rp)**1.5 * (Tv0-Gamma*zz)*Rd    &
-                                                                 / (1 + 2*Rd*(Tv0-Gamma*zz)*zz/(g*zp**2)        &
-                                                                      - pb/Deltap * exp((r/rp)**1.5 + (zz/zp)**2 )))
-          ELSE
-            vt = 0
-          ENDIF
-          d1=sin(latc)*cos(lat)-cos(latc)*sin(lat)*cos(lon-lonc)
-          d2=cos(latc)*sin(lon-lonc)
+          d1=sin(latc)*cos(lat(ij))-cos(latc)*sin(lat(ij))*cos(lon(ij)-lonc)
+          d2=cos(latc)*sin(lon(ij)-lonc)
           d=MAX(1e-25,sqrt(d1**2+d2**2))
           ulon(ij+u_right,l)=vt*d1/d
           ulat(ij+u_right,l)=vt*d2/d
-
-
-          
-          CALL xyz2lonlat(xyz_e(ij+u_lup,:),lon,lat)
-          zz=0.5*(z(ij,l)+z(ij+t_lup,l))
-          r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc))
-          IF (zz<=zt) THEN
-            vt=-omega*sin(latc)*r+sqrt( (omega*sin(latc)*r)**2 -  1.5*(r/rp)**1.5 * (Tv0-Gamma*zz)*Rd    &
-                                                                 / (1 + 2*Rd*(Tv0-Gamma*zz)*zz/(g*zp**2)        &
-                                                                      - pb/Deltap * exp((r/rp)**1.5 + (zz/zp)**2 )))
-
-         ELSE
-            vt = 0
-          ENDIF
-          d1=sin(latc)*cos(lat)-cos(latc)*sin(lat)*cos(lon-lonc)
-          d2=cos(latc)*sin(lon-lonc)
-          d=MAX(1e-25,sqrt(d1**2+d2**2))
-          ulon(ij+u_lup,l)=vt*d1/d
-          ulat(ij+u_lup,l)=vt*d2/d
-          
-          
-          
-          
-          CALL xyz2lonlat(xyz_e(ij+u_ldown,:),lon,lat)
-          zz=0.5*(z(ij,l)+z(ij+t_ldown,l))
-          r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc))
-          IF (zz<=zt) THEN
-            vt=-omega*sin(latc)*r+sqrt( (omega*sin(latc)*r)**2 -  1.5*(r/rp)**1.5 * (Tv0-Gamma*zz)*Rd    &
-                                                                 / (1 + 2*Rd*(Tv0-Gamma*zz)*zz/(g*zp**2)        &
-                                                                      - pb/Deltap * exp((r/rp)**1.5 + (zz/zp)**2 )))
-
-          ELSE
-            vt = 0
-          ENDIF
-          d1=sin(latc)*cos(lat)-cos(latc)*sin(lat)*cos(lon-lonc)
-          d2=cos(latc)*sin(lon-lonc)
-          d=MAX(1e-25,sqrt(d1**2+d2**2))
-          ulon(ij+u_ldown,l)=vt*d1/d
-          ulat(ij+u_ldown,l)=vt*d2/d
-          
-
-
-          CALL xyz2lonlat(xyz_i(ij,:),lon,lat)
-          zz=z(ij,l)
-          r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc))
-          IF (zz<=zt) THEN
-            vt=-omega*sin(latc)*r+sqrt( (omega*sin(latc)*r)**2 -  1.5*(r/rp)**1.5 * (Tv0-Gamma*zz)*Rd    &
-                                                                 / (1 + 2*Rd*(Tv0-Gamma*zz)*zz/(g*zp**2)        &
-                                                                      - pb/Deltap * exp((r/rp)**1.5 + (zz/zp)**2 )))
-          ELSE
-            vt = 0
-          ENDIF
-          d1=sin(latc)*cos(lat)-cos(latc)*sin(lat)*cos(lon-lonc)
-          d2=cos(latc)*sin(lon-lonc)
-          d=MAX(1e-25,sqrt(d1**2+d2**2))
-!          ulon(ij+u_ldown,l)=vt*d1/d
-!          ulat(ij+u_ldown,l)=vt*d2/d
-!          out_i(ij,l)=sqrt((vt*d1/d)**2+(vt*d2/d)**2)
-        ENDDO
-      ENDDO
+       END DO
     ENDDO
-    
-    CALL compute_wind_perp_from_lonlat_compound(ulon, ulat, ue)
-    CALL compute_temperature2theta_rhodz(ps,T,theta_rhodz,0)
-
-  END SUBROUTINE compute_etat0_dcmip5
+  END SUBROUTINE compute_etat0
   
   SUBROUTINE compute_z(pmodel,ps,r,z)
@@ -257,13 +94,10 @@
      REAL(rstd),INTENT(OUT) :: z
      
-     REAL(rstd) :: Tv0,Tvt,pt
-     REAL(rstd) :: pave, pp, dpdz
+     REAL(rstd) :: pt, pave, pp, dpdz
      REAL(rstd) :: znp1
      REAL(rstd) :: epsilon
      REAL(rstd) :: p
-     INTEGER  :: n
-          
-     Tv0=T0/(1+0.608*q0)
-     Tvt=Tv0-Gamma*zt
+     INTEGER  :: n     
+
      pt=pb*(Tvt/Tv0)**(g/(Rd*Gamma))
 
@@ -274,5 +108,4 @@
      ENDIF
      
-!     PRINT*,pmodel,pt,z,zt
      IF (z>zt .OR. r>1000000) return
      
@@ -297,6 +130,4 @@
        
        epsilon=ABS( (znp1-z)/znp1 )
-!       PRINT *,"epsilon",epsilon,r,z,znp1
-!       PRINT *,"----",n,pmodel,ps,pave,dpdz
        z=znp1
        n=n+1
Index: codes/icosagcm/trunk/src/etat0_jablonowsky06.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0_jablonowsky06.f90	(revision 202)
+++ codes/icosagcm/trunk/src/etat0_jablonowsky06.f90	(revision 203)
@@ -11,242 +11,10 @@
   REAL(rstd),PARAMETER :: Gamma=0.005
   REAL(rstd),PARAMETER :: up0=1
-  PUBLIC  test_etat0_jablonowsky06, etat0, compute_etat0_jablonowsky06, compute_etat0_new
+
+  PUBLIC compute_etat0
+
 CONTAINS
   
-  SUBROUTINE test_etat0_jablonowsky06
-  USE icosa
-  USE kinetic_mod
-  USE pression_mod
-  USE exner_mod
-  USE geopotential_mod
-  USE vorticity_mod
-  IMPLICIT NONE
-    TYPE(t_field),POINTER,SAVE :: f_ps(:)
-    TYPE(t_field),POINTER,SAVE :: f_phis(:)
-    TYPE(t_field),POINTER,SAVE :: f_theta_rhodz(:)
-    TYPE(t_field),POINTER,SAVE :: f_u(:)
-    TYPE(t_field),POINTER,SAVE :: f_Ki(:)
-    TYPE(t_field),POINTER,SAVE :: f_temp(:)
-    TYPE(t_field),POINTER,SAVE :: f_p(:)
-    TYPE(t_field),POINTER,SAVE :: f_pks(:)
-    TYPE(t_field),POINTER,SAVE :: f_pk(:)
-    TYPE(t_field),POINTER,SAVE :: f_phi(:)
-    TYPE(t_field),POINTER,SAVE :: f_vort(:)
-    TYPE(t_field),POINTER,SAVE :: f_q(:)
-  
-    REAL(rstd),POINTER :: Ki(:,:)
-    REAL(rstd),POINTER :: temp(:)
-    INTEGER :: ind
-    
-    
-    CALL allocate_field(f_ps,field_t,type_real)
-    CALL allocate_field(f_phis,field_t,type_real)
-    CALL allocate_field(f_theta_rhodz,field_t,type_real,llm)
-    CALL allocate_field(f_u,field_u,type_real,llm)
-    CALL allocate_field(f_p,field_t,type_real,llm+1)
-    CALL allocate_field(f_Ki,field_t,type_real,llm)
-    CALL allocate_field(f_pks,field_t,type_real)
-    CALL allocate_field(f_pk,field_t,type_real,llm)
-    CALL allocate_field(f_phi,field_t,type_real,llm)
-    CALL allocate_field(f_temp,field_t,type_real)
-    CALL allocate_field(f_vort,field_z,type_real,llm)
-    
-    CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q)
-
-    CALL kinetic(f_u,f_Ki)
-    CALL vorticity(f_u,f_vort)
-
-    CALL pression(f_ps,f_p)     
-    CALL exner(f_ps,f_p,f_pks,f_pk)  
-    CALL geopotential(f_phis,f_pks,f_pk,f_theta_rhodz,f_phi)
-
-    CALL writefield('ps',f_ps)
-    CALL writefield('phis',f_phis)
-    CALL writefield('theta',f_theta_rhodz)
-    CALL writefield('f_phi',f_phi)
-
-    CALL writefield('Ki',f_Ki)
-    CALL writefield('vort',f_vort)
-    
-  END SUBROUTINE test_etat0_jablonowsky06
-   
-    
-  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
-  USE icosa
-  IMPLICIT NONE
-    TYPE(t_field),POINTER :: f_ps(:)
-    TYPE(t_field),POINTER :: f_phis(:)
-    TYPE(t_field),POINTER :: f_theta_rhodz(:)
-    TYPE(t_field),POINTER :: f_u(:)
-    TYPE(t_field),POINTER :: f_q(:)
-  
-    REAL(rstd),POINTER :: ps(:)
-    REAL(rstd),POINTER :: phis(:)
-    REAL(rstd),POINTER :: theta_rhodz(:,:)
-    REAL(rstd),POINTER :: u(:,:)
-    REAL(rstd),POINTER :: q(:,:,:)
-    INTEGER :: ind
-    
-    DO ind=1,ndomain
-      IF (.NOT. assigned_domain(ind)) CYCLE
-      CALL swap_dimensions(ind)
-      CALL swap_geometry(ind)
-      ps=f_ps(ind)
-      phis=f_phis(ind)
-      theta_rhodz=f_theta_rhodz(ind)
-      u=f_u(ind)
-      q=f_q(ind)
-      q=0
-      CALL compute_etat0_jablonowsky06(ps, phis, theta_rhodz, u)
-    ENDDO
-
-  END SUBROUTINE etat0
-  
-  SUBROUTINE compute_etat0_jablonowsky06(ps, phis, theta_rhodz, u)
-  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)
-  
-  INTEGER :: i,j,l,ij
-  REAL(rstd) :: theta(iim*jjm,llm)
-  REAL(rstd) :: eta(llm)
-  REAL(rstd) :: etav(llm)
-  REAL(rstd) :: etas, etavs
-  REAL(rstd) :: lon,lat
-  REAL(rstd) :: ulon(3)
-  REAL(rstd) :: ep(3), norm_ep
-  REAL(rstd) :: Tave, T
-  REAL(rstd) :: phis_ave
-  REAL(rstd) :: V0(3)
-  REAL(rstd) :: r2
-  REAL(rstd) :: utot
-  REAL(rstd) :: lonx,latx
- 
-    DO l=1,llm
-      eta(l)= 0.5 *( ap(l)/ps0+bp(l) + ap(l+1)/ps0+bp(l+1) )
-      etav(l)=(eta(l)-eta0)*Pi/2
-    ENDDO
-    etas=ap(1)+bp(1)
-    etavs=(etas-eta0)*Pi/2
-
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
-        ps(ij)=ps0
-      ENDDO
-    ENDDO
-    
-    
-    CALL lonlat2xyz(Pi/9,2*Pi/9,V0)
-
-    u(:,:)=1e10      
-    DO l=1,llm
-      DO j=jj_begin,jj_end
-        DO i=ii_begin,ii_end
-          ij=(j-1)*iim+i
-          
-          CALL xyz2lonlat(xyz_e(ij+u_right,:)/radius,lon,lat)
-          CALL cross_product2(V0,xyz_e(ij+u_right,:)/radius,ep)
-          r2=(asin(sqrt(sum(ep*ep))))**2
-          utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-r2/0.01)
-
-          ulon(1) = -sin(lon) * utot
-          ulon(2) = cos(lon)  * utot
-          ulon(3) = 0         * utot
-          
-                    
-          CALL cross_product2(xyz_v(ij+z_rdown,:)/radius,xyz_v(ij+z_rup,:)/radius,ep)
-          norm_ep=sqrt(sum(ep(:)**2))
-          IF (norm_ep>1e-30) THEN
-            ep=-ep/norm_ep*ne(ij,right)
-            u(ij+u_right,l)=sum(ep(:)*ulon(:))
-          ENDIF
-
-
-          
-          CALL xyz2lonlat(xyz_e(ij+u_lup,:)/radius,lon,lat)
-          CALL cross_product2(V0,xyz_e(ij+u_lup,:)/radius,ep)
-          r2=(asin(sqrt(sum(ep*ep))))**2
-          utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-r2/0.01)
-          ulon(1) = -sin(lon) * utot
-          ulon(2) = cos(lon)  * utot
-          ulon(3) = 0         * utot
-          
-                    
-          CALL cross_product2(xyz_v(ij+z_up,:)/radius,xyz_v(ij+z_lup,:)/radius,ep)
-          norm_ep=sqrt(sum(ep(:)**2))
-          IF (norm_ep>1e-30) THEN
-            ep=-ep/norm_ep*ne(ij,lup)
-            u(ij+u_lup,l)=sum(ep(:)*ulon(:))
-          ENDIF
-
-
-
-
-
-                    
-          CALL xyz2lonlat(xyz_e(ij+u_ldown,:)/radius,lon,lat)
-          CALL cross_product2(V0,xyz_e(ij+u_ldown,:)/radius,ep)
-          r2=(asin(sqrt(sum(ep*ep))))**2
-          utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-r2/0.01)
-          ulon(1) = -sin(lon) * utot
-          ulon(2) = cos(lon)  * utot
-          ulon(3) = 0         * utot
-          
-                    
-          CALL cross_product2(xyz_v(ij+z_ldown,:)/radius,xyz_v(ij+z_down,:)/radius,ep)
-          norm_ep=sqrt(sum(ep(:)**2))
-          IF (norm_ep>1e-30) THEN
-            ep=-ep/norm_ep*ne(ij,ldown)
-            u(ij+u_ldown,l)=sum(ep(:)*ulon(:))
-          ENDIF          
-
-       ENDDO
-     ENDDO
-    ENDDO 
-     
-     
-     DO l=1,llm
-       Tave=T0*eta(l)**(Rd*Gamma/g)
-       IF (etat>eta(l)) Tave=Tave+DeltaT*(etat-eta(l))**5
-       DO j=jj_begin,jj_end
-         DO i=ii_begin,ii_end
-           ij=(j-1)*iim+i
-           CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat)
-             
-            T=Tave+ 0.75*(eta(l)*Pi*u0/Rd)*sin(etav(l))*cos(etav(l))**0.5              &
-                             * ( (-2*sin(lat)**6*(cos(lat)**2+1./3)+10./63)*2*u0*cos(etav(l))**1.5 &
-                                + (8./5*cos(lat)**3*(sin(lat)**2+2./3)-Pi/4)*radius*Omega)
-            
-            theta(ij,l)=T*eta(l)**(-kappa)
-
-          ENDDO
-       ENDDO
-     ENDDO
-     
-     
-     phis_ave=T0*g/Gamma*(1-etas**(Rd*Gamma/g))
-     DO j=jj_begin,jj_end
-       DO i=ii_begin,ii_end
-         ij=(j-1)*iim+i
-         CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat)
-         phis(ij)=phis_ave+u0*cos(etavs)**1.5*( (-2*sin(lat)**6 * (cos(lat)**2+1./3) + 10./63 )*u0*cos(etavs)**1.5  &
-                                           +(8./5*cos(lat)**3 * (sin(lat)**2 + 2./3) - Pi/4)*radius*Omega )
-       ENDDO
-     ENDDO
-
-    CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,0)
-
-  END SUBROUTINE compute_etat0_jablonowsky06
-
-  SUBROUTINE compute_etat0_new(ngrid,lat,lon, phis, ps, temp, ulon, ulat)
+  SUBROUTINE compute_etat0(ngrid,lon,lat, phis, ps, temp, ulon, ulat)
     USE disvert_mod
     IMPLICIT NONE  
@@ -283,5 +51,5 @@
        IF (etat>eta(l)) Tave=Tave+DeltaT*(etat-eta(l))**5
        DO ij=1,ngrid
-          r2 = arc(Pi/9,2*Pi/9, lon(ij),lat(ij))**2
+          r2 = arc(Pi/9.,2.*Pi/9., lon(ij),lat(ij))**2
           temp(ij,l) = Tave + 0.75*(eta(l)*Pi*u0/Rd)*sin(etav(l))*cos(etav(l))**0.5  &
                * ( (-2*sin(lat(ij))**6*(cos(lat(ij))**2+1./3)+10./63)*2*u0*cos(etav(l))**1.5 &
@@ -293,5 +61,5 @@
     
     
-  END SUBROUTINE compute_etat0_new
+  END SUBROUTINE compute_etat0
   
 END MODULE etat0_jablonowsky06_mod
Index: codes/icosagcm/trunk/src/vector.f90
===================================================================
--- codes/icosagcm/trunk/src/vector.f90	(revision 202)
+++ codes/icosagcm/trunk/src/vector.f90	(revision 203)
@@ -48,17 +48,7 @@
     END SUBROUTINE cross_product2
 
-    FUNCTION arc(lon1,lat1, lon2,lat2)
-      REAL(rstd) :: lon1,lat1, lon2,lat2, arc2, &
-           x1,x2,x3,y1,y2,y3,z1,z2,z3
-      x1=cos(lat1)*cos(lon1)
-      x2=cos(lat2)*cos(lon2)
-      y1=cos(lat1)*sin(lon1)
-      y2=cos(lat2)*sin(lon2)
-      z1=sin(lat1)
-      z2=sin(lat2)
-      x3=y1*z2-y2*z1
-      y3=z1*x2-z2*x1
-      z3=x1*y2-x2*y1
-      arc=ASIN(SQRT(x3*x3+y3*y3+z3*z3))
+    FUNCTION arc(lon,lat, lonc,latc)
+      REAL(rstd) :: lon,lat, lonc,latc, arc
+      arc=ACOS(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc))
     END FUNCTION arc
 
