Index: codes/icosagcm/trunk/src/caldyn_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 285)
+++ codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 286)
@@ -328,14 +328,11 @@
   INTEGER :: ij
   DO ij=ij_begin_ext,ij_end_ext
-     CALL xyz2lonlat(xyz_e(ij+u_right,:),lon,lat)
-     ulon(ij+u_right)=a*omega*cos(lat)
+     ulon(ij+u_right)=a*omega*cos(lat_e(ij+u_right))
      ulat(ij+u_right)=0
 
-     CALL xyz2lonlat(xyz_e(ij+u_lup,:),lon,lat)
-     ulon(ij+u_lup)=a*omega*cos(lat)
+     ulon(ij+u_lup)=a*omega*cos(lat_e(ij+u_lup))
      ulat(ij+u_lup)=0
 
-     CALL xyz2lonlat(xyz_e(ij+u_ldown,:),lon,lat)
-     ulon(ij+u_ldown)=a*omega*cos(lat)
+     ulon(ij+u_ldown)=a*omega*cos(lat_e(ij+u_ldown))
      ulat(ij+u_ldown)=0
   END DO
Index: codes/icosagcm/trunk/src/check_conserve.f90
===================================================================
--- codes/icosagcm/trunk/src/check_conserve.f90	(revision 285)
+++ codes/icosagcm/trunk/src/check_conserve.f90	(revision 286)
@@ -256,11 +256,5 @@
 	  ij=(j-1)*iim+i
           IF (domain(ind)%own(i,j)) THEN
-!            CALL xyz2lonlat(xyz_i(ij,:),lon,lat) 
-! very bad inlining with intel compiler (>14 ?)
-             xyz(:)=xyz_i(ij,:)
-             xyz(:)=xyz(:)/sqrt(sum(xyz(:)**2))
-             lat=asin(xyz(3))
-             lon=atan2(xyz(2),xyz(1))
-  
+            lat=lat_i(ij)  
             ang=ang+rad*cos(lat)*masse(ij,l)*(ulon(ij,l)+ radomeg*cos(lat))
           END IF 
Index: codes/icosagcm/trunk/src/dissip_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/dissip_gcm.f90	(revision 285)
+++ codes/icosagcm/trunk/src/dissip_gcm.f90	(revision 286)
@@ -540,7 +540,6 @@
 
            nn = ij+shift_u
-           CALL xyz2lonlat(xyz_e(nn,:),lon,lat)
            zcoords = 1 ; hybrid_eta = .FALSE. ! use z instead of p or hyam/hybm
-           CALL test2_schaer_mountain(lon,lat,p,z,zcoords,hybrid_eta, &
+           CALL test2_schaer_mountain(lon_e(nn),lat_e(nn),p,z,zcoords,hybrid_eta, &
                 hyam,hybm,shear,ulon,ulat,w,t,phis,ps,rho,q)
 !           u3d = ulon*elon_e(nn,:) + ulat*elat_e(nn,:)
Index: codes/icosagcm/trunk/src/etat0.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0.f90	(revision 285)
+++ codes/icosagcm/trunk/src/etat0.f90	(revision 286)
@@ -158,12 +158,8 @@
     REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm,nqtot)
 
-    REAL(rstd) :: lon_i(iim*jjm)
-    REAL(rstd) :: lat_i(iim*jjm)
     REAL(rstd) :: temp_i(iim*jjm,llm)
     REAL(rstd) :: ulon_i(iim*jjm,llm)
     REAL(rstd) :: ulat_i(iim*jjm,llm)
 
-    REAL(rstd) :: lon_e(3*iim*jjm)
-    REAL(rstd) :: lat_e(3*iim*jjm)
     REAL(rstd) :: ps_e(3*iim*jjm)
     REAL(rstd) :: mass_e(3*iim*jjm,llm)
@@ -175,16 +171,4 @@
 
     INTEGER :: l,i,j,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 xyz2lonlat(xyz_i(ij,:)/radius,lon_i(ij),lat_i(ij))
