Index: codes/icosagcm/trunk/src/etat0_temperature.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0_temperature.f90	(revision 328)
+++ codes/icosagcm/trunk/src/etat0_temperature.f90	(revision 328)
@@ -0,0 +1,81 @@
+MODULE etat0_temperature_mod
+  USE prec
+  USE icosa, ONLY: llm,nqtot
+  IMPLICIT NONE
+  REAL(rstd), SAVE, ALLOCATABLE :: t_profile(:)
+!$OMP THREADPRIVATE(t_profile)
+CONTAINS
+  
+  SUBROUTINE getin_etat0
+    USE getin_mod, ONLY: getin
+    USE mpipara, ONLY: is_mpi_root
+    USE omp_para, ONLY: omp_in_parallel
+    USE transfert_omp_mod, ONLY: bcast_omp
+    USE free_unit_mod, ONLY : free_unit
+    IMPLICIT NONE
+    INTEGER :: unit,ok
+    INTEGER :: l
+    CHARACTER(len=255) :: temperature_file
+  
+    temperature_file="profile.in" ! default file name
+    ! but users may want to use some other file name
+    CALL getin("temperature_profile_file",temperature_file)
+    
+    ALLOCATE(t_profile(llm))
+ 
+ 
+ !$OMP MASTER
+    unit=free_unit()
+    OPEN(unit,file=temperature_file,status="old",action="read",iostat=ok)
+    IF (ok/=0) THEN
+      WRITE(*,*) "getin_etat0 error: input file ",trim(temperature_file)," not found!"
+      CALL ABORT
+    ENDIF
+    ! read in t_profile() line by line, starting from first atmospheric
+    ! layer, up to model top
+    DO l=1,llm
+      READ(unit,fmt=*,iostat=ok) t_profile(l)
+      IF (ok/=0) THEN
+        WRITE(*,*) "getin_etat0 error: failed reading t_profile(l) for l=",l
+        CALL ABORT
+      ENDIF
+    ENDDO
+    
+    CLOSE(unit)
+    IF (is_mpi_root) THEN
+      WRITE(*,*) "Using input temperature profile from file ",trim(temperature_file),":"
+      DO l=1,llm
+         WRITE(*,*) "  TEMP(l=",l,")=",t_profile(l)
+      ENDDO
+    ENDIF
+!$OMP END MASTER
+
+    IF (omp_in_parallel()) THEN
+      CALL bcast_omp(t_profile)
+    ENDIF
+    
+  END SUBROUTINE getin_etat0
+
+  SUBROUTINE compute_etat0(ngrid, phis, ps, temp, ulon, ulat, q)
+    USE earth_const, ONLY: preff
+    IMPLICIT NONE  
+    INTEGER, INTENT(IN)    :: 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
+    phis(:)=0
+    ps(:)=preff
+    DO l=1,llm
+      temp(:,l)=t_profile(l)
+    ENDDO
+    ulon(:,:)=0
+    ulat(:,:)=0
+    q(:,:,:)=0
+
+  END SUBROUTINE compute_etat0
+
+END MODULE etat0_temperature_mod
Index: codes/icosagcm/trunk/src/physics_lmdz_generic.F90
===================================================================
--- codes/icosagcm/trunk/src/physics_lmdz_generic.F90	(revision 328)
+++ codes/icosagcm/trunk/src/physics_lmdz_generic.F90	(revision 328)
@@ -0,0 +1,613 @@
+MODULE physics_lmdz_generic_mod
+  USE field_mod
+  USE transfert_mod
+  
+  INTEGER,SAVE :: nbp_phys
+  TYPE(t_message) :: req_u
+
+  TYPE(t_field),POINTER :: f_p(:) 
+  TYPE(t_field),POINTER :: f_pks(:)  
+  TYPE(t_field),POINTER :: f_pk(:)  
+  TYPE(t_field),POINTER :: f_p_layer(:)    
+  TYPE(t_field),POINTER :: f_theta(:)   
+  TYPE(t_field),POINTER :: f_phi(:)   
+  TYPE(t_field),POINTER :: f_Temp(:)    
+  TYPE(t_field),POINTER :: f_ulon(:)    
+  TYPE(t_field),POINTER :: f_ulat(:)   
+  TYPE(t_field),POINTER :: f_dulon(:)
+  TYPE(t_field),POINTER :: f_dulat(:)
+  TYPE(t_field),POINTER :: f_dTemp(:)
+  TYPE(t_field),POINTER :: f_dq(:)
+  TYPE(t_field),POINTER :: f_dps(:)
+  TYPE(t_field),POINTER :: f_duc(:)
+
+  INTEGER :: start_clock
+  INTEGER :: stop_clock
+  INTEGER :: count_clock=0
+  
+  REAL :: start_day
+  REAL :: day_length
+  REAL :: year_length
+  
+ INTEGER,ALLOCATABLE,SAVE :: domain_offset(:)
+
+  
+
+CONTAINS
+
+  SUBROUTINE init_physics
+  USE icosa
+  USE domain_mod
+  USE dimensions
+  USE mpi_mod
+  USE mpipara
+  USE disvert_mod
+  USE xios_mod
+  
+  IMPLICIT NONE
+  INTEGER  :: distrib(0:mpi_size-1)
+  INTEGER  :: ind,i,j,ij,pos
+
+  REAL(rstd),ALLOCATABLE :: latfi(:)
+  REAL(rstd),ALLOCATABLE :: lonfi(:)
+  REAL(rstd),ALLOCATABLE :: airefi(:)
+  REAL(rstd),ALLOCATABLE :: bounds_latfi(:,:)
+  REAL(rstd),ALLOCATABLE :: bounds_lonfi(:,:)
+
+    start_day=0
+    day_length=86400
+    year_length=86400*365.25
+    
+    CALL getin('start_day',start_day)
+    CALL getin('day_length',day_length)
+    CALL getin('year_length',year_length)
+
+!$OMP PARALLEL
+    CALL allocate_field(f_p,field_t,type_real,llm+1)
+    CALL allocate_field(f_pks,field_t,type_real)
+    CALL allocate_field(f_pk,field_t,type_real,llm)
+    CALL allocate_field(f_p_layer,field_t,type_real,llm)
+    CALL allocate_field(f_theta,field_t,type_real,llm)
+    CALL allocate_field(f_phi,field_t,type_real,llm)
+    CALL allocate_field(f_Temp,field_t,type_real,llm)
+    CALL allocate_field(f_ulon,field_t,type_real,llm)
+    CALL allocate_field(f_ulat,field_t,type_real,llm)
+    CALL allocate_field(f_dulon,field_t,type_real,llm)
+    CALL allocate_field(f_dulat,field_t,type_real,llm)
+    CALL allocate_field(f_dTemp,field_t,type_real,llm)
+    CALL allocate_field(f_dq,field_t,type_real,llm,nqtot)
+    CALL allocate_field(f_dps,field_t,type_real)
+    CALL allocate_field(f_duc,field_t,type_real,3,llm)    
+!$OMP END PARALLEL    
+
+    ALLOCATE(domain_offset(ndomain))
+    nbp_phys=0
+    domain_offset(1)=0
+    DO ind=1,ndomain
+      CALL swap_dimensions(ind)
+      IF (ind<ndomain) THEN 
+        domain_offset(ind+1)=domain_offset(ind)+ii_nb*jj_nb
+      ENDIF
+      nbp_phys=nbp_phys+ii_nb*jj_nb
+    ENDDO
+    
+    CALL MPI_ALLGATHER(nbp_phys,1,MPI_INTEGER,distrib,1,MPI_INTEGER,comm_icosa,ierr)
+    
+    ALLOCATE(latfi(nbp_phys))
+    ALLOCATE(lonfi(nbp_phys))
+    ALLOCATE(airefi(nbp_phys))
+    ALLOCATE(bounds_latfi(nbp_phys,6))
+    ALLOCATE(bounds_lonfi(nbp_phys,6))
+    
+    pos=0
+    DO ind=1,ndomain
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          pos=pos+1
+          CALL xyz2lonlat(xyz_i(ij,:),lonfi(pos),latfi(pos))
+          CALL xyz2lonlat(xyz_v(ij+z_rup,:), bounds_lonfi(pos,1), bounds_latfi(pos,1))
+          CALL xyz2lonlat(xyz_v(ij+z_up,:), bounds_lonfi(pos,2), bounds_latfi(pos,2))
+          CALL xyz2lonlat(xyz_v(ij+z_lup,:), bounds_lonfi(pos,3), bounds_latfi(pos,3))
+          CALL xyz2lonlat(xyz_v(ij+z_ldown,:), bounds_lonfi(pos,4), bounds_latfi(pos,4))
+          CALL xyz2lonlat(xyz_v(ij+z_down,:), bounds_lonfi(pos,5), bounds_latfi(pos,5))
+          CALL xyz2lonlat(xyz_v(ij+z_rdown,:), bounds_lonfi(pos,6), bounds_latfi(pos,6))
+          airefi(pos)=Ai(ij)  
+        ENDDO
+      ENDDO
+    ENDDO
+
+
+#ifdef CPP_PHYSICS_LMDZ_GENERIC
+    CALL initialize_unstructured_physics(nbp_phys,llm, comm_icosa, mpi_size,distrib,            &
+                                         day_length,start_day,itau_physics*dt,                  &
+                                         6,latfi,lonfi,airefi,bounds_lonfi,bounds_latfi,         &
+                                         radius,g, gas_constant/mu, cpp,                        &
+                                         preff, ap, bp)                                         
+#endif
+    
+    CALL xios_set_context    
+    
+!    CALL init_phys_lmdz(128,97,llm, comm_icosa, mpi_size, distrib)
+    
+#ifdef CPP_PHYSICS_LMDZ_GENERIC
+    CALL init_gcm_lmdz(nbp_phys,mpi_size,distrib,latfi,lonfi,airefi)
+#endif
+!    CALL init_phys_lmdz(128,97,llm, comm_icosa, mpi_size, distrib)
+!    CALL iniphysiq(llm,day_length,start_day,itau_physics*dt,     &
+!                   latfi,lonfi,airefi,radius,g, gas_constant/mu, cpp,     &
+!                   ap, bp)
+    
+    
+    
+  END SUBROUTINE init_physics
+  
+  SUBROUTINE physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_wflux, f_q)
+  USE ICOSA
+  USE time_mod
+  USE disvert_mod
+  USE transfert_mod
+  USE mpipara
+  USE xios_mod
+  USE trace
+  IMPLICIT NONE
+    INTEGER,INTENT(IN)    :: it
+    TYPE(t_field),POINTER :: f_phis(:)
+    TYPE(t_field),POINTER :: f_ps(:)
+    TYPE(t_field),POINTER :: f_theta_rhodz(:)
+    TYPE(t_field),POINTER :: f_u(:)
+    TYPE(t_field),POINTER :: f_wflux(:)
+    TYPE(t_field),POINTER :: f_q(:)
+  
+    REAL(rstd),POINTER :: phis(:)
+    REAL(rstd),POINTER :: ps(:)
+    REAL(rstd),POINTER :: theta_rhodz(:,:)
+    REAL(rstd),POINTER :: u(:,:)
+    REAL(rstd),POINTER :: wflux(:,:)
+    REAL(rstd),POINTER :: q(:,:,:)
+    REAL(rstd),POINTER :: p(:,:)
+    REAL(rstd),POINTER :: pks(:)
+    REAL(rstd),POINTER :: pk(:,:)
+    REAL(rstd),POINTER :: p_layer(:,:)
+    REAL(rstd),POINTER :: theta(:,:)
+    REAL(rstd),POINTER :: phi(:,:)
+    REAL(rstd),POINTER :: Temp(:,:)
+    REAL(rstd),POINTER :: ulon(:,:)
+    REAL(rstd),POINTER :: ulat(:,:)
+    REAL(rstd),POINTER :: dulon(:,:)
+    REAL(rstd),POINTER :: dulat(:,:)
+    REAL(rstd),POINTER :: dTemp(:,:)
+    REAL(rstd),POINTER :: dq(:,:,:)
+    REAL(rstd),POINTER :: dps(:)
+    REAL(rstd),POINTER :: duc(:,:,:)
+
+
+    INTEGER :: ind
+    
+    REAL(rstd),ALLOCATABLE,SAVE :: ps_phy(:)
+    REAL(rstd),ALLOCATABLE,SAVE :: p_phy(:,:)
+    REAL(rstd),ALLOCATABLE,SAVE :: p_layer_phy(:,:)
+    REAL(rstd),ALLOCATABLE,SAVE :: Temp_phy(:,:)
+    REAL(rstd),ALLOCATABLE,SAVE :: phis_phy(:)
+    REAL(rstd),ALLOCATABLE,SAVE :: phi_phy(:,:)
+    REAL(rstd),ALLOCATABLE,SAVE :: ulon_phy(:,:)
+    REAL(rstd),ALLOCATABLE,SAVE :: ulat_phy(:,:)
+    REAL(rstd),ALLOCATABLE,SAVE :: q_phy(:,:,:)
+    REAL(rstd),ALLOCATABLE,SAVE :: wflux_phy(:,:)
+    REAL(rstd),ALLOCATABLE,SAVE :: dulon_phy(:,:)
+    REAL(rstd),ALLOCATABLE,SAVE :: dulat_phy(:,:)
+    REAL(rstd),ALLOCATABLE,SAVE :: dTemp_phy(:,:)
+    REAL(rstd),ALLOCATABLE,SAVE :: dq_phy(:,:,:)
+    REAL(rstd),ALLOCATABLE,SAVE :: dps_phy(:)
+    REAL(rstd) :: dtphy 
+    REAL(rstd) :: time
+    REAL(rstd) :: day
+    REAL(rstd) :: real_time
+    INTEGER    :: offset
+    LOGICAL :: lafin=.FALSE.
+    LOGICAL,SAVE :: first=.TRUE.
+!$OMP THREADPRIVATE(first)
+
+    
+    IF (first) THEN
+      first=.FALSE.
+      CALL init_message(f_u,req_e1_vect,req_u)
+!$OMP MASTER
+      ALLOCATE(ps_phy(nbp_phys))
+      ALLOCATE(p_phy(nbp_phys,llm+1))
+      ALLOCATE(p_layer_phy(nbp_phys,llm))
+      ALLOCATE(Temp_phy(nbp_phys,llm))
+      ALLOCATE(phis_phy(nbp_phys))
+      ALLOCATE(phi_phy(nbp_phys,llm))
+      ALLOCATE(ulon_phy(nbp_phys,llm))
+      ALLOCATE(ulat_phy(nbp_phys,llm))
+      ALLOCATE(q_phy(nbp_phys,llm,nqtot))
+      ALLOCATE(wflux_phy(nbp_phys,llm))
+      ALLOCATE(dulon_phy(nbp_phys,llm))
+      ALLOCATE(dulat_phy(nbp_phys,llm))
+      ALLOCATE(dTemp_phy(nbp_phys,llm))
+      ALLOCATE(dq_phy(nbp_phys,llm,nqtot))
+      ALLOCATE(dps_phy(nbp_phys))
+      
+!$OMP END MASTER
+!$OMP BARRIER
+    ENDIF
+
+    IF(it==itaumax) THEN 
+      lafin=.TRUE.
+    ELSE 
+      lafin=.FALSE.
+    ENDIF
+
+!$OMP MASTER        
+!    CALL update_calendar(it)
+!$OMP END MASTER
+!$OMP BARRIER
+    dtphy=itau_physics*dt
+    real_time=start_day*day_length+it*dt
+    day  = INT( MODULO(real_time,year_length) / day_length) 
+    time = MODULO(real_time,day_length) / day_length
+
+!$OMP MASTER    
+    IF (is_mpi_root) PRINT*,"Enterring in physic : day", day, "  time : ",time 
+!$OMP END MASTER    
+    
+    
+    
+    
+    CALL transfert_message(f_u,req_u)
+    
+    DO ind=1,ndomain
+      CALL swap_dimensions(ind)
+      IF (assigned_domain(ind)) THEN
+        offset=domain_offset(ind)
+        CALL swap_geometry(ind)
+      
+        phis=f_phis(ind)
+        ps=f_ps(ind)
+        theta_rhodz=f_theta_rhodz(ind)
+        u=f_u(ind)
+        q=f_q(ind)
+        wflux=f_wflux(ind)
+        p=f_p(ind)
+        pks=f_pks(ind)
+        pk=f_pk(ind)
+        p_layer=f_p_layer(ind)
+        theta=f_theta(ind)
+        phi=f_phi(ind)
+        Temp=f_Temp(ind)
+        ulon=f_ulon(ind)
+        ulat=f_ulat(ind)
+            
+        CALL grid_icosa_to_physics
+
+      ENDIF
+    ENDDO
+
+!$OMP BARRIER
+!$OMP MASTER
+    CALL SYSTEM_CLOCK(start_clock)
+!$OMP END MASTER
+    CALL trace_start("physic")
+!    CALL trace_off()
+
+#ifdef CPP_PHYSICS_LMDZ_GENERIC
+    CALL calfis_icosa(dtphy, lafin, day, time, presnivs,      &
+                      p_phy, p_layer_phy, phi_phy, phis_phy, ulon_phy, ulat_phy, Temp_phy, q_phy, wflux_phy, &
+                      dulon_phy, dulat_phy, dTemp_phy, dq_phy, dps_phy  )
+#endif
+!    CALL trace_on()
+    CALL trace_end("physic")
+!$OMP MASTER
+    CALL SYSTEM_CLOCK(stop_clock)
+    count_clock=count_clock+stop_clock-start_clock
+!$OMP END MASTER
+
+!$OMP BARRIER                       
+
+    DO ind=1,ndomain
+      CALL swap_dimensions(ind)
+      IF (assigned_domain(ind)) THEN
+        CALL swap_geometry(ind)
+        offset=domain_offset(ind)
+      
+        theta_rhodz=f_theta_rhodz(ind)
+        u=f_u(ind)
+        q=f_q(ind)
+        ps=f_ps(ind)
+        dulon=f_dulon(ind)
+        dulat=f_dulat(ind)
+        Temp=f_temp(ind)
+        dTemp=f_dTemp(ind)
+        dq=f_dq(ind)
+        dps=f_dps(ind)
+        duc=f_duc(ind)
+        p=f_p(ind)
+        pks=f_pks(ind)
+        pk=f_pk(ind)
+      
+        CALL grid_physics_to_icosa
+      ENDIF
+    ENDDO
+
+!$OMP BARRIER
+    CALL xios_set_context    
+     
+  CONTAINS
+   
+    SUBROUTINE grid_physics_to_icosa
+    USE theta2theta_rhodz_mod
+    USE omp_para
+    IMPLICIT NONE
+      INTEGER :: i,j,ij,l,iq
+          
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          offset=offset+1
+          
+          dulon(ij,ll_begin:ll_end)=dulon_phy(offset,ll_begin:ll_end)
+          dulat(ij,ll_begin:ll_end)=dulat_phy(offset,ll_begin:ll_end)
+          dTemp(ij,ll_begin:ll_end)=dTemp_phy(offset,ll_begin:ll_end)
+          Temp(ij,ll_begin:ll_end) = Temp_phy(offset,ll_begin:ll_end)
+          dq(ij,ll_begin:ll_end,:)=dq_phy(offset,ll_begin:ll_end,:)
+          if (is_omp_first_level) dps(ij)=dps_phy(offset)
+        ENDDO
+      ENDDO
+    
+      DO l=ll_begin,ll_end
+        DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+            ij=(j-1)*iim+i
+            duc(ij,:,l)=dulon(ij,l)*elon_i(ij,:)+dulat(ij,l)*elat_i(ij,:)
+          ENDDO
+        ENDDO
+      ENDDO
+
+      DO l=ll_begin,ll_end
+        DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+            ij=(j-1)*iim+i
+            u(ij+u_right,l) = u(ij+u_right,l) + dtphy * sum( 0.5*(duc(ij,:,l) + duc(ij+t_right,:,l))*ep_e(ij+u_right,:) )
+            u(ij+u_lup,l) = u(ij+u_lup,l) + dtphy * sum( 0.5*(duc(ij,:,l) + duc(ij+t_lup,:,l))*ep_e(ij+u_lup,:) )
+            u(ij+u_ldown,l) = u(ij+u_ldown,l) + dtphy*sum( 0.5*(duc(ij,:,l) + duc(ij+t_ldown,:,l))*ep_e(ij+u_ldown,:) )
+          ENDDO
+        ENDDO
+      ENDDO          
+
+      DO l=ll_begin,ll_end
+        DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+            ij=(j-1)*iim+i
+            Temp(ij,l)=Temp(ij,l)+ dtphy * dTemp(ij,l)
+          ENDDO
+        ENDDO
+      ENDDO          
+      
+      DO iq=1,nqtot
+        DO l=ll_begin,ll_end
+          DO j=jj_begin,jj_end
+            DO i=ii_begin,ii_end
+              ij=(j-1)*iim+i
+              q(ij,l,iq)=q(ij,l,iq)+ dtphy * dq(ij,l,iq)
+            ENDDO
+          ENDDO
+        ENDDO 
+      ENDDO
+
+!$OMP BARRIER
+      
+     IF (is_omp_first_level) THEN
+       DO j=jj_begin,jj_end
+         DO i=ii_begin,ii_end
+           ij=(j-1)*iim+i
+           ps(ij)=ps(ij)+ dtphy * dps(ij)
+          ENDDO
+       ENDDO
+     ENDIF
+
+!     CALL compute_temperature2theta_rhodz(ps,Temp,theta_rhodz,0)
+
+! compute pression
+!$OMP BARRIER
+      DO    l    = ll_begin,ll_endp1
+        DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+            ij=(j-1)*iim+i
+            p(ij,l) = ap(l) + bp(l) * ps(ij)
+          ENDDO
+        ENDDO
+      ENDDO
+
+!$OMP BARRIER
+
+! compute exner
+       
+       IF (is_omp_first_level) THEN
+         DO j=jj_begin,jj_end
+            DO i=ii_begin,ii_end
+               ij=(j-1)*iim+i
+               pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
+            ENDDO
+         ENDDO
+       ENDIF
+
+       ! 3D : pk
+       DO l = ll_begin,ll_end
+          DO j=jj_begin,jj_end
+             DO i=ii_begin,ii_end
+                ij=(j-1)*iim+i
+                pk(ij,l) = cpp * ((.5/preff)*(p(ij,l)+p(ij,l+1))) ** kappa
+             ENDDO
+          ENDDO
+       ENDDO
+
+!$OMP BARRIER
+
+!   compute theta, temperature and pression at layer
+    DO    l    = ll_begin, ll_end
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          theta_rhodz(ij,l) = temp(ij,l) * ((p(ij,l)-p(ij,l+1))/g) / (pk(ij,l) / cpp )
+        ENDDO
+      ENDDO
+    ENDDO
+    
+    END SUBROUTINE grid_physics_to_icosa
+    
+    
+    SUBROUTINE grid_icosa_to_physics
+    USE pression_mod
+    USE exner_mod
+    USE theta2theta_rhodz_mod
+    USE geopotential_mod
+    USE wind_mod
+    USE omp_para
+    IMPLICIT NONE
+    
+    REAL(rstd) :: uc(3)
+    INTEGER :: i,j,ij,l
+    
+
+! compute pression
+
+      DO    l    = ll_begin,ll_endp1
+        DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+            ij=(j-1)*iim+i
+            p(ij,l) = ap(l) + bp(l) * ps(ij)
+          ENDDO
+        ENDDO
+      ENDDO
+
+!$OMP BARRIER
+
+! compute exner
+       
+       IF (is_omp_first_level) THEN
+         DO j=jj_begin,jj_end
+            DO i=ii_begin,ii_end
+               ij=(j-1)*iim+i
+               pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
+            ENDDO
+         ENDDO
+       ENDIF
+
+       ! 3D : pk
+       DO l = ll_begin,ll_end
+          DO j=jj_begin,jj_end
+             DO i=ii_begin,ii_end
+                ij=(j-1)*iim+i
+                pk(ij,l) = cpp * ((.5/preff)*(p(ij,l)+p(ij,l+1))) ** kappa
+             ENDDO
+          ENDDO
+       ENDDO
+
+!$OMP BARRIER
+
+!   compute theta, temperature and pression at layer
+    DO    l    = ll_begin, ll_end
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          theta(ij,l) = theta_rhodz(ij,l) / ((p(ij,l)-p(ij,l+1))/g)
+          Temp(ij,l) = theta(ij,l) * pk(ij,l) / cpp
+          p_layer(ij,l)=preff*(pk(ij,l)/cpp)**(1./kappa) 
+        ENDDO
+      ENDDO
+    ENDDO
+
+
+!!! Compute geopotential
+       
+  ! for first layer
+  IF (is_omp_first_level) THEN
+    DO j=jj_begin,jj_end
+      DO i=ii_begin,ii_end
+        ij=(j-1)*iim+i
+        phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) )
+      ENDDO
+    ENDDO
+  ENDIF
+!!-> implicit flush on phi(:,1)
+          
+  ! for other layers
+  DO l = ll_beginp1, ll_end
+    DO j=jj_begin,jj_end
+      DO i=ii_begin,ii_end
+        ij=(j-1)*iim+i
+        phi(ij,l) =  0.5 * ( theta(ij,l)  + theta(ij,l-1) )  & 
+                         * (  pk(ij,l-1) -  pk(ij,l)    )
+      ENDDO
+    ENDDO
+  ENDDO       
+
+!$OMP BARRIER
+
+
+  IF (is_omp_first_level) THEN
+    DO l = 2, llm
+      DO j=jj_begin,jj_end
+! ---> Bug compilo intel ici en openmp
+! ---> Couper la boucle
+       IF (j==jj_end+1) PRINT*,"this message must not be printed"
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          phi(ij,l) = phi(ij,l)+ phi(ij,l-1)
+        ENDDO
+      ENDDO
+    ENDDO
+! --> IMPLICIT FLUSH on phi --> non
+  ENDIF 
+
+! compute wind centered lon lat compound
+    DO l=ll_begin,ll_end
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          uc(:)=1/Ai(ij)*                                                                                                &
+                        ( ne(ij,right)*u(ij+u_right,l)*le(ij+u_right)*((xyz_v(ij+z_rdown,:)+xyz_v(ij+z_rup,:))/2-centroid(ij,:))  &
+                         + ne(ij,rup)*u(ij+u_rup,l)*le(ij+u_rup)*((xyz_v(ij+z_rup,:)+xyz_v(ij+z_up,:))/2-centroid(ij,:))          &
+                         + ne(ij,lup)*u(ij+u_lup,l)*le(ij+u_lup)*((xyz_v(ij+z_up,:)+xyz_v(ij+z_lup,:))/2-centroid(ij,:))          &
+                         + ne(ij,left)*u(ij+u_left,l)*le(ij+u_left)*((xyz_v(ij+z_lup,:)+xyz_v(ij+z_ldown,:))/2-centroid(ij,:))    &
+                         + ne(ij,ldown)*u(ij+u_ldown,l)*le(ij+u_ldown)*((xyz_v(ij+z_ldown,:)+xyz_v(ij+z_down,:))/2-centroid(ij,:))&
+                         + ne(ij,rdown)*u(ij+u_rdown,l)*le(ij+u_rdown)*((xyz_v(ij+z_down,:)+xyz_v(ij+z_rdown,:))/2-centroid(ij,:)))
+          ulon(ij,l)=sum(uc(:)*elon_i(ij,:))
+          ulat(ij,l)=sum(uc(:)*elat_i(ij,:)) 
+        ENDDO
+      ENDDO
+    ENDDO
+
+!$OMP BARRIER
+
+
+      
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          offset=offset+1
+
+          IF (is_omp_first_level) ps_phy(offset) = ps(ij)
+          p_phy(offset,ll_begin:ll_endp1) = p(ij,ll_begin:ll_endp1)
+          p_layer_phy(offset,ll_begin:ll_end) = p_layer(ij,ll_begin:ll_end)
+          Temp_phy(offset,ll_begin:ll_end) = Temp(ij,ll_begin:ll_end)
+          IF (is_omp_first_level) phis_phy(offset) = phis(ij)
+          phi_phy(offset,ll_begin:ll_end) = phi(ij,ll_begin:ll_end)-phis(ij)
+          ulon_phy(offset,ll_begin:ll_end) = ulon(ij,ll_begin:ll_end)
+          ulat_phy(offset,ll_begin:ll_end) = ulat(ij,ll_begin:ll_end)
+          q_phy(offset,ll_begin:ll_end,:) = q(ij,ll_begin:ll_end,:)
+          wflux_phy(offset,ll_begin:ll_end) = wflux(ij,ll_begin:ll_end)
+        ENDDO
+      ENDDO
+     
+     END SUBROUTINE grid_icosa_to_physics
+
+  END SUBROUTINE physics
+  
+  
+END MODULE physics_lmdz_generic_mod
+    
+  
Index: codes/icosagcm/trunk/src/sponge.f90
===================================================================
--- codes/icosagcm/trunk/src/sponge.f90	(revision 328)
+++ codes/icosagcm/trunk/src/sponge.f90	(revision 328)
@@ -0,0 +1,148 @@
+MODULE sponge_mod
+  USE icosa
+
+  PRIVATE
+  
+  REAL,SAVE :: tau_sponge !inverse of charactericstic relaxation time scale at the topmost layer (Hz)
+  INTEGER,SAVE :: iflag_sponge !0 --> for no sponge
+			       !1 --> for sponge over 4 topmost layers
+			       !2 --> for sponge from top to ~1% of top layer pressure
+  INTEGER,SAVE :: mode_sponge !1 --> u and v relax towards 0
+                              !2 --> u and v relax towards their zonal mean
+                              !3 --> u,v and pot. temp. relax towards their zonal mean
+!$OMP THREADPRIVATE(tau_sponge,iflag_sponge,mode_sponge)
+  REAL,ALLOCATABLE,SAVE :: rdamp(:) ! quenching coefficient
+  REAL,ALLOCATABLE,SAVE:: lambda(:) ! inverse or quenching time scale (Hz)
+
+  PUBLIC sponge, init_sponge, iflag_sponge
+
+CONTAINS
+
+  SUBROUTINE init_sponge
+  USE icosa
+  USE disvert_mod
+  USE omp_para
+  USE mpipara, ONLY: is_mpi_master
+  IMPLICIT NONE
+    INTEGER :: l
+
+    tau_sponge = 1.e-4
+    CALL getin("tau_sponge",tau_sponge)
+    PRINT*,'tau_sponge = ',tau_sponge
+
+    iflag_sponge = 0
+    CALL getin("iflag_sponge",iflag_sponge)
+    PRINT*,'iflag_sponge = ',iflag_sponge
+    
+    mode_sponge = 1
+    CALL getin("mode_sponge",mode_sponge)
+    PRINT*,'mode_sponge = ',mode_sponge
+
+    IF (iflag_sponge == 0) THEN
+      PRINT*,'init_sponge: no sponge'
+      RETURN
+    ENDIF
+
+!$OMP MASTER
+    ALLOCATE(rdamp(llm))
+    ALLOCATE(lambda(llm))
+
+    IF (iflag_sponge == 1) THEN
+! sponge quenching over the topmost 4 atmospheric layers
+        lambda(:)=0.
+        lambda(llm)=tau_sponge
+        lambda(llm-1)=tau_sponge/2.
+        lambda(llm-2)=tau_sponge/4.
+        lambda(llm-3)=tau_sponge/8.
+    ELSE IF (iflag_sponge == 2) THEN
+! sponge quenching over topmost layers down to pressures which are
+! higher than 100 times the topmost layer pressure
+        lambda(:)=tau_sponge*max(presnivs(llm)/presnivs(:)-0.01,0.)
+    ELSE
+      PRINT*,'Bad selector for variable iflag_sponge : <',iflag_sponge,             &
+            '> options are 0,1,2'
+      STOP
+    ENDIF
+
+! quenching coefficient rdamp(:)
+!         rdamp(:)=dt*lambda(:) ! Explicit Euler approx.
+!      rdamp(:)=1.-exp(-lambda(:)*dt)
+    rdamp(:)=itau_dissip*lambda(:)
+
+
+    IF (is_mpi_master) THEN
+      PRINT*,'lambda = '
+      DO l=ll_begin,ll_end
+        PRINT*,l,lambda(l)
+      ENDDO
+      PRINT*,'rdamp = '
+      DO l=ll_begin,ll_end
+        PRINT*,l,rdamp(l)
+      ENDDO
+    ENDIF
+!$OMP END MASTER
+!$OMP BARRIER
+
+  END SUBROUTINE init_sponge
+
+  SUBROUTINE sponge(f_ue,f_due,f_theta_rhodz,f_dtheta_rhodz)
+  USE icosa
+  USE theta2theta_rhodz_mod
+  USE pression_mod
+  USE exner_mod
+  USE geopotential_mod
+  USE trace
+  USE time_mod
+  USE omp_para
+  IMPLICIT NONE
+    TYPE(t_field),POINTER :: f_ue(:)
+    TYPE(t_field),POINTER :: f_due(:)
+    TYPE(t_field),POINTER :: f_theta_rhodz(:)
+    TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
+
+    REAL(rstd),POINTER         :: due(:,:)!,due_sponge(:,:)
+    REAL(rstd),POINTER         :: ue(:,:)
+    REAL(rstd),POINTER         :: theta_rhodz(:,:)
+    REAL(rstd),POINTER         :: dtheta_rhodz(:,:)!,dtheta_sponge(:,:)
+
+    INTEGER :: ind
+    INTEGER :: l,ij
+
+!$OMP BARRIER
+    
+    CALL trace_start("sponge")
+
+    IF (mode_sponge == 1) THEN
+      DO ind=1,ndomain
+        IF (.NOT. assigned_domain(ind)) CYCLE
+        CALL swap_dimensions(ind)
+        CALL swap_geometry(ind)
+        ue=f_ue(ind) 
+        due=f_due(ind) 
+        theta_rhodz=f_theta_rhodz(ind)
+        dtheta_rhodz=f_dtheta_rhodz(ind)
+
+        DO l=ll_begin,ll_end
+!$SIMD
+          DO ij=ij_begin,ij_end
+
+            due(ij+u_right,l) = -rdamp(l)*ue(ij+u_right,l)
+            due(ij+u_lup,l)   = -rdamp(l)*ue(ij+u_lup,l)
+            due(ij+u_ldown,l) = -rdamp(l)*ue(ij+u_ldown,l)
+
+            dtheta_rhodz(ij,l) = 0.0
+          ENDDO
+        ENDDO
+      END DO
+    ELSE
+      PRINT*,'Bad selector for variable mode_sponge : <',mode_sponge,             &
+            '> options 2 and 3 not available for the moment!'
+      STOP
+    ENDIF
+
+   CALL trace_end("sponge")
+
+!$OMP BARRIER
+  END SUBROUTINE sponge
+
+END MODULE sponge_mod