-             CALL xyz2lonlat(xyz_e(ij+u_right,:)/radius,lon_e(ij+u_right),lat_e(ij+u_right))
-             CALL xyz2lonlat(xyz_e(ij+u_lup,:)/radius,lon_e(ij+u_lup),lat_e(ij+u_lup))
-             CALL xyz2lonlat(xyz_e(ij+u_ldown,:)/radius,lon_e(ij+u_ldown),lat_e(ij+u_ldown))
-          END DO
-       END DO
-    END DO
 
     SELECT CASE (TRIM(etat0_type))
Index: codes/icosagcm/trunk/src/etat0_academic.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0_academic.f90	(revision 285)
+++ codes/icosagcm/trunk/src/etat0_academic.f90	(revision 286)
@@ -109,6 +109,5 @@
            DO i=ii_begin-1,ii_end+1
              ij=(j-1)*iim+i
-             CALL xyz2lonlat(xyz_i(ij,:),lon,lat)
-             ddsin=sin(lat)-sin(pi/20.)
+             ddsin=sin(lat_i(ij))-sin(pi/20.)
              theta(ij,l)=300.*(1+1./18.*(1.-3.*ddsin*ddsin)+thetarappel)
             ENDDO
@@ -144,25 +143,18 @@
          ij=(j-1)*iim+i
 
-         CALL xyz2lonlat(xyz_e(ij+u_right,:),lon,lat)
+         lat=lat_e(ij+u_right)
          IF (abs(sin(lat))<1.e-4) lat=1e-4
          x=cos(lat) ; x=x*x ; x=x*x ; x=x*x
          fact(ij+u_right)=(1.-x)/(2.*omega*sin(lat))
 
-         CALL xyz2lonlat(xyz_e(ij+u_lup,:),lon,lat)
+         lat=lat_e(ij+u_lup)
          IF (abs(sin(lat))<1.e-4) lat=1e-4
          x=cos(lat) ; x=x*x ; x=x*x ; x=x*x
          fact(ij+u_lup)=(1.-x)/(2.*omega*sin(lat))
 
-         CALL xyz2lonlat(xyz_e(ij+u_ldown,:),lon,lat)
+         lat=lat_e(ij+u_ldown)
          IF (abs(sin(lat))<1.e-4) lat=1e-4
          x=cos(lat) ; x=x*x ; x=x*x ; x=x*x
          fact(ij+u_ldown)=(1.-x)/(2.*omega*sin(lat))
-
-!         fact(ij+u_right)=-1
-!         fact(ij+u_lup)=-1
-!         fact(ij+u_ldown)=-1
-
-!         CALL xyz2lonlat(xyz_e(ij+u_ldown,:),lon,lat)
-!         ps(ij)=lat
        ENDDO
      ENDDO
Index: codes/icosagcm/trunk/src/etat0_dcmip1.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0_dcmip1.f90	(revision 285)
+++ codes/icosagcm/trunk/src/etat0_dcmip1.f90	(revision 286)
@@ -188,6 +188,5 @@
           DO i=ii_begin-1,ii_end+1
             n=(j-1)*iim+i
-            CALL xyz2lonlat(xyz_i(n,:),lon,lat)
-            CALL dist_lonlat(lon0,lat0,lon,lat,rr1)  ! GC distance from center 
+            CALL dist_lonlat(lon0,lat0,lon_i(n),lat_i(n),rr1)  ! GC distance from center 
             rr1 = radius*rr1
             IF ( rr1 .LT. R0 ) then 
@@ -212,8 +211,7 @@
           DO i=ii_begin-1,ii_end+1
             n=(j-1)*iim+i
-            CALL xyz2lonlat(xyz_i(n,:),lon,lat)
-            CALL dist_lonlat(lonc1,latc1,lon,lat,rr1)  ! GC distance from center 
+            CALL dist_lonlat(lonc1,latc1,lon_i(n),lat_i(n),rr1)  ! GC distance from center 
             rr1 = radius*rr1
-            CALL dist_lonlat(lonc2,latc2,lon,lat,rr2)  ! GC distance from center 
+            CALL dist_lonlat(lonc2,latc2,lon_i(n),lat_i(n),rr2)  ! GC distance from center 
             rr2 = radius*rr2
             dd1t1 = rr1/rt
@@ -244,8 +242,7 @@
           DO i=ii_begin-1,ii_end+1
             n=(j-1)*iim+i
-            CALL xyz2lonlat(xyz_i(n,:),lon,lat)
-            CALL dist_lonlat(lonc1,latc1,lon,lat,rr1)  ! GC distance from center 
+            CALL dist_lonlat(lonc1,latc1,lon_i(n),lat_i(n),rr1)  ! GC distance from center 
             rr1 = radius*rr1
-            CALL dist_lonlat(lonc2,latc2,lon,lat,rr2)  ! GC distance from center 
+            CALL dist_lonlat(lonc2,latc2,lon_i(n),lat_i(n),rr2)  ! GC distance from center 
             rr2 = radius*rr2
             dd1t1 = rr1/rt
Index: codes/icosagcm/trunk/src/etat0_dcmip2.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0_dcmip2.f90	(revision 285)
+++ codes/icosagcm/trunk/src/etat0_dcmip2.f90	(revision 286)
@@ -86,5 +86,5 @@
           DO i=ii_begin,ii_end
              ij=(j-1)*iim+i
-             CALL comp_all(xyz_i(ij,:), ps(ij),phis(ij),temp(ij,l), dummy1,dummy2)
+             CALL comp_all(lon_i(ij),lat_i(ij), ps(ij),phis(ij),temp(ij,l), dummy1,dummy2)
           END DO
        END DO
@@ -100,7 +100,10 @@
           DO i=ii_begin,ii_end
              ij=(j-1)*iim+i
-             CALL comp_all(xyz_e(ij+u_right,:), dummy1,dummy2,dummy3, ulon(ij+u_right,l),ulat(ij+u_right,l))
-             CALL comp_all(xyz_e(ij+u_lup,:), dummy1,dummy2,dummy3, ulon(ij+u_lup,l),ulat(ij+u_lup,l))
-             CALL comp_all(xyz_e(ij+u_ldown,:), dummy1,dummy2,dummy3, ulon(ij+u_ldown,l),ulat(ij+u_ldown,l))
+             CALL comp_all(lon_e(ij+u_right), lat_e(ij+u_right), &
+                  dummy1,dummy2,dummy3, ulon(ij+u_right,l),ulat(ij+u_right,l))
+             CALL comp_all(lon_e(ij+u_lup), lat_e(ij+u_lup), &
+                  dummy1,dummy2,dummy3, ulon(ij+u_lup,l),ulat(ij+u_lup,l))
+             CALL comp_all(lon_e(ij+u_ldown), lat_e(ij+u_ldown), &
+                  dummy1,dummy2,dummy3, ulon(ij+u_ldown,l),ulat(ij+u_ldown,l))
           END DO
        END DO
@@ -112,11 +115,10 @@
     CONTAINS
       
-      SUBROUTINE comp_all(xyz, psj,phisj,tempj, ulonj,ulatj)
+      SUBROUTINE comp_all(lon,lat, psj,phisj,tempj, ulonj,ulatj)
         USE dcmip_initial_conditions_test_1_2_3
-        REAL(rstd), INTENT(IN) :: xyz(3)
+        REAL(rstd), INTENT(IN) :: lon, lat
         REAL(rstd), INTENT(OUT) :: psj,phisj,tempj,ulonj,ulatj
-        REAL :: lon,lat, dummy
+        REAL :: dummy
         dummy=0.
-        CALL xyz2lonlat(xyz,lon,lat)
         SELECT CASE (icase)
         CASE(0)
Index: codes/icosagcm/trunk/src/etat0_dcmip3.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0_dcmip3.f90	(revision 285)
+++ codes/icosagcm/trunk/src/etat0_dcmip3.f90	(revision 286)
@@ -83,5 +83,5 @@
   REAL(rstd) :: Rd        ! gas constant of dry air, P=rho.Rd.T
   REAL(rstd) :: alpha, rm
-  REAL(rstd) :: lon,lat, C0, C1, GG
+  REAL(rstd) :: C0, C1, GG
   REAL(rstd) :: p0psk, pspsk,r,zz, thetab, thetap
   REAL(rstd) :: dummy, pp
@@ -103,10 +103,9 @@
      DO i=ii_begin,ii_end
         ij=(j-1)*iim+i
-        CALL xyz2lonlat(xyz_i(ij,:),lon,lat)
 
         IF(use_dcmip_routine) THEN
-           CALL test3_gravity_wave(lon,lat,pp,dummy,0, dummy,dummy,dummy,dummy,phis(ij),ps(ij),dummy,dummy)
+           CALL test3_gravity_wave(lon_i(ij),lat_i(ij),pp,dummy,0, dummy,dummy,dummy,dummy,phis(ij),ps(ij),dummy,dummy)
         ELSE
-           C1=C0*(cos(2*lat)-1)
+           C1=C0*(cos(2*lat_i(ij))-1)
            
            !--- GROUND GEOPOTENTIAL
@@ -120,5 +119,5 @@
            
            
-           r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc))
+           r=radius*acos(sin(latc)*sin(lat_i(ij))+cos(latc)*cos(lat_i(ij))*cos(lon_i(ij)-lonc))
            s(ij)= d**2/(d**2+r**2)
         END IF
@@ -134,6 +133,6 @@
            pp=0.5*(p(ij,l+1)+p(ij,l))  ! full-layer pressure
            IF(use_dcmip_routine) THEN
-              CALL xyz2lonlat(xyz_i(ij,:),lon,lat)
-              CALL test3_gravity_wave(lon,lat,pp,dummy,0,dummy,dummy,dummy,T(ij,l),dummy,dummy,dummy,dummy)
+              CALL test3_gravity_wave(lon_i(ij),lat_i(ij),pp,dummy,0, &
+                   dummy,dummy,dummy,T(ij,l),dummy,dummy,dummy,dummy)
            ELSE
               pspsk=(pp/ps(ij))**kappa
@@ -162,28 +161,20 @@
            ij=(j-1)*iim+i
            IF(use_dcmip_routine) THEN
-              CALL xyz2lonlat(xyz_e(ij+u_right,:),lon,lat)
-              CALL test3_gravity_wave(lon,lat,pp,0.,0, &
-                   ulon(ij+u_right,l),ulat(ij+u_right,l),&
+              CALL test3_gravity_wave(lon_e(ij+u_right),lat_e(ij+u_right), &
+                   pp,0.,0, ulon(ij+u_right,l),ulat(ij+u_right,l),&
                    dummy,dummy,dummy,dummy,dummy,dummy)
-              CALL xyz2lonlat(xyz_e(ij+u_lup,:),lon,lat)
-              CALL test3_gravity_wave(lon,lat,pp,0.,0,&
-                   ulon(ij+u_lup,l),ulat(ij+u_lup,l),&
+              CALL test3_gravity_wave(lon_e(ij+u_lup),lat_e(ij+u_lup), &
+                   pp,0.,0, ulon(ij+u_lup,l),ulat(ij+u_lup,l),&
                    dummy,dummy,dummy,dummy,dummy,dummy)
-              CALL xyz2lonlat(xyz_e(ij+u_ldown,:),lon,lat)
-              CALL test3_gravity_wave(lon,lat,pp,0.,0,&
-                   ulon(ij+u_ldown,l),ulat(ij+u_ldown,l),&
+              CALL test3_gravity_wave(lon_e(ij+u_ldown),lat_e(ij+u_ldown), &
+                   pp,0.,0, ulon(ij+u_ldown,l),ulat(ij+u_ldown,l),&
                    dummy,dummy,dummy,dummy,dummy,dummy)
            ELSE
-              CALL xyz2lonlat(xyz_e(ij+u_right,:),lon,lat)
-              ulon(ij+u_right,l)=u0*cos(lat)
-              ulat(ij+u_right,l)=0
-              
-              CALL xyz2lonlat(xyz_e(ij+u_lup,:),lon,lat)
-              ulon(ij+u_lup,l)=u0*cos(lat)
-              ulat(ij+u_lup,l)=0
-              
-              CALL xyz2lonlat(xyz_e(ij+u_ldown,:),lon,lat)
-              ulon(ij+u_ldown,l)=u0*cos(lat)
-              ulat(ij+u_ldown,l)=0
+              ulon(ij+u_right,l) = u0*cos(lat_e(ij+u_right))
+              ulat(ij+u_right,l) = 0
+              ulon(ij+u_lup,l)   = u0*cos(lat_e(ij+u_lup))
+              ulat(ij+u_lup,l)   = 0
+              ulon(ij+u_ldown,l) = u0*cos(lat_e(ij+u_ldown))
+              ulat(ij+u_ldown,l) = 0
            END IF
         ENDDO
Index: codes/icosagcm/trunk/src/etat0_dcmip4.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0_dcmip4.f90	(revision 285)
+++ codes/icosagcm/trunk/src/etat0_dcmip4.f90	(revision 286)
@@ -121,5 +121,5 @@
           ij=(j-1)*iim+i
           
-          CALL xyz2lonlat(xyz_e(ij+u_right,:)/radius,lon,lat)
+          lon=lon_e(ij+u_right) ; lat=lat_e(ij+u_right)
           K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)
           r=radius*acos(K)
@@ -127,6 +127,5 @@
           u(ij+u_right,l) = utot * sum(elon_e(ij+u_right,:) * ep_e(ij+u_right,:))
 
-
-          CALL xyz2lonlat(xyz_e(ij+u_lup,:)/radius,lon,lat)
+          lon=lon_e(ij+u_lup) ; lat=lat_e(ij+u_lup)
           K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)
           r=radius*acos(K)
@@ -134,5 +133,5 @@
           u(ij+u_lup,l) = utot * sum(elon_e(ij+u_lup,:) * ep_e(ij+u_lup,:))
 
-          CALL xyz2lonlat(xyz_e(ij+u_ldown,:)/radius,lon,lat)
+          lon=lon_e(ij+u_ldown) ; lat=lat_e(ij+u_ldown)
           K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)
           r=radius*acos(K)
@@ -151,6 +150,5 @@
          DO i=ii_begin-1,ii_end+1
            ij=(j-1)*iim+i
-           CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat)
-            
+           lat=lat_i(ij)
             Y(ij,l)=((-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)
@@ -168,5 +166,5 @@
        DO i=ii_begin,ii_end
          ij=(j-1)*iim+i
-         CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat)
+         lat=lat_i(ij)
          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 )
@@ -191,5 +189,5 @@
            DO i=ii_begin,ii_end
              ij=(j-1)*iim+i
-             CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat)
+             lon=lon_i(ij) ; lat=lat_i(ij)
              dthetaodeta=dthetaodeta_ave + 3./4. * Pi * u0/Rd*(1-kappa)*eta(l)**(-kappa)*sin(etav(l))*cos(etav(l))**0.5 * Y(ij,l) & 
                                         + 3/8. * Pi**2*u0/Rd * eta(l)**(1-kappa) * cos(etav(l))**1.5 * Y(ij,l)                    & 
@@ -220,6 +218,5 @@
            DO i=ii_begin,ii_end
              ij=(j-1)*iim+i
-             CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat)
-             q(ij,l,1)=q0*exp(-(lat/latw)**4)*exp(-((eta(l)-1)*preff/pw)**2)
+             q(ij,l,1)=q0*exp(-(lat_i(ij)/latw)**4)*exp(-((eta(l)-1)*preff/pw)**2)
            ENDDO
          ENDDO
Index: codes/icosagcm/trunk/src/etat0_heldsz.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0_heldsz.f90	(revision 285)
+++ codes/icosagcm/trunk/src/etat0_heldsz.f90	(revision 286)
@@ -6,5 +6,4 @@
   TYPE(t_field),POINTER :: f_theta_eq(:)
   TYPE(t_field),POINTER :: f_theta(:)
-  TYPE(t_field),POINTER :: f_clat(:) ! FIXME, duplication
 
   REAL(rstd),ALLOCATABLE,SAVE :: knewt_t(:),kfrict(:)
@@ -65,5 +64,4 @@
     REAL(rstd),POINTER :: u(:,:)
     REAL(rstd),POINTER :: q(:,:,:)
-    REAL(rstd),POINTER :: clat(:) 
     REAL(rstd),POINTER :: theta_eq(:,:) 
     REAL(rstd),POINTER :: theta(:,:) 
@@ -78,6 +76,5 @@
 
        theta_eq=f_theta_eq(ind) 
-       clat=f_clat(ind) 
-       CALL compute_Teq(clat,theta_eq) ! FIXME : already done by Init_Teq
+       CALL compute_Teq(lat_i,theta_eq) ! FIXME : already done by Init_Teq
 
        ps=f_ps(ind)
@@ -111,5 +108,4 @@
        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)) 
 
@@ -147,7 +143,6 @@
           CALL swap_dimensions(ind)
           CALL swap_geometry(ind)
-          clat=f_clat(ind)
           theta_eq=f_theta_eq(ind)
-          CALL compute_Teq(clat,theta_eq)
+          CALL compute_Teq(lat_i,theta_eq)
        ENDDO
 
@@ -159,23 +154,13 @@
   END SUBROUTINE init_Teq
 
-  SUBROUTINE compute_Teq(clat,theta_eq)
+  SUBROUTINE compute_Teq(lat,theta_eq)
     USE icosa
     USE disvert_mod
     IMPLICIT NONE  
-    REAL(rstd),INTENT(OUT) :: clat(iim*jjm)
+    REAL(rstd),INTENT(IN) :: lat(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) 
+    REAL(rstd) :: r, zsig, ddsin, tetastrat, tetajl
     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
@@ -185,5 +170,5 @@
           DO i=ii_begin-1,ii_end+1
              ij=(j-1)*iim+i
-	     ddsin=slat(ij) 
+	     ddsin=SIN(lat(ij)) 
 	     tetajl=teta0-delt_y*ddsin*ddsin+eps*ddsin &
                   -delt_z*(1.-ddsin*ddsin)*log(zsig)
@@ -244,10 +229,9 @@
        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) 
+       CALL compute_heldsz(ps,theta_eq,lat_i, theta_rhodz,u, theta) 
     ENDDO
   END SUBROUTINE held_suarez
 
-  SUBROUTINE compute_heldsz(ps,theta_eq,clat, theta_rhodz,u, theta) 
+  SUBROUTINE compute_heldsz(ps,theta_eq,lat, theta_rhodz,u, theta) 
     USE icosa
     USE theta2theta_rhodz_mod
@@ -255,5 +239,5 @@
     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(IN)    :: lat(iim*jjm) 
     REAL(rstd),INTENT(INOUT) :: theta_rhodz(iim*jjm,llm)
     REAL(rstd),INTENT(INOUT) :: u(3*iim*jjm,llm)
@@ -268,5 +252,5 @@
              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 )
+                  (knewt_g+knewt_t(l)*COS(lat(ij))**4 )
           ENDDO
        ENDDO
Index: codes/icosagcm/trunk/src/geometry.f90
===================================================================
--- codes/icosagcm/trunk/src/geometry.f90	(revision 285)
+++ codes/icosagcm/trunk/src/geometry.f90	(revision 286)
@@ -3,8 +3,12 @@
   
   TYPE t_geometry
+    TYPE(t_field),POINTER :: centroid(:)
     TYPE(t_field),POINTER :: xyz_i(:)
-    TYPE(t_field),POINTER :: centroid(:)
     TYPE(t_field),POINTER :: xyz_e(:)
     TYPE(t_field),POINTER :: xyz_v(:)
+    TYPE(t_field),POINTER :: lon_i(:)
+    TYPE(t_field),POINTER :: lon_e(:)
+    TYPE(t_field),POINTER :: lat_i(:)
+    TYPE(t_field),POINTER :: lat_e(:)
     TYPE(t_field),POINTER :: ep_e(:)
     TYPE(t_field),POINTER :: et_e(:)
@@ -30,10 +34,20 @@
   REAL(rstd),POINTER :: Ai(:)          ! area of a cell
 !$OMP THREADPRIVATE(Ai)
+  REAL(rstd),POINTER :: centroid(:,:)  ! coordinate of the centroid of the cell
+!$OMP THREADPRIVATE(centroid)
   REAL(rstd),POINTER :: xyz_i(:,:)     ! coordinate of the center of the cell (voronoi)
 !$OMP THREADPRIVATE(xyz_i)
-  REAL(rstd),POINTER :: centroid(:,:)  ! coordinate of the centroid of the cell
-!$OMP THREADPRIVATE(centroid)
   REAL(rstd),POINTER :: xyz_e(:,:)     ! coordinate of a wind point on the cell on a edge
 !$OMP THREADPRIVATE(xyz_e)
+  REAL(rstd),POINTER :: xyz_v(:,:)     ! coordinate of a vertex (center of the dual mesh)
+!$OMP THREADPRIVATE(xyz_v)
+  REAL(rstd),POINTER :: lon_i(:)       ! longitude of the center of the cell (voronoi)
+!$OMP THREADPRIVATE(lon_i)
+  REAL(rstd),POINTER :: lon_e(:)       ! longitude of a wind point on the cell on a edge
+!$OMP THREADPRIVATE(lon_e)
+  REAL(rstd),POINTER :: lat_i(:)       ! latitude of the center of the cell (voronoi)
+!$OMP THREADPRIVATE(lat_i)
+  REAL(rstd),POINTER :: lat_e(:)       ! latitude of a wind point on the cell on a edge
+!$OMP THREADPRIVATE(lat_e)
   REAL(rstd),POINTER :: ep_e(:,:)      ! perpendicular unit vector of a edge (outsider)
 !$OMP THREADPRIVATE(ep_e)
@@ -48,6 +62,4 @@
   REAL(rstd),POINTER :: elat_e(:,:)    ! unit latitude vector on a wind point
 !$OMP THREADPRIVATE(elat_e)
-  REAL(rstd),POINTER :: xyz_v(:,:)     ! coordinate of a vertex (center of the dual mesh)
-!$OMP THREADPRIVATE(xyz_v)
   REAL(rstd),POINTER :: Av(:)          ! area of dual mesk cell
 !$OMP THREADPRIVATE(Av)
@@ -84,13 +96,20 @@
   
     CALL allocate_field(geom%Ai,field_t,type_real,name='Ai')
+
     CALL allocate_field(geom%xyz_i,field_t,type_real,3)
+    CALL allocate_field(geom%lon_i,field_t,type_real)
+    CALL allocate_field(geom%lat_i,field_t,type_real)
+    CALL allocate_field(geom%elon_i,field_t,type_real,3)
+    CALL allocate_field(geom%elat_i,field_t,type_real,3)
     CALL allocate_field(geom%centroid,field_t,type_real,3)
+
     CALL allocate_field(geom%xyz_e,field_u,type_real,3)
+    CALL allocate_field(geom%lon_e,field_u,type_real)
+    CALL allocate_field(geom%lat_e,field_u,type_real)
+    CALL allocate_field(geom%elon_e,field_u,type_real,3)
+    CALL allocate_field(geom%elat_e,field_u,type_real,3)
     CALL allocate_field(geom%ep_e,field_u,type_real,3)
     CALL allocate_field(geom%et_e,field_u,type_real,3)
-    CALL allocate_field(geom%elon_i,field_t,type_real,3)
-    CALL allocate_field(geom%elat_i,field_t,type_real,3)
-    CALL allocate_field(geom%elon_e,field_u,type_real,3)
-    CALL allocate_field(geom%elat_e,field_u,type_real,3)
+
     CALL allocate_field(geom%xyz_v,field_z,type_real,3)
     CALL allocate_field(geom%de,field_u,type_real)
@@ -119,4 +138,8 @@
     ep_e=geom%ep_e(ind)
     et_e=geom%et_e(ind)
+    lon_i=geom%lon_i(ind)
+    lat_i=geom%lat_i(ind)
+    lon_e=geom%lon_e(ind)
+    lat_e=geom%lat_e(ind)
     elon_i=geom%elon_i(ind)
     elat_i=geom%elat_i(ind)
@@ -394,4 +417,6 @@
 
           CALL xyz2lonlat(xyz_i(n,:),lon,lat)
+          lon_i(n)=lon
+          lat_i(n)=lat
           elon_i(n,1) = -sin(lon)
           elon_i(n,2) = cos(lon)
@@ -403,4 +428,6 @@
           
           CALL xyz2lonlat(xyz_e(n+u_right,:),lon,lat)
+          lon_e(n+u_right)=lon
+          lat_e(n+u_right)=lat
           elon_e(n+u_right,1) = -sin(lon)
           elon_e(n+u_right,2) = cos(lon)
@@ -411,4 +438,6 @@
           
           CALL xyz2lonlat(xyz_e(n+u_lup,:),lon,lat)
+          lon_e(n+u_lup)=lon
+          lat_e(n+u_lup)=lat
           elon_e(n+u_lup,1) = -sin(lon)
           elon_e(n+u_lup,2) = cos(lon)
@@ -419,4 +448,6 @@
  
           CALL xyz2lonlat(xyz_e(n+u_ldown,:),lon,lat)
+          lon_e(n+u_ldown)=lon
+          lat_e(n+u_ldown)=lat
           elon_e(n+u_ldown,1) = -sin(lon)
           elon_e(n+u_ldown,2) = cos(lon)
Index: codes/icosagcm/trunk/src/physics_interface.f90
===================================================================
--- codes/icosagcm/trunk/src/physics_interface.f90	(revision 285)
+++ codes/icosagcm/trunk/src/physics_interface.f90	(revision 286)
@@ -194,20 +194,8 @@
        CALL swap_geometry(ind)
        CALL pack_domain_2D(pack_info(ind), Ai, physics_inout%Ai)
-       CALL pack_lonlat(pack_info(ind))
+       CALL pack_domain_2D(pack_info(ind), lon_i, physics_inout%lon)
+       CALL pack_domain_2D(pack_info(ind), lat_i, physics_inout%lat)
     END DO
   END SUBROUTINE init_pack_after
-
-  SUBROUTINE pack_lonlat(info)
-    USE icosa
-    IMPLICIT NONE
-    TYPE(t_pack_info) :: info
-    REAL(rstd) :: lon(iim*jjm), lat(iim*jjm)
-    INTEGER :: ij
-    DO ij=1,iim*jjm
-       CALL xyz2lonlat(xyz_i(ij,:),lon(ij),lat(ij)) 
-    ENDDO
-    CALL pack_domain_2D(info, lon, physics_inout%lon)
-    CALL pack_domain_2D(info, lat, physics_inout%lat)
-  END SUBROUTINE pack_lonlat
 
 !-------------------------------- Pack / Unpack 2D ---------------------------
