Index: codes/icosagcm/trunk/src/advect_tracer.f90
===================================================================
--- codes/icosagcm/trunk/src/advect_tracer.f90	(revision 148)
+++ codes/icosagcm/trunk/src/advect_tracer.f90	(revision 149)
@@ -66,5 +66,5 @@
     CALL transfert_request(f_rhodz,req_i1)
 	
-    IF (is_mpi_root) PRINT *, 'Advection scheme'
+!    IF (is_mpi_root) PRINT *, 'Advection scheme'
 
 !    DO ind=1,ndomain
Index: codes/icosagcm/trunk/src/caldyn_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 148)
+++ codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 149)
@@ -11,5 +11,5 @@
   TYPE(t_field),POINTER :: f_buf_v(:), f_buf_s(:), f_buf_p(:)
 
-  PUBLIC init_caldyn, caldyn, write_output_fields
+  PUBLIC init_caldyn, caldyn, write_output_fields,un2ulonlat
 
   INTEGER :: caldyn_hydrostat, caldyn_conserv
@@ -722,15 +722,15 @@
     
     CALL writefield("ps",f_ps)
-    CALL writefield("dps",f_dps)
-    CALL writefield("phis",f_phis)
-    CALL vorticity(f_u,f_buf_v)
-    CALL writefield("vort",f_buf_v)
-
-    CALL w_omega(f_ps, f_u, f_buf_i)
-    CALL writefield('omega', f_buf_i)
-    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
-      CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
-      CALL writefield("omega"//TRIM(str_pression),f_buf_s)
-    ENDIF
+!    CALL writefield("dps",f_dps)
+!    CALL writefield("phis",f_phis)
+!    CALL vorticity(f_u,f_buf_v)
+!    CALL writefield("vort",f_buf_v)
+
+!    CALL w_omega(f_ps, f_u, f_buf_i)
+!    CALL writefield('omega', f_buf_i)
+!    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
+!      CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
+!      CALL writefield("omega"//TRIM(str_pression),f_buf_s)
+!    ENDIF
     
     ! Temperature
@@ -767,8 +767,8 @@
     ! geopotential
     CALL thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_buf_s,f_buf_p,f_buf1_i,f_buf2_i,f_buf_i)
-    CALL writefield("p",f_buf_p)
-    CALL writefield("phi",f_buf_i)
-    CALL writefield("theta",f_buf1_i) ! potential temperature
-    CALL writefield("pk",f_buf2_i) ! Exner pressure
+!    CALL writefield("p",f_buf_p)
+!    CALL writefield("phi",f_buf_i)
+     CALL writefield("theta",f_buf1_i) ! potential temperature
+!    CALL writefield("pk",f_buf2_i) ! Exner pressure
   
   
Index: codes/icosagcm/trunk/src/check_conserve.f90
===================================================================
--- codes/icosagcm/trunk/src/check_conserve.f90	(revision 149)
+++ codes/icosagcm/trunk/src/check_conserve.f90	(revision 149)
@@ -0,0 +1,314 @@
+MODULE check_conserve_mod
+	USE icosa 
+	IMPLICIT NONE 
+
+	PRIVATE 
+     TYPE(t_field),POINTER,SAVE :: f_pk(:), f_pks(:), f_p(:)
+     TYPE(t_field),POINTER,SAVE :: f_vort(:)
+     TYPE(t_field),POINTER,SAVE :: f_rhodz(:) 
+   
+     PUBLIC compute_conserve 
+     REAL(rstd),SAVE:: mtot0,ztot0,etot0,ang0,stot0,rmsv0
+     REAL(rstd),SAVE::rmsdpdt,etot,ang,stot,rmsv
+     REAL(rstd),SAVE:: ztot,mtot
+	
+     
+Contains 
+
+!---------------------------------------------------------------------
+
+   SUBROUTINE allocate_check_conserve
+    USE icosa
+    IMPLICIT NONE
+    CALL allocate_field(f_pk,field_t,type_real,llm)
+    CALL allocate_field(f_p,field_t,type_real,llm+1)
+    CALL allocate_field(f_pks,field_t,type_real)
+    CALL allocate_field(f_rhodz,field_t,type_real,llm)
+    CALL allocate_field(f_vort,field_z,type_real,llm)
+   END SUBROUTINE allocate_check_conserve
+
+!---------------------------------------------------------------------
+
+   SUBROUTINE check_mass_conserve(f_ps,f_dps)
+  USE icosa
+  IMPLICIT NONE
+    TYPE(t_field),POINTER :: f_ps(:)
+    TYPE(t_field),POINTER :: f_dps(:)
+    REAL(rstd),POINTER :: ps(:),dps(:) 
+    INTEGER :: ind,i,j,ij  
+
+   DO ind=1,ndomain
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      dps=f_dps(ind) 
+      ps=f_ps(ind)
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          IF (domain(ind)%own(i,j)) THEN
+            mtot=mtot+ps(ij)*Ai(ij)
+	       rmsdpdt=rmsdpdt+dps(ij)*dps(ij)
+	     ENDIF
+        ENDDO
+      ENDDO
+    ENDDO
+  END SUBROUTINE check_mass_conserve
+
+!---------------------------------------------------------------------
+
+   SUBROUTINE check_EN(f_ue,f_theta_rhodz,f_phis) 
+    USE icosa
+    USE pression_mod
+    USE vorticity_mod
+    IMPLICIT NONE
+    TYPE(t_field), POINTER :: f_ue(:)
+    TYPE(t_field), POINTER :: f_theta_rhodz(:)
+    TYPE(t_field), POINTER :: f_phis(:)
+    REAL(rstd), POINTER :: ue(:,:)
+    REAL(rstd), POINTER :: p(:,:) 
+    REAL(rstd), POINTER :: theta_rhodz(:,:) 
+    REAL(rstd), POINTER :: pk(:,:) 
+    REAL(rstd), POINTER :: phis(:) 
+    REAL(rstd), POINTER :: rhodz(:,:) 
+    INTEGER :: ind
+
+
+      Do ind=1,ndomain 
+        CALL swap_dimensions(ind)
+        CALL swap_geometry(ind)
+	 ue=f_ue(ind) 
+         p=f_p(ind)
+         rhodz=f_rhodz(ind) 
+         theta_rhodz=f_theta_rhodz(ind) 
+	    pk=f_pk(ind) 
+	    phis=f_phis(ind) 
+        CALL compute_en(ind,ue,p,rhodz,theta_rhodz,pk,phis)
+	 END DO 
+  END SUBROUTINE CHECK_EN 
+
+   SUBROUTINE compute_en(ind,u,p,rhodz,theta_rhodz,pk,phis)
+	 USE icosa
+      USE disvert_mod
+      USE wind_mod
+     IMPLICIT NONE
+    INTEGER,INTENT(IN)::ind 
+    REAL(rstd),INTENT(IN) :: u(3*iim*jjm,llm)
+    REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1)
+    REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(IN) :: theta_rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(IN) :: pk(iim*jjm,llm)
+    REAL(rstd),INTENT(IN) :: phis(iim*jjm)
+    REAL(rstd):: theta(iim*jjm,llm)
+    REAL(rstd)::KE(iim*jjm,llm) 
+    REAL(rstd) :: ucenter(iim*jjm,3,llm)
+    REAL(rstd) :: ulon(iim*jjm,llm)
+    REAL(rstd) :: ulat(iim*jjm,llm)
+
+    REAL(rstd)::masse(iim*jjm,llm)  
+    REAL(rstd)::rad,radomeg,lat,lon
+    INTEGER :: i,j,ij,l,ij2
+
+  
+       DO l = 1, llm
+        DO j=jj_begin,jj_end
+         DO i=ii_begin,ii_end
+	        ij=(j-1)*iim+i
+	         masse(ij,l) = rhodz(ij,l)*Ai(ij) 
+              theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
+             IF (domain(ind)%own(i,j)) THEN
+	         stot = stot + Ai(ij)*theta_rhodz(ij,l)
+             END IF 
+          END DO 
+         END DO
+	   END DO 
+
+       DO l = 1, llm
+        DO j=jj_begin,jj_end
+         DO i=ii_begin,ii_end
+	        ij=(j-1)*iim+i
+           IF (domain(ind)%own(i,j)) THEN
+	        KE(ij,l)=   1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +   &
+                                     le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +           &
+                                     le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +           &
+                                     le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +        &
+                                     le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +     &
+                                     le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )  
+             rmsv = rmsv + 2*masse(ij,l)*KE(ij,l) 
+            END IF 
+	     END DO
+         END DO 
+        END DO  
+
+       DO l = 1, llm
+        DO j=jj_begin,jj_end
+         DO i=ii_begin,ii_end
+	        ij=(j-1)*iim+i
+           IF (domain(ind)%own(i,j)) THEN
+	      etot = etot + masse(ij,l)*(phis(ij)+theta(ij,l)*pk(ij,l)+KE(ij,l))
+           END IF 
+          END DO 
+        END DO 
+       END DO 
+
+!------------------------------ ANGULAR VELOCITY 
+	CALL compute_wind_centered(u,ucenter)
+     CALL compute_wind_centered_lonlat_compound(ucenter,ulon,ulat)
+	 rad = radius 
+      radomeg = rad*omega
+
+      DO l = 1, llm
+        DO j=jj_begin,jj_end
+         DO i=ii_begin,ii_end
+	        ij=(j-1)*iim+i
+           IF (domain(ind)%own(i,j)) THEN
+             CALL xyz2lonlat(xyz_i(ij,:),lon,lat) 
+                ang=ang+rad*cos(lat)*masse(ij,l)*(ulon(ij,l)+ radomeg*cos(lat))
+           END IF 
+         END DO
+       END DO
+      END DO    
+          
+	END SUBROUTINE compute_en
+
+!---------------------------------------------------------------------
+
+   SUBROUTINE check_PV 
+    USE icosa
+    IMPLICIT NONE
+
+    REAL(rstd), POINTER :: vort(:,:)
+    REAL(rstd), POINTER :: rhodz(:,:) 
+    INTEGER :: ind
+
+    Do ind=1,ndomain 
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      vort=f_vort(ind)
+      rhodz=f_rhodz(ind) 
+      CALL compute_PV(vort,rhodz) 
+    ENDDO
+ 
+   END SUBROUTINE check_PV 
+
+	SUBROUTINE compute_PV(vort,rhodz) 
+	  USE icosa
+      USE disvert_mod
+     IMPLICIT NONE 
+    REAL(rstd),INTENT(IN) :: vort(2*iim*jjm,llm)
+    REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm) 
+    REAL(rstd)::qv1,qv2 
+    REAL(rstd)::hv1,hv2  
+    INTEGER :: i,j,ij,l,ij2
+
+	hv1 = 0.0 ; hv2 = 0.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
+            hv1 =  Riv2(ij,vup)        * rhodz(ij,l)           &
+              + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     &
+              + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
+      
+            qv1 = ( vort(ij+z_up,l)+fv(ij+z_up) )/hv1
+            
+            hv2  =  Riv2(ij,vdown)     * rhodz(ij,l)          &
+              + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  &
+              + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
+                       
+            qv2 =( vort(ij+z_down,l)+fv(ij+z_down) )/hv2
+
+            ztot = ztot + (qv1*qv1*hv1 + qv2*qv2*hv2) 
+	  ENDDO
+         ENDDO
+        ENDDO
+    END SUBROUTINE compute_PV
+!---------------------------------------------------------------------
+
+    SUBROUTINE compute_rhodz(p,rhodz) 
+	USE icosa 
+	IMPLICIT NONE 
+     REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1)
+     REAL(rstd),INTENT(OUT):: rhodz(iim*jjm,llm)
+     REAL(rstd) ::mass(iim*jjm,llm)
+     INTEGER :: i,j,ij,l 
+	
+     DO l = 1, llm
+      DO j=jj_begin,jj_end
+       DO i=ii_begin,ii_end
+         ij=(j-1)*iim+i
+         mass(ij,l) = ( p(ij,l) - p(ij,l+1) ) * Ai(ij)/g
+         rhodz(ij,l) = mass(ij,l) / Ai(ij)
+       ENDDO
+      ENDDO
+     ENDDO
+
+	END SUBROUTINE compute_rhodz
+!====================================================================
+
+	SUBROUTINE compute_conserve(f_ps,f_dps,f_ue,f_theta_rhodz,f_phis,it) 
+	 USE icosa 
+      USE pression_mod
+      USE vorticity_mod
+      USE caldyn_gcm_mod
+      USE exner_mod 
+      USE mpipara
+     TYPE(t_field),POINTER :: f_ps(:)
+     TYPE(t_field),POINTER :: f_dps(:)
+     TYPE(t_field),POINTER :: f_ue(:)
+     TYPE(t_field),POINTER :: f_theta_rhodz(:)
+     TYPE(t_field),POINTER :: f_phis(:)
+     INTEGER::it 
+     REAL(rstd),POINTER :: p(:,:),rhodz(:,:) 
+     INTEGER::ind 
+
+     etot=0.0; ang=0.0;stot=0.0;rmsv=0.0 
+     mtot=0.0 ; rmsdpdt=0.0 ; ztot = 0.0 
+
+     CALL allocate_check_conserve
+
+     CALL pression(f_ps,f_p)
+
+     Do ind=1,ndomain 
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+         p=f_p(ind) 
+        rhodz=f_rhodz(ind)
+       CALL compute_rhodz(p,rhodz) 
+      END DO 
+
+     CALL vorticity(f_ue,f_vort)
+     CALL check_mass_conserve(f_ps,f_dps)
+     CALL check_PV 
+     CALL exner(f_ps,f_p,f_pks,f_pk)
+     CALL check_EN(f_ue,f_theta_rhodz,f_phis) 
+    	
+	   IF ( it == 0  ) Then 
+	   ztot0 = ztot
+           mtot0 = mtot
+           etot0 = etot 
+           ang0  = ang 
+           stot0 = stot 
+        END IF 
+
+      rmsv=SQRT(rmsv/mtot) 
+      ztot=ztot/ztot0 ; mtot=mtot/mtot0 
+      etot=etot/etot0 ; ang=ang/ang0 ; stot=stot/stot0 
+      rmsdpdt= daysec*1.e-2*sqrt(rmsdpdt/ncell_glo)  
+
+
+	  IF (is_mpi_root) THEN 
+       OPEN(134,file="checkconsicosa.txt",position='append')
+       write(134,4000)mtot,rmsdpdt,etot,ztot,stot,rmsv,ang
+       write(134,*)mtot,rmsdpdt,etot,ztot,stot,rmsv,ang 
+       write(134,*)"=================================================="
+       write(*,4000)mtot,rmsdpdt,etot,ztot,stot,rmsv,ang
+
+4000   FORMAT(10x,'masse',4x,'rmsdpdt',7x,'energie',2x,'enstrophie'  &
+      ,2x,'entropie',3x,'rmsv',4x,'mt.ang',/,'GLOB  '                &
+       ,f10.6,e13.6,5f10.3/)     
+       close(134) 
+       END IF 
+    END SUBROUTINE compute_conserve
+
+ END MODULE check_conserve_mod
Index: codes/icosagcm/trunk/src/check_conserve.f90~
===================================================================
--- codes/icosagcm/trunk/src/check_conserve.f90~	(revision 149)
+++ codes/icosagcm/trunk/src/check_conserve.f90~	(revision 149)
@@ -0,0 +1,315 @@
+MODULE check_conserve_mod
+	USE icosa 
+	IMPLICIT NONE 
+
+	PRIVATE 
+     TYPE(t_field),POINTER,SAVE :: f_pk(:), f_pks(:), f_p(:)
+     TYPE(t_field),POINTER,SAVE :: f_vort(:)
+     TYPE(t_field),POINTER,SAVE :: f_rhodz(:) 
+   
+     PUBLIC compute_conserve 
+     REAL(rstd),SAVE:: mtot0,ztot0,etot0,ang0,stot0,rmsv0
+     REAL(rstd),SAVE::rmsdpdt,etot,ang,stot,rmsv
+     REAL(rstd),SAVE:: ztot,mtot
+	
+     
+Contains 
+
+!---------------------------------------------------------------------
+
+   SUBROUTINE allocate_check_conserve
+    USE icosa
+    IMPLICIT NONE
+    CALL allocate_field(f_pk,field_t,type_real,llm)
+    CALL allocate_field(f_p,field_t,type_real,llm+1)
+    CALL allocate_field(f_pks,field_t,type_real)
+    CALL allocate_field(f_rhodz,field_t,type_real,llm)
+    CALL allocate_field(f_vort,field_z,type_real,llm)
+   END SUBROUTINE allocate_check_conserve
+
+!---------------------------------------------------------------------
+
+   SUBROUTINE check_mass_conserve(f_ps,f_dps)
+  USE icosa
+  IMPLICIT NONE
+    TYPE(t_field),POINTER :: f_ps(:)
+    TYPE(t_field),POINTER :: f_dps(:)
+    REAL(rstd),POINTER :: ps(:),dps(:) 
+    INTEGER :: ind,i,j,ij  
+
+   DO ind=1,ndomain
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      dps=f_dps(ind) 
+      ps=f_ps(ind)
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          IF (domain(ind)%own(i,j)) THEN
+            mtot=mtot+ps(ij)*Ai(ij)
+	       rmsdpdt=rmsdpdt+dps(ij)*dps(ij)
+	     ENDIF
+        ENDDO
+      ENDDO
+    ENDDO
+  END SUBROUTINE check_mass_conserve
+
+!---------------------------------------------------------------------
+
+   SUBROUTINE check_EN(f_ue,f_theta_rhodz,f_phis) 
+    USE icosa
+    USE pression_mod
+    USE vorticity_mod
+    IMPLICIT NONE
+    TYPE(t_field), POINTER :: f_ue(:)
+    TYPE(t_field), POINTER :: f_theta_rhodz(:)
+    TYPE(t_field), POINTER :: f_phis(:)
+    REAL(rstd), POINTER :: ue(:,:)
+    REAL(rstd), POINTER :: p(:,:) 
+    REAL(rstd), POINTER :: theta_rhodz(:,:) 
+    REAL(rstd), POINTER :: pk(:,:) 
+    REAL(rstd), POINTER :: phis(:) 
+    REAL(rstd), POINTER :: rhodz(:,:) 
+    INTEGER :: ind
+
+
+      Do ind=1,ndomain 
+        CALL swap_dimensions(ind)
+        CALL swap_geometry(ind)
+	    ue=f_ue(ind) 
+         p=f_p(ind)
+         rhodz=f_rhodz(ind) 
+         theta_rhodz=f_theta_rhodz(ind) 
+	    pk=f_pk(ind) 
+	    phis=f_phis(ind) 
+        CALL compute_en(ind,ue,p,rhodz,theta_rhodz,pk,phis)
+	 END DO 
+  END SUBROUTINE CHECK_EN 
+
+   SUBROUTINE compute_en(ind,u,p,rhodz,theta_rhodz,pk,phis)
+	 USE icosa
+      USE disvert_mod
+      USE wind_mod
+     IMPLICIT NONE
+    INTEGER,INTENT(IN)::ind 
+    REAL(rstd),INTENT(IN) :: u(3*iim*jjm,llm)
+    REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1)
+    REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(IN) :: theta_rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(IN) :: pk(iim*jjm,llm)
+    REAL(rstd),INTENT(IN) :: phis(iim*jjm)
+    REAL(rstd):: theta(iim*jjm,llm)
+    REAL(rstd)::KE(iim*jjm,llm) 
+    REAL(rstd) :: ucenter(iim*jjm,3,llm)
+    REAL(rstd) :: ulon(iim*jjm,llm)
+    REAL(rstd) :: ulat(iim*jjm,llm)
+
+    REAL(rstd)::masse(iim*jjm,llm)  
+    REAL(rstd)::rad,radomeg,lat,lon
+    INTEGER :: i,j,ij,l,ij2
+
+  
+       DO l = 1, llm
+        DO j=jj_begin,jj_end
+         DO i=ii_begin,ii_end
+	        ij=(j-1)*iim+i
+	         masse(ij,l) = rhodz(ij,l)*Ai(ij) 
+              theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
+             IF (domain(ind)%own(i,j)) THEN
+	         stot = stot + Ai(ij)*theta_rhodz(ij,l)
+             END IF 
+          END DO 
+         END DO
+	   END DO 
+
+       DO l = 1, llm
+        DO j=jj_begin,jj_end
+         DO i=ii_begin,ii_end
+	        ij=(j-1)*iim+i
+           IF (domain(ind)%own(i,j)) THEN
+	        KE(ij,l)=   1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +   &
+                                     le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +           &
+                                     le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +           &
+                                     le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +        &
+                                     le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +     &
+                                     le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )  
+             rmsv = rmsv + 2*masse(ij,l)*KE(ij,l) 
+            END IF 
+	     END DO
+         END DO 
+        END DO  
+
+       DO l = 1, llm
+        DO j=jj_begin,jj_end
+         DO i=ii_begin,ii_end
+	        ij=(j-1)*iim+i
+           IF (domain(ind)%own(i,j)) THEN
+	      etot = etot + masse(ij,l)*(phis(ij)+theta(ij,l)*pk(ij,l)+KE(ij,l))
+           END IF 
+          END DO 
+        END DO 
+       END DO 
+
+!------------------------------ ANGULAR VELOCITY 
+	CALL compute_wind_centered(u,ucenter)
+     CALL compute_wind_centered_lonlat_compound(ucenter,ulon,ulat)
+	 rad = radius 
+      radomeg = rad*omega
+
+      DO l = 1, llm
+        DO j=jj_begin,jj_end
+         DO i=ii_begin,ii_end
+	        ij=(j-1)*iim+i
+           IF (domain(ind)%own(i,j)) THEN
+             CALL xyz2lonlat(xyz_i(ij,:),lon,lat) 
+                ang=ang+rad*cos(lat)*masse(ij,l)*(ulon(ij,l)+ radomeg*cos(lat))
+           END IF 
+         END DO
+       END DO
+      END DO    
+          
+	END SUBROUTINE compute_en
+
+!---------------------------------------------------------------------
+
+   SUBROUTINE check_PV 
+    USE icosa
+    IMPLICIT NONE
+
+    REAL(rstd), POINTER :: vort(:,:)
+    REAL(rstd), POINTER :: rhodz(:,:) 
+    INTEGER :: ind
+
+    Do ind=1,ndomain 
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      vort=f_vort(ind)
+      rhodz=f_rhodz(ind) 
+      CALL compute_PV(vort,rhodz) 
+    ENDDO
+ 
+   END SUBROUTINE check_PV 
+
+	SUBROUTINE compute_PV(vort,rhodz) 
+	  USE icosa
+      USE disvert_mod
+     IMPLICIT NONE 
+    REAL(rstd),INTENT(IN) :: vort(2*iim*jjm,llm)
+    REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm) 
+    REAL(rstd)::qv1,qv2 
+    REAL(rstd)::hv1,hv2  
+    INTEGER :: i,j,ij,l,ij2
+
+	hv1 = 0.0 ; hv2 = 0.0 
+     
+     DO l = 1,llm
+      DO j=jj_begin,jj_end
+       DO i=ii_begin,ii_end
+           ij=(j-1)*iim+i
+            hv1 =  Riv2(ij,vup)       * rhodz(ij,l)           &
+              + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     &
+              + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
+      
+            qv1 = ( vort(ij+z_up,l)+fv(ij+z_up) )/hv1
+            
+            hv2 =  Riv2(ij,vdown)     * rhodz(ij,l)          &
+              + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  &
+              + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
+                       
+            qv2 =( vort(ij+z_down,l)+fv(ij+z_down) )/hv2
+
+            ztot= ztot + (qv1*qv1*hv1 + qv2*qv2*hv2)
+	    ENDDO
+       ENDDO
+      ENDDO    
+    END SUBROUTINE compute_PV
+!---------------------------------------------------------------------
+
+    SUBROUTINE compute_rhodz(p,rhodz) 
+	USE icosa 
+	IMPLICIT NONE 
+     REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1)
+     REAL(rstd),INTENT(OUT):: rhodz(iim*jjm,llm)
+     REAL(rstd) ::mass(iim*jjm,llm)
+     INTEGER :: i,j,ij,l 
+	
+     DO l = 1, llm
+      DO j=jj_begin,jj_end
+       DO i=ii_begin,ii_end
+         ij=(j-1)*iim+i
+         mass(ij,l) = ( p(ij,l) - p(ij,l+1) ) * Ai(ij)/g
+         rhodz(ij,l) = mass(ij,l) / Ai(ij)
+       ENDDO
+      ENDDO
+     ENDDO
+
+	END SUBROUTINE compute_rhodz
+
+!=============================================================================
+
+	SUBROUTINE compute_conserve(f_ps,f_dps,f_ue,f_theta_rhodz,f_phis,it) 
+	 USE icosa 
+      USE pression_mod
+      USE vorticity_mod
+      USE caldyn_gcm_mod
+      USE exner_mod 
+      USE mpipara
+     TYPE(t_field),POINTER :: f_ps(:)
+     TYPE(t_field),POINTER :: f_dps(:)
+     TYPE(t_field),POINTER :: f_ue(:)
+     TYPE(t_field),POINTER :: f_theta_rhodz(:)
+     TYPE(t_field),POINTER :: f_phis(:)
+     INTEGER::it 
+     REAL(rstd),POINTER :: p(:,:),rhodz(:,:) 
+     INTEGER::ind 
+
+     etot=0.0; ang=0.0;stot=0.0;rmsv=0.0 
+     mtot=0.0 ; rmsdpdt=0.0 ; ztot = 0.0 
+
+     CALL allocate_check_conserve
+
+     CALL pression(f_ps,f_p)
+
+     Do ind=1,ndomain 
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+         p=f_p(ind) 
+        rhodz=f_rhodz(ind)
+       CALL compute_rhodz(p,rhodz) 
+      END DO 
+
+     CALL vorticity(f_ue,f_vort)
+     CALL check_mass_conserve(f_ps,f_dps)
+     CALL check_PV 
+     CALL exner(f_ps,f_p,f_pks,f_pk)
+     CALL check_EN(f_ue,f_theta_rhodz,f_phis) 
+    	
+	   IF ( it == 0  ) Then 
+		 ztot0 = ztot
+           mtot0 = mtot
+           etot0 = etot 
+           ang0  = ang 
+           stot0 = stot 
+        END IF 
+
+      rmsv=SQRT(rmsv/mtot) 
+	 ztot=ztot/ztot0 ; mtot=mtot/mtot0 
+      etot=etot/etot0 ; ang=ang/ang0 ; stot=stot/stot0 
+	 rmsdpdt= daysec*1.e-2*sqrt(rmsdpdt/ncell_glo)  
+
+
+	  IF (is_mpi_root) THEN 
+       OPEN(134,file="checkconsicosa.txt",position='append')
+       write(134,4000)mtot,rmsdpdt,etot,ztot,stot,rmsv,ang
+       write(134,*)mtot,rmsdpdt,etot,ztot,stot,rmsv,ang 
+       write(134,*)"=================================================="
+       write(*,4000)mtot,rmsdpdt,etot,ztot,stot,rmsv,ang
+
+4000   FORMAT(10x,'masse',4x,'rmsdpdt',7x,'energie',2x,'enstrophie'  &
+      ,2x,'entropie',3x,'rmsv',4x,'mt.ang',/,'GLOB  '                &
+       ,f10.6,e13.6,5f10.3/)     
+       close(134) 
+       END IF 
+    END SUBROUTINE compute_conserve
+
+ END MODULE check_conserve_mod
Index: codes/icosagcm/trunk/src/diffff
===================================================================
--- codes/icosagcm/trunk/src/diffff	(revision 149)
+++ codes/icosagcm/trunk/src/diffff	(revision 149)
@@ -0,0 +1,493 @@
+Index: time.f90
+===================================================================
+--- time.f90	(revision 142)
++++ time.f90	(working copy)
+@@ -10,11 +10,23 @@
+   REAL(rstd),SAVE :: write_period
+   INTEGER,SAVE    :: itau_out, itau_adv, itau_dissip, itau_physics, itaumax
+   
++  INTEGER,SAVE :: day_step,ndays
++  REAL(rstd),SAVE :: jD_ref,jH_ref
++  INTEGER,SAVE :: day_ini,day_end,annee_ref,day_ref
++  REAL(rstd),SAVE::start_time
++  CHARACTER(LEN=255) :: time_style
++  INTEGER,SAVE:: an, mois, jour
++  REAL(rstd),SAVE:: heure
++  CHARACTER (LEN=10):: calend
++
+   PUBLIC create_time_counter_header, update_time_counter, close_time_counter, init_time,  &
+-         dt, write_period, itau_out, itau_adv, itau_dissip, itau_physics, itaumax
++         dt, write_period, itau_out, itau_adv, itau_dissip, itau_physics, itaumax, & 
++day_step,ndays,jD_ref,jH_ref,day_ini,day_end,annee_ref,day_ref,an, mois, jour,heure, &
++            calend,time_style
+ 
+ 
+ 
++
+ CONTAINS
+   
+   SUBROUTINE init_time
+@@ -23,22 +35,18 @@
+   USE mpipara
+   IMPLICIT NONE
+   REAL(rstd) :: run_length
+-  
++
++
++     time_style='dcmip'
++   CALL getin('time_style',time_style)
++
++   IF (TRIM(time_style)=='dcmip')  Then
+     dt=90.
+     CALL getin('dt',dt)
+ 
+     itaumax=100
+     CALL getin('itaumax',itaumax)
+ 
+-    itau_adv=1
+-    CALL getin('itau_adv',itau_adv)
+-
+-    itau_dissip=1
+-    CALL getin('itau_dissip',itau_dissip)
+-
+-    itau_physics=1
+-    CALL getin('itau_physics',itau_physics)
+-
+     run_length=dt*itaumax
+     CALL getin('run_length',run_length)
+     itaumax=run_length/dt
+@@ -53,6 +61,17 @@
+     write_period=write_period/scale_factor
+     itau_out=FLOOR(.5+write_period/dt)
+     IF (is_mpi_root) PRINT *, 'Output frequency (scaled) set to ',write_period, ' : itau_out = ',itau_out
++    ENDIF 
++
++    itau_adv=1
++    CALL getin('itau_adv',itau_adv)
++
++    itau_dissip=1
++    CALL getin('itau_dissip',itau_dissip)
++
++    itau_physics=1
++    CALL getin('itau_physics',itau_physics)
++
+     
+     CALL create_time_counter_header
+     
+Index: caldyn_gcm.f90
+===================================================================
+--- caldyn_gcm.f90	(revision 142)
++++ caldyn_gcm.f90	(working copy)
+@@ -10,7 +10,7 @@
+   TYPE(t_field),POINTER :: f_buf_i(:), f_buf_ulon(:), f_buf_ulat(:), f_buf_u3d(:)
+   TYPE(t_field),POINTER :: f_buf_v(:), f_buf_s(:), f_buf_p(:)
+ 
+-  PUBLIC init_caldyn, caldyn, write_output_fields
++  PUBLIC init_caldyn, caldyn, write_output_fields,un2ulonlat
+ 
+   INTEGER :: caldyn_hydrostat, caldyn_conserv
+ 
+@@ -708,17 +708,17 @@
+     str_pression=ADJUSTL(str_pression)
+     
+     CALL writefield("ps",f_ps)
+-    CALL writefield("dps",f_dps)
+-    CALL writefield("phis",f_phis)
+-    CALL vorticity(f_u,f_buf_v)
+-    CALL writefield("vort",f_buf_v)
++!    CALL writefield("dps",f_dps)
++!    CALL writefield("phis",f_phis)
++!    CALL vorticity(f_u,f_buf_v)
++!    CALL writefield("vort",f_buf_v)
+ 
+-    CALL w_omega(f_ps, f_u, f_buf_i)
+-    CALL writefield('omega', f_buf_i)
+-    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
+-      CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
+-      CALL writefield("omega"//TRIM(str_pression),f_buf_s)
+-    ENDIF
++!    CALL w_omega(f_ps, f_u, f_buf_i)
++!    CALL writefield('omega', f_buf_i)
++!    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
++!      CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
++!      CALL writefield("omega"//TRIM(str_pression),f_buf_s)
++!    ENDIF
+     
+     ! Temperature
+     CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i) ;
+@@ -753,10 +753,10 @@
+     
+     ! geopotential
+     CALL thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_buf_s,f_buf_p,f_buf1_i,f_buf2_i,f_buf_i)
+-    CALL writefield("p",f_buf_p)
+-    CALL writefield("phi",f_buf_i)
+-    CALL writefield("theta",f_buf1_i) ! potential temperature
+-    CALL writefield("pk",f_buf2_i) ! Exner pressure
++!    CALL writefield("p",f_buf_p)
++!    CALL writefield("phi",f_buf_i)
++     CALL writefield("theta",f_buf1_i) ! potential temperature
++!    CALL writefield("pk",f_buf2_i) ! Exner pressure
+   
+   
+   END SUBROUTINE write_output_fields
+Index: advect_tracer.f90
+===================================================================
+--- advect_tracer.f90	(revision 142)
++++ advect_tracer.f90	(working copy)
+@@ -57,7 +57,7 @@
+     CALL transfert_request(f_q,req_i1)
+     CALL transfert_request(f_rhodz,req_i1)
+ 	
+-    IF (is_mpi_root) PRINT *, 'Advection scheme'
++!    IF (is_mpi_root) PRINT *, 'Advection scheme'
+ 
+ !    DO ind=1,ndomain
+ !       CALL swap_dimensions(ind)
+Index: dissip_gcm.f90
+===================================================================
+--- dissip_gcm.f90	(revision 142)
++++ dissip_gcm.f90	(working copy)
+@@ -93,7 +93,8 @@
+       rayleigh_friction_type=2
+       IF (is_mpi_root) PRINT *, 'Rayleigh friction : Schaer-like mountain with shear DCMIP2.2'
+    CASE DEFAULT
+-      IF (is_mpi_root) PRINT *, 'Bad selector : rayleigh_friction_type =', TRIM(rayleigh_friction_key), ' in dissip_gcm.f90/init_dissip'
++      IF (is_mpi_root) PRINT *, 'Bad selector : rayleigh_friction_type =', TRIM(rayleigh_friction_key), & 
++        ' in dissip_gcm.f90/init_dissip'
+       STOP
+    END SELECT
+ 
+Index: physics.f90
+===================================================================
+--- physics.f90	(revision 142)
++++ physics.f90	(working copy)
+@@ -7,7 +7,8 @@
+ 
+   SUBROUTINE init_physics
+   USE icosa
+-  USE physics_dcmip_mod, init_physics_dcmip=>init_physics
++  USE physics_dcmip_mod,init_physics_dcmip=>init_physics
++  USE physics_dry_mod
+   IMPLICIT NONE
+     
+     CALL getin("physics",physics_type)
+@@ -17,36 +18,57 @@
+     
+       CASE ('dcmip')
+         CALL init_physics_dcmip
++
++      CASE ('lmd')
++	CALL init_physics_dry
+       
+       CASE DEFAULT
+-         PRINT*, 'Bad selector for variable physics <',physics_type, &
++         PRINT*, 'Bad selector for variable physics init <',physics_type, &
+               '> options are <none>, <dcmip>,'
+-         STOP
++
+     END SELECT
+     
+   END SUBROUTINE init_physics
+ 
+-  SUBROUTINE physics(it,f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
++  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
+     INTEGER, INTENT(IN)   :: it
++    REAL(rstd),INTENT(IN)::jD_cur,jH_cur
+     TYPE(t_field),POINTER :: f_phis(:)
+     TYPE(t_field),POINTER :: f_ps(:)
+     TYPE(t_field),POINTER :: f_theta_rhodz(:)
+     TYPE(t_field),POINTER :: f_ue(:)
+     TYPE(t_field),POINTER :: f_q(:)
++    LOGICAL:: firstcall,lastcall
+     
+     SELECT CASE(TRIM(physics_type))
+       CASE ('none')
++
++	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)
++	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)
++
++      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
++  STOP
+     END SELECT
+     
+   END SUBROUTINE physics
+Index: etat0.f90
+===================================================================
+--- etat0.f90	(revision 142)
++++ etat0.f90	(working copy)
+@@ -1,4 +1,5 @@
+ MODULE etat0_mod
++    CHARACTER(len=255),SAVE :: etat0_type
+ 
+ CONTAINS
+   
+@@ -11,6 +12,10 @@
+     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 
++    USE dynetat0_hz_mod,  ONLY : dynetat0_hz=>etat0 
++
+     IMPLICIT NONE
+     TYPE(t_field),POINTER :: f_ps(:)
+     TYPE(t_field),POINTER :: f_phis(:)
+@@ -18,7 +23,6 @@
+     TYPE(t_field),POINTER :: f_u(:)
+     TYPE(t_field),POINTER :: f_q(:)
+     
+-    CHARACTER(len=255) :: etat0_type
+     etat0_type='jablonowsky06'
+     CALL getin("etat0",etat0_type)
+     
+@@ -27,6 +31,9 @@
+        CALL etat0_jablonowsky06(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
+     CASE ('academic')
+        CALL etat0_academic(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
++    CASE ('heldsz')
++	print*,"heldsz test case"
++       CALL etat0_heldsz(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
+     CASE ('dcmip1')
+        CALL etat0_dcmip1(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
+     CASE ('dcmip2_mountain','dcmip2_schaer_noshear','dcmip2_schaer_shear')
+@@ -37,6 +44,12 @@
+        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"	
++       CALL dynetat0_start(f_ps,f_phis,f_theta_rhodz,f_u,f_q) 
++	CASE ('readnf_hz') 
++	  print*,"readnf_hz used"
++       CALL dynetat0_hz(f_ps,f_phis,f_theta_rhodz,f_u,f_q) 
+    CASE DEFAULT
+        PRINT*, 'Bad selector for variable etat0 <',etat0_type, &
+             '> options are <jablonowsky06>, <academic>, <dcmip[1-4]> '
+Index: timeloop_gcm.f90
+===================================================================
+--- timeloop_gcm.f90	(revision 142)
++++ timeloop_gcm.f90	(working copy)
+@@ -17,7 +17,11 @@
+   USE advect_tracer_mod
+   USE physics_mod
+   USE mpipara
+-  
++!-------------------------------
++  USE IOIPSL
++  USE maxicosa
++  USE check_conserve_mod
++!------------------------------
+   IMPLICIT NONE
+   TYPE(t_field),POINTER :: f_phis(:)
+ !  TYPE(t_field),POINTER :: f_theta(:)
+@@ -46,6 +50,10 @@
+   INTEGER :: it,i,j,n, nb_stage, stage, matsuno_period, scheme
+   CHARACTER(len=255) :: scheme_name
+   LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time
++  CHARACTER(len=7) :: first 
++  REAL(rstd),SAVE :: jD_cur, jH_cur
++  REAL(rstd),SAVE :: start_time
++
+ !  INTEGER :: itaumax
+ !  REAL(rstd) ::write_period
+ !  INTEGER    :: itau_out
+@@ -84,6 +92,58 @@
+   CALL allocate_field(f_hfluxt,field_u,type_real,llm)   ! mass "fluxes" accumulated in time
+   CALL allocate_field(f_wfluxt,field_t,type_real,llm+1)
+ 
++!----------------------------------------------------
++  IF (TRIM(time_style)=='lmd')  Then
++
++   day_step=180
++   CALL getin('day_step',day_step)
++
++   ndays=1
++   CALL getin('ndays',ndays)
++
++   dt = daysec/REAL(day_step)
++   itaumax = ndays*day_step
++
++   calend = 'earth_360d'
++   CALL getin('calend', calend)
++
++   day_ini = 0
++   CALL getin('day_ini',day_ini)
++
++   day_end = 0
++   CALL getin('day_end',day_end)
++
++   annee_ref = 1998
++   CALL getin('annee_ref',annee_ref)
++
++   start_time = 0
++   CALL getin('start_time',start_time) 
++
++   write_period=0
++   CALL getin('write_period',write_period)
++      
++   write_period=write_period/scale_factor
++   itau_out=FLOOR(write_period/dt)
++   
++   PRINT *, 'Output frequency (scaled) set to ',write_period, ' : itau_out = ',itau_out 
++
++  mois = 1 ; heure = 0.
++  call ymds2ju(annee_ref, mois, day_ref, heure, jD_ref)
++  jH_ref = jD_ref - int(jD_ref) 
++  jD_ref = int(jD_ref) 
++  
++  CALL ioconf_startdate(INT(jD_ref),jH_ref) 
++  write(*,*)'annee_ref, mois, day_ref, heure, jD_ref'
++  write(*,*)annee_ref, mois, day_ref, heure, jD_ref
++  write(*,*)"ndays,day_step,itaumax,dt======>"
++  write(*,*)ndays,day_step,itaumax,dt
++  call ju2ymds(jD_ref+jH_ref,an, mois, jour, heure)
++  write(*,*)'jD_ref+jH_ref,an, mois, jour, heure'
++  write(*,*)jD_ref+jH_ref,an, mois, jour, heure
++  day_end = day_ini + ndays 
++ END IF 
++!----------------------------------------------------
++
+   scheme_name='runge_kutta'
+   CALL getin('scheme',scheme_name)
+ 
+@@ -122,7 +182,6 @@
+ !  CALL allocate_field(f_dum2,field_u,type_real,llm)
+ !  CALL allocate_field(f_dtheta_rhodzm1,field_t,type_real,llm)
+ !  CALL allocate_field(f_dtheta_rhodzm2,field_t,type_real,llm)
+-
+ !  CALL allocate_field(f_theta,field_t,type_real,llm)
+ !  CALL allocate_field(f_dtheta,field_t,type_real,llm)
+ 
+@@ -132,7 +191,11 @@
+   CALL init_advect_tracer
+   CALL init_physics
+   
+-  CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
++!============================ SET 0 to all field 
++  CALL initial0(f_ps,f_phis,f_theta_rhodz,f_u,f_q)
++! CALL dynetat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q)
++  CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q)
++!================================================
+   CALL writefield("phis",f_phis,once=.TRUE.)
+   CALL transfert_request(f_q,req_i1) 
+ 
+@@ -153,17 +216,18 @@
+      rhodz=f_rhodz(ind); ps=f_ps(ind)
+      CALL compute_rhodz(.FALSE., ps, rhodz)   
+   END DO
++
+   
+   DO it=0,itaumax
+ 
+-    IF (is_mpi_root) PRINT *,"It No :",It,"   t :",dt*It
+     IF (mod(it,itau_out)==0 ) THEN
+-      CALL writefield("q",f_q)
++!      IF (is_mpi_root) PRINT *,"It No :",It,"   t :",dt*It
+       CALL update_time_counter(dt*it)
++      CALL compute_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it)  
+     ENDIF
+-    
+-    CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q)
+ 
++      CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q)
++
+     DO stage=1,nb_stage
+        CALL caldyn((stage==1) .AND. (MOD(it,itau_out)==0), &
+             f_phis,f_ps,f_theta_rhodz,f_u, f_q, &
+@@ -177,11 +241,11 @@
+           CALL  leapfrog_matsuno_scheme(stage)
+           
+           !      CASE ('leapfrog')
+-          !        CALL leapfrog_scheme
++          !      CALL leapfrog_scheme
+           !
+           !      CASE ('adam_bashforth')
+-          !        CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz)
+-          !        CALL adam_bashforth_scheme
++          !      CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz)
++          !      CALL adam_bashforth_scheme
+        CASE DEFAULT
+           STOP
+        END SELECT
+@@ -192,7 +256,7 @@
+ 
+     IF(MOD(it+1,itau_adv)==0) THEN
+        CALL transfert_request(f_wfluxt,req_i1) ! FIXME
+-!       CALL transfert_request(f_hfluxt,req_e1) ! FIXME
++!      CALL transfert_request(f_hfluxt,req_e1) ! FIXME
+ 
+        CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz)  ! update q and rhodz after RK step
+        fluxt_zero=.TRUE.
+@@ -207,11 +271,17 @@
+           CALL swap_geometry(ind)
+           rhodz=f_rhodz(ind); ps=f_ps(ind); dps=f_dps(ind); 
+           wflux=f_wflux(ind); wfluxt=f_wfluxt(ind)
+-          CALL compute_rhodz(.FALSE., ps, rhodz)   
++!!!!!          CALL compute_rhodz(.FALSE., ps, rhodz)   
+        END DO
+-
+     END IF
+-
++!----------------------------------------------------
++  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)
++!	print*,"Just b4 phys"
++    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)
+     ENDDO
+  
+@@ -399,7 +469,7 @@
+        DO j=jj_begin-dd,jj_end+dd
+           DO i=ii_begin-dd,ii_end+dd
+              ij=(j-1)*iim+i
+-             m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g 
++             m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g
+              IF(comp) THEN
+                 rhodz(ij,l) = m
+              ELSE
+@@ -414,7 +484,7 @@
+           PRINT *, 'Discrepancy between ps and rhodz detected', err
+           STOP
+        ELSE
+-!          PRINT *, 'No discrepancy between ps and rhodz detected'
++          PRINT *, 'No discrepancy between ps and rhodz detected'
+        END IF
+     END IF
+ 
+Index: icosa_gcm.f90
+===================================================================
+--- icosa_gcm.f90	(revision 142)
++++ icosa_gcm.f90	(working copy)
+@@ -24,6 +24,7 @@
+   CALL compute_domain
+   CALL init_transfert
+   CALL init_writefield
++
+ !  CALL allocate_field(sum_ne,field_T,type_real)
+ !  CALL allocate_field_glo(sum_ne_glo,field_T,type_real)
+ !  
Index: codes/icosagcm/trunk/src/dimphys.f90
===================================================================
--- codes/icosagcm/trunk/src/dimphys.f90	(revision 149)
+++ codes/icosagcm/trunk/src/dimphys.f90	(revision 149)
@@ -0,0 +1,27 @@
+MODULE dimphys_mod
+        USE ICOSA
+      INTEGER::nsoilmx=10,offset
+      INTEGER::ngridmx,nlayermx
+      LOGICAL::callrad,calldifv,calladj,callcond,callsoil
+      LOGICAL::season,diurnal,lverbose
+      INTEGER:: iradia
+      REAL:: period_sort,timeperi
+      REAL:: aphelie,periheli,year_day,peri_day,obliquit
+      REAL::e_elips,p_elips,unitastr
+      REAL:: coefvis,coefir
+
+      REAL,ALLOCATABLE:: albedo(:),emissiv(:),inertie(:)
+      REAL,ALLOCATABLE:: tsurf(:),tsoil(:,:)
+      REAL,ALLOCATABLE:: z0(:),rnatur(:)
+      REAL,ALLOCATABLE:: dtrad(:,:) ,fluxrad(:)
+      REAL,ALLOCATABLE:: capcal(:)  ,fluxgrd(:)
+      REAL,ALLOCATABLE:: q2(:,:),q2l(:,:)
+
+
+!      REAL,allocatable:: albedo(:),emissiv(:),inertie(:)
+!      REAL,allocatable:: z0(:),rnatur(:)
+
+!      REAL,allocatable::tsurf(:),tsoil(:,:)
+!      REAL,allocatable::fluxrad(:),fluxgrd(:)
+!      REAL,ALLOCATABLE::capcal(:)
+END MODULE dimphys_mod
Index: codes/icosagcm/trunk/src/dissip_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/dissip_gcm.f90	(revision 148)
+++ codes/icosagcm/trunk/src/dissip_gcm.f90	(revision 149)
@@ -96,5 +96,6 @@
       IF (is_mpi_root) PRINT *, 'Rayleigh friction : Schaer-like mountain with shear DCMIP2.2'
    CASE DEFAULT
-      IF (is_mpi_root) PRINT *, 'Bad selector : rayleigh_friction_type =', TRIM(rayleigh_friction_key), ' in dissip_gcm.f90/init_dissip'
+      IF (is_mpi_root) PRINT *, 'Bad selector : rayleigh_friction_type =', TRIM(rayleigh_friction_key), & 
+        ' in dissip_gcm.f90/init_dissip'
       STOP
    END SELECT
Index: codes/icosagcm/trunk/src/domain.f90
===================================================================
--- codes/icosagcm/trunk/src/domain.f90	(revision 148)
+++ codes/icosagcm/trunk/src/domain.f90	(revision 149)
@@ -271,5 +271,6 @@
             nf2=domain_glo(edge_glo(e)%assign_domain)%face
             d%edge_assign_sign(k,i,j)=1-2*MOD(12+tab_index(nf,nf2,0),2)
-            IF (MOD(6+k+tab_index(nf,nf2,0),6)/=edge_glo(e)%assign_pos .AND. MOD(6+k+tab_index(nf,nf2,0),6) /= MOD(edge_glo(e)%assign_pos+3,6)) THEN
+            IF (MOD(6+k+tab_index(nf,nf2,0),6)/=edge_glo(e)%assign_pos .AND. MOD(6+k+tab_index(nf,nf2,0),6) & 
+                /= MOD(edge_glo(e)%assign_pos+3,6)) THEN
               d%edge_assign_sign(k,i,j)=-d%edge_assign_sign(k,i,j)
             ENDIF
Index: codes/icosagcm/trunk/src/dynetat0_gcm_mod.f90
===================================================================
--- codes/icosagcm/trunk/src/dynetat0_gcm_mod.f90	(revision 149)
+++ codes/icosagcm/trunk/src/dynetat0_gcm_mod.f90	(revision 149)
@@ -0,0 +1,367 @@
+MODULE dynetat0_gcm_mod 
+  USE genmod
+  USE icosa
+  USE caldyn_gcm_mod 
+ 	IMPLICIT NONE 
+	  PRIVATE
+
+       PUBLIC  etat0
+	 INTEGER,SAVE::ncell
+	 TYPE(t_field),POINTER:: f_iu(:)
+      TYPE(t_field),POINTER:: f_iv(:) 
+	 TYPE(t_field),POINTER:: f_iue(:)
+      TYPE(t_field),POINTER:: f_ive(:) 
+	 REAL(rstd),POINTER :: iu(:,:),iv(:,:)
+      REAL(rstd),POINTER :: iue(:,:),ive(:,:) 
+
+	  
+CONTAINS
+
+	SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q) 
+  USE icosa
+  USE caldyn_mod
+  USE write_field  
+  USE maxicosa
+	IMPLICIT NONE
+	TYPE(t_domain),POINTER :: d 
+	TYPE(t_field),POINTER:: f_ps(:)
+	TYPE(t_field),POINTER:: f_phis(:)
+	TYPE(t_field),POINTER:: f_u(:)
+	TYPE(t_field),POINTER:: f_q(:)
+	TYPE(t_field),POINTER:: f_theta_rhodz(:) 
+     TYPE(t_field),POINTER::  f_buf_i3(:), f_buf1_i(:), f_buf2_i(:)
+     REAL(rstd),POINTER :: ps(:)
+     REAL(rstd),POINTER :: phis(:)
+     REAL(rstd),POINTER :: theta_rhodz(:,:)
+     REAL(rstd),POINTER :: u(:,:) 
+     REAL(rstd),POINTER :: q(:,:,:)
+	REAL(rstd):: maxff,minff,maxuu,minuu
+	INTEGER :: ind
+ 
+	CALL allocate_field(f_iu,field_t,type_real,llm) 
+     CALL allocate_field(f_iv,field_t,type_real,llm)
+	CALL allocate_field(f_iue,field_u,type_real,llm) 
+     CALL allocate_field(f_ive,field_u,type_real,llm)
+	CALL allocate_field(f_u,field_u,type_real,llm) 
+	CALL allocate_field(f_buf1_i,field_t,type_real,llm) 
+	CALL allocate_field(f_buf2_i,field_t,type_real,llm) 
+	CALL allocate_field(f_buf_i3,field_u,type_real,3,llm) 
+  
+	PRINT*,"IN NETCDF READ"
+!------------------------------------zero 
+	DO ind=1,ndomain
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+	  iu = f_iu(ind) 
+	  iv = f_iv(ind) 
+      iue = f_iue(ind) 
+	 ive = f_ive(ind) 	
+	  iu = 0.0 
+	  iv = 0.0 
+        u = f_u(ind) 
+	   u = 0.0 
+	 iue = 0.0	
+	 ive = 0.0 
+     END DO 
+!--------------------------------------------
+	 ncell = 0
+     DO ind=1,ndomain
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+	 d => domain_glo(ind)
+      ps=f_ps(ind)
+      phis=f_phis(ind)
+      theta_rhodz=f_theta_rhodz(ind)
+      q=f_q(ind)
+      iu=f_iu(ind) 
+	 iv=f_iv(ind) 
+      CALL compute_dynetat0(ind,d,ps,phis,theta_rhodz,iu,iv,q)
+	ENDDO
+
+     CALL transfert_request(f_ps,req_i1)
+	CALL transfert_request(f_phis,req_i1)
+	CALL transfert_request(f_theta_rhodz,req_i1)
+	CALL transfert_request(f_q,req_i1)
+	CALL transfert_request(f_iu,req_i1)
+	CALL transfert_request(f_iv,req_i1)
+!------------------------------------------
+	DO ind=1,ndomain
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+	 u=f_u(ind)
+	 iu=f_iu(ind) 
+	 iv=f_iv(ind) 
+	 iue=f_iue(ind) 
+	 ive=f_ive(ind) 
+  	 CALL compute_dynetatu(iu,iv,iue,ive,u)
+	ENDDO
+!----------------------------------------------------
+!------------- OUTPUT OF Variables 
+	CALL un2ulonlat(f_u,f_buf_i3,f_buf1_i,f_buf2_i) 
+	CALL writefield("buf1",f_buf1_i)
+  END SUBROUTINE etat0
+
+!==================================================================
+  SUBROUTINE compute_dynetat0(ind,d,ps,phis,theta_rhodz,iu,iv,q) 
+   use icosa
+   use netcdf
+   use wind_mod 
+   USE disvert_mod
+ 	IMPLICIT NONE 
+   TYPE(t_domain),POINTER :: d 
+   CHARACTER*20::dimname 
+   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) :: q(iim*jjm,llm,nqtot)
+   REAL(rstd),ALLOCATABLE :: mass(:,:)   ! mass    
+   REAL(rstd),ALLOCATABLE :: rhodz(:,:)   ! mass density  
+   REAL(rstd),ALLOCATABLE :: theta(:,:) 
+   REAL(rstd),ALLOCATABLE :: p(:,:)  ! pression
+   REAL(rstd),POINTER :: iu(:,:),iv(:,:)
+   REAL(rstd),POINTER :: icops(:)
+   REAL(rstd),POINTER :: icophis(:)
+   REAL(rstd),POINTER :: icou(:,:),icov(:,:)
+   REAL(rstd),POINTER :: icotheta(:,:)
+   REAL(rstd),POINTER :: icoq(:,:,:)
+
+   INTEGER length,iq,ind,l
+   PARAMETER (length = 100)
+   REAL tab_cntrl(length) ! tableau des parametres du run
+   INTEGER::ierr,nid,ncid,nvarid,dimid,nind
+   INTEGER::ncells  
+   INTEGER::halo_size,i,j,k,ij
+   LOGICAL::single 
+   INTEGER::nDims,nVars,nGlobalAtts,unlimdimid
+   INTEGER:: len
+   CHARACTER(LEN=200):: iqq 
+   
+!	OPEN NETCDF FILE 
+	 ierr = NF90_OPEN ("start_icosa25.nc",NF90_NOWRITE,nid)
+      IF (ierr .NE. NF90_NOERR) THEN
+        write(*,*)'dynetat0: with file start_icosa.nc'
+        write(*,*)' ierr = ', ierr
+        STOP 
+      ENDIF
+
+     ierr= nf90_inquire(nid,nDims,nVars,nGlobalAtts,unlimdimid)
+	IF (ierr .NE. NF90_NOERR) THEN
+        write(*,*)'Problem in inquire'
+        write(*,*)' ierr = ', ierr
+        STOP 
+      ENDIF
+
+!	PRINT*,"nDims,nVars,nGlobalAtts,unlimdimid"
+!	PRINT*,nDims,nVars,nGlobalAtts,unlimdimid
+
+
+	 ierr = NF90_INQ_DIMID(nid,"ncells",dimid)
+	 IF (ierr .NE. NF90_NOERR ) THEN 
+	   write(*,*)'ncells is not present in start_icosa.nc'
+        write(*,*)' ierr = ', ierr
+	   STOP 
+	 ENDIF 
+	   
+	 ierr = nf90_inquire_dimension(nid,dimid,dimname,ncells)
+	 IF (ierr .NE. NF90_NOERR ) THEN 
+	   write(*,*)'ncells  in start_icosa.nc'
+        write(*,*)' ierr = ', ierr
+	   STOP 
+	 ENDIF
+
+	  ALLOCATE(icops(ncells))
+	  ALLOCATE(icophis(ncells))
+	  ALLOCATE(icou(ncells,llm))
+	  ALLOCATE(icov(ncells,llm))
+	  ALLOCATE(icotheta(ncells,llm))
+	  ALLOCATE(icoq(ncells,llm,nqtot))
+	  ALLOCATE(p(iim*jjm,llm+1))
+	  ALLOCATE(theta(iim*jjm,llm)) 
+       ALLOCATE(mass(iim*jjm,llm))   ! mass    
+       ALLOCATE(rhodz(iim*jjm,llm))   ! mass density   
+!============================================================
+      ierr = NF90_INQ_VARID(nid, "phisinit", nvarid)
+      IF (ierr .NE. NF90_NOERR) THEN
+        write(*,*)"dynetat0: phisinit is absent"
+	   write(*,*)' ierr = ', ierr
+         STOP       
+      ENDIF
+
+      ierr = NF90_GET_VAR(nid, nvarid, icophis)
+      IF (ierr .NE. NF90_NOERR) THEN
+         write(*,*)"dynetat0: PROBLEM IN PHIS"
+         STOP
+      ENDIF
+!==============================================================
+	  ierr = NF90_INQ_VARID(nid, "ps", nvarid)
+      IF (ierr .NE. NF90_NOERR) THEN
+        write(*,*)"dynetat0: ps is absent"
+	   write(*,*)' ierr = ', ierr
+         STOP       
+      ENDIF
+
+      ierr = NF90_GET_VAR(nid, nvarid, icops)
+      IF (ierr .NE. NF90_NOERR) THEN
+         write(*,*)"dynetat0: PROBLEM IN PS"
+         STOP
+      ENDIF
+!================================================================
+	  ierr = NF90_INQ_VARID(nid, "teta", nvarid)
+      IF (ierr .NE. NF90_NOERR) THEN
+        write(*,*)"dynetat0: teta is not available in start.nc"
+	   write(*,*)' ierr = ', ierr
+         STOP       
+      ENDIF
+
+      ierr = NF90_GET_VAR(nid, nvarid,icotheta)
+      IF (ierr .NE. NF90_NOERR) THEN
+         write(*,*)"dynetat0: PROBLEM IN Teta"
+         STOP
+      ENDIF
+!================================================================
+	DO iq = 1,nqtot   
+		write(iqq,*)INT(iq)
+		iqq=ADJUSTL(iqq)  
+     	ierr = NF90_INQ_VARID(nid,"q"//iqq, nvarid)
+      IF (ierr .NE. NF90_NOERR) THEN
+        write(*,*)"dynetat0: ","q"//iqq,"not here"
+	   write(*,*)' ierr = ', ierr
+!         STOP       
+      ENDIF
+
+      ierr = NF90_GET_VAR(nid, nvarid,icoq(:,:,iq))
+      IF (ierr .NE. NF90_NOERR) THEN
+         write(*,*)"dynetat0: PROBLEM IN Q"
+!         STOP
+      ENDIF
+	END DO
+!================================================================
+	GO TO 121
+	  ierr = NF90_INQ_VARID(nid, "q01", nvarid)
+      IF (ierr .NE. NF90_NOERR) THEN
+        write(*,*)"dynetat0: q1 is not available in start.nc"
+	   write(*,*)' ierr = ', ierr
+         STOP       
+      ENDIF
+
+      ierr = NF90_GET_VAR(nid, nvarid,icoq(:,:,1))
+      IF (ierr .NE. NF90_NOERR) THEN
+         write(*,*)"dynetat0: PROBLEM IN Q01"
+         STOP
+      ENDIF
+121	CONTINUE 
+!================================================================
+	  ierr = NF90_INQ_VARID(nid, "ucov", nvarid)
+      IF (ierr .NE. NF90_NOERR) THEN
+        write(*,*)"dynetat0: ucov is not available in start.nc"
+	   write(*,*)' ierr = ', ierr
+         STOP       
+      ENDIF
+
+      ierr = NF90_GET_VAR(nid, nvarid,icou)
+      IF (ierr .NE. NF90_NOERR) THEN
+         write(*,*)"dynetat0: PROBLEM IN ucov"
+         STOP
+      ENDIF
+	PRINT*,"UCOV is read using start_icosa.nc" 
+!================================================================
+	  ierr = NF90_INQ_VARID(nid, "vcov", nvarid)
+      IF (ierr .NE. NF90_NOERR) THEN
+        write(*,*)"dynetat0: PROBLEM in VCOV"
+	   write(*,*)' ierr = ', ierr
+         STOP       
+      ENDIF
+
+      ierr = NF90_GET_VAR(nid, nvarid,icov)
+      IF (ierr .NE. NF90_NOERR) THEN
+         write(*,*)"dynetat0: PROBLEM IN vcov"
+         STOP
+      ENDIF
+!================================================================
+			iu = 0.0 ; iv = 0.0 
+             DO j=d%jj_begin,d%jj_end
+              DO i=d%ii_begin,d%ii_end
+                  k=d%iim*(j-1)+i
+                 IF (d%assign_domain(i,j)==ind ) THEN
+                     ncell=ncell+1
+		           phis(k)=  0.0  !icophis(ncell) 
+	                ps(k)= icops(ncell) 
+	                theta(k,:) = icotheta(ncell,:) 
+	                q(k,:,1)= icoq(ncell,:,1) 
+	                iu(k,:) = icou(ncell,:)
+				 iv(k,:) = icov(ncell,:) 
+	         ENDIF
+                ENDDO
+              ENDDO
+
+    DO    l    = 1, llm+1
+      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
+
+   DO l = 1, llm
+     DO j=jj_begin,jj_end
+       DO i=ii_begin,ii_end
+         ij=(j-1)*iim+i
+         mass(ij,l) = ( p(ij,l) - p(ij,l+1) )*Ai(ij)/g
+         rhodz(ij,l) = mass(ij,l) / Ai(ij)
+       ENDDO
+     ENDDO
+   ENDDO
+
+    DO    l    = 1, llm
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          theta_rhodz(ij,l) = theta(ij,l)*rhodz(ij,l)
+        ENDDO
+      ENDDO
+    ENDDO
+
+	  DEALLOCATE(icops)
+	  DEALLOCATE(icophis)
+	  DEALLOCATE(icou)
+	  DEALLOCATE(icov)
+	  DEALLOCATE(icotheta)
+	  DEALLOCATE(p)
+	  DEALLOCATE(theta) 
+       DEALLOCATE(mass)   ! mass    
+       DEALLOCATE(rhodz)   !
+	END SUBROUTINE compute_dynetat0
+
+!==================================================================
+	  SUBROUTINE compute_dynetatu(iu,iv,iue,ive,u) 
+   use icosa
+   use wind_mod 
+    	IMPLICIT NONE 
+   CHARACTER*20::dimname 
+   REAL(rstd),INTENT(OUT):: u(3*iim*jjm,llm)
+   REAL(rstd) :: iu(iim*jjm,llm),iv(iim*jjm,llm)
+   REAL(rstd) :: iue(3*iim*jjm,llm),ive(3*iim*jjm,llm)
+   INTEGER::halo_size,i,j,k,ij,l
+
+
+  Do l = 1, llm
+   DO j=jj_begin-1,jj_end+1
+      DO i=ii_begin-1,ii_end+1
+         k=iim*(j-1)+i
+	   iue(k+u_right,l)=0.5*(iu(k,l)+iu(k+t_right,l))
+        iue(k+u_lup,l)  =0.5*(iu(k,l)+iu(k+t_lup,l))
+        iue(k+u_ldown,l)=0.5*(iu(k,l)+iu(k+t_ldown,l)) 
+!------------------------------------------------------ 
+	   ive(k+u_right,l)=0.5*(iv(k,l)+iv(k+t_right,l))
+        ive(k+u_lup,l)  =0.5*(iv(k,l)+iv(k+t_lup,l))
+        ive(k+u_ldown,l)=0.5*(iv(k,l)+iv(k+t_ldown,l)) 
+	 END DO
+    END DO 
+  END DO 
+	CALL compute_wind_perp_from_lonlat_compound(iue,ive,u) 
+
+	END SUBROUTINE compute_dynetatu
+
+
+ END MODULE dynetat0_gcm_mod 
Index: codes/icosagcm/trunk/src/dynetat0_hz_mod.f90
===================================================================
--- codes/icosagcm/trunk/src/dynetat0_hz_mod.f90	(revision 149)
+++ codes/icosagcm/trunk/src/dynetat0_hz_mod.f90	(revision 149)
@@ -0,0 +1,348 @@
+MODULE dynetat0_hz_mod 
+  USE genmod
+  USE icosa
+  USE caldyn_gcm_mod 
+  IMPLICIT NONE 
+  PRIVATE
+
+       PUBLIC  etat0
+       INTEGER,SAVE::ncell
+       TYPE(t_field),POINTER:: f_iu(:)
+       TYPE(t_field),POINTER:: f_iv(:) 
+       TYPE(t_field),POINTER:: f_iue(:)
+       TYPE(t_field),POINTER:: f_ive(:) 
+       REAL(rstd),POINTER :: iu(:,:),iv(:,:)
+       REAL(rstd),POINTER :: iue(:,:),ive(:,:) 
+
+ 
+CONTAINS
+
+    SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q) 
+  USE icosa
+  USE caldyn_mod
+  USE write_field  
+  USE maxicosa
+        IMPLICIT NONE
+     TYPE(t_domain),POINTER :: d 
+     TYPE(t_field),POINTER:: f_ps(:)
+     TYPE(t_field),POINTER:: f_phis(:)
+     TYPE(t_field),POINTER:: f_u(:)
+     TYPE(t_field),POINTER:: f_q(:)
+     TYPE(t_field),POINTER:: f_theta_rhodz(:) 
+     TYPE(t_field),POINTER::  f_buf_i3(:), f_buf1_i(:), f_buf2_i(:)
+     REAL(rstd),POINTER :: ps(:)
+     REAL(rstd),POINTER :: phis(:)
+     REAL(rstd),POINTER :: theta_rhodz(:,:)
+     REAL(rstd),POINTER :: u(:,:) 
+     REAL(rstd),POINTER :: q(:,:,:)
+     REAL(rstd):: maxff,minff,maxuu,minuu
+     INTEGER :: ind
+ 
+     CALL allocate_field(f_iu,field_t,type_real,llm) 
+     CALL allocate_field(f_iv,field_t,type_real,llm)
+     CALL allocate_field(f_iue,field_u,type_real,llm) 
+     CALL allocate_field(f_ive,field_u,type_real,llm)
+     CALL allocate_field(f_u,field_u,type_real,llm) 
+     CALL allocate_field(f_buf1_i,field_t,type_real,llm) 
+     CALL allocate_field(f_buf2_i,field_t,type_real,llm) 
+     CALL allocate_field(f_buf_i3,field_u,type_real,3,llm) 
+  
+     PRINT*,"IN NETCDF READ"
+!------------------------------------zero 
+	DO ind=1,ndomain
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+	  iu = f_iu(ind) 
+	  iv = f_iv(ind) 
+      iue = f_iue(ind) 
+	 ive = f_ive(ind) 	
+	   u = f_u(ind) 
+	  iu = 0.0 
+	  iv = 0.0 
+	   u = 0.0 
+	 iue = 0.0	
+	 ive = 0.0 
+     END DO 
+!--------------------------------------------
+     ncell = 0
+     DO ind=1,ndomain
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      d => domain_glo(ind)
+      ps=f_ps(ind)
+      phis=f_phis(ind)
+      theta_rhodz=f_theta_rhodz(ind)
+      q=f_q(ind)
+      iu=f_iu(ind) 
+      iv=f_iv(ind) 
+      CALL compute_dynetat0(ind,d,ps,phis,theta_rhodz,iu,iv,q)
+      ENDDO
+
+     CALL transfert_request(f_ps,req_i1)
+	CALL transfert_request(f_phis,req_i1)
+	CALL transfert_request(f_theta_rhodz,req_i1)
+	CALL transfert_request(f_q,req_i1)
+	CALL transfert_request(f_iu,req_i1)
+	CALL transfert_request(f_iv,req_i1)
+!------------------------------------------
+	DO ind=1,ndomain
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+	 u=f_u(ind)
+	 iu=f_iu(ind) 
+	 iv=f_iv(ind) 
+	 iue=f_iue(ind) 
+	 ive=f_ive(ind) 
+  	 CALL compute_dynetatu(iu,iv,iue,ive,u)
+	ENDDO
+!------------- OUTPUT OF Variables 
+	CALL un2ulonlat(f_u,f_buf_i3,f_buf1_i,f_buf2_i) 
+	CALL writefield("buf1",f_buf1_i)
+  END SUBROUTINE etat0
+
+!==================================================================
+  SUBROUTINE compute_dynetat0(ind,d,ps,phis,theta_rhodz,iu,iv,q) 
+   use icosa
+   use netcdf
+   use wind_mod 
+   USE disvert_mod
+ 	IMPLICIT NONE 
+   TYPE(t_domain),POINTER :: d 
+   CHARACTER*20::dimname 
+   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) :: q(iim*jjm,llm,nqtot)
+   REAL(rstd),ALLOCATABLE :: mass(:,:)   ! mass    
+   REAL(rstd),ALLOCATABLE :: rhodz(:,:)   ! mass density  
+   REAL(rstd),ALLOCATABLE :: theta(:,:) 
+   REAL(rstd),ALLOCATABLE :: p(:,:)  ! pression
+   REAL(rstd),POINTER :: iu(:,:),iv(:,:)
+   REAL(rstd),POINTER :: icops(:)
+   REAL(rstd),POINTER :: icophis(:)
+   REAL(rstd),POINTER :: icou(:,:),icov(:,:)
+   REAL(rstd),POINTER :: icotheta(:,:)
+   REAL(rstd),POINTER :: icoq(:,:,:)
+   
+   INTEGER length,iq,ind,l
+   PARAMETER (length = 1)
+   CHARACTER(LEN=200):: iqq ! tableau des parametres du run
+   INTEGER::ierr,nid,ncid,nvarid,dimid,nind
+   INTEGER::ncells  
+   INTEGER::halo_size,i,j,k,ij
+   LOGICAL::single 
+   INTEGER::nDims,nVars,nGlobalAtts,unlimdimid
+   INTEGER:: len
+   
+!	OPEN NETCDF FILE 
+	 ierr = NF90_OPEN ("icosa_hz30.nc",NF90_NOWRITE,nid)
+      IF (ierr .NE. NF90_NOERR) THEN
+        write(*,*)'dynetat0: with file icosa_hz30.nc'
+        write(*,*)' ierr = ', ierr
+        STOP 
+      ENDIF
+
+     ierr= nf90_inquire(nid,nDims,nVars,nGlobalAtts,unlimdimid)
+	IF (ierr .NE. NF90_NOERR) THEN
+        write(*,*)'Problem in inquire'
+        write(*,*)' ierr = ', ierr
+        STOP 
+      ENDIF
+
+	PRINT*,"nDims,nVars,nGlobalAtts,unlimdimid"
+	PRINT*,nDims,nVars,nGlobalAtts,unlimdimid
+
+
+	 ierr = NF90_INQ_DIMID(nid,"ncells",dimid)
+	 IF (ierr .NE. NF90_NOERR ) THEN 
+	   write(*,*)'ncells is not present in hzdy_30.nc'
+        write(*,*)' ierr = ', ierr
+	   STOP 
+	 ENDIF 
+	   
+	 ierr = nf90_inquire_dimension(nid,dimid,dimname,ncells)
+	 IF (ierr .NE. NF90_NOERR ) THEN 
+	   write(*,*)'ncells  in hzdy_30.nc'
+        write(*,*)' ierr = ', ierr
+	   STOP 
+	 ENDIF
+
+	  ALLOCATE(icops(ncells))
+	  ALLOCATE(icophis(ncells))
+	  ALLOCATE(icou(ncells,llm))
+	  ALLOCATE(icov(ncells,llm))
+	  ALLOCATE(icotheta(ncells,llm))
+	  ALLOCATE(icoq(ncells,llm,nqtot))
+       ALLOCATE(p(iim*jjm,llm+1))
+	  ALLOCATE(theta(iim*jjm,llm)) 
+       ALLOCATE(mass(iim*jjm,llm))   ! mass    
+       ALLOCATE(rhodz(iim*jjm,llm))   ! mass density   
+!============================================================
+      ierr = NF90_INQ_VARID(nid, "phis", nvarid)
+      IF (ierr .NE. NF90_NOERR) THEN
+        write(*,*)"dynetat0: phis is absent"
+	   write(*,*)' ierr = ', ierr
+         STOP       
+      ENDIF
+
+      ierr = NF90_GET_VAR(nid, nvarid, icophis)
+      IF (ierr .NE. NF90_NOERR) THEN
+         write(*,*)"dynetat0: PROBLEM IN PHIS"
+         STOP
+      ENDIF
+!==============================================================
+	  ierr = NF90_INQ_VARID(nid, "ps", nvarid)
+      IF (ierr .NE. NF90_NOERR) THEN
+        write(*,*)"dynetat0: ps is absent"
+	   write(*,*)' ierr = ', ierr
+         STOP       
+      ENDIF
+
+      ierr = NF90_GET_VAR(nid, nvarid, icops)
+      IF (ierr .NE. NF90_NOERR) THEN
+         write(*,*)"dynetat0: PROBLEM IN PS"
+         STOP
+      ENDIF
+!================================================================
+	  ierr = NF90_INQ_VARID(nid, "theta", nvarid)
+      IF (ierr .NE. NF90_NOERR) THEN
+        write(*,*)"dynetat0: teta is not available in start.nc"
+	   write(*,*)' ierr = ', ierr
+         STOP       
+      ENDIF
+
+      ierr = NF90_GET_VAR(nid, nvarid,icotheta)
+      IF (ierr .NE. NF90_NOERR) THEN
+         write(*,*)"dynetat0: PROBLEM IN Teta"
+         STOP
+      ENDIF
+!================================================================
+	DO iq = 1,nqtot   
+		write(iqq,*)INT(iq)
+		iqq=ADJUSTL(iqq)  
+     	ierr = NF90_INQ_VARID(nid,"q"//iqq, nvarid)
+      IF (ierr .NE. NF90_NOERR) THEN
+        write(*,*)"dynetat0: ","q"//iqq,"not here"
+	   write(*,*)' ierr = ', ierr
+         STOP       
+      ENDIF
+
+      ierr = NF90_GET_VAR(nid, nvarid,icoq(:,:,iq))
+      IF (ierr .NE. NF90_NOERR) THEN
+         write(*,*)"dynetat0: PROBLEM IN Q"
+         STOP
+      ENDIF
+	END DO
+!===============================================================
+	 ierr = NF90_INQ_VARID(nid, "ulon", nvarid)
+      IF (ierr .NE. NF90_NOERR) THEN
+        write(*,*)"dynetat0: ulon is not in file"
+	   write(*,*)' ierr = ', ierr
+         STOP       
+      ENDIF
+	
+      ierr = NF90_GET_VAR(nid, nvarid,icou)
+      IF (ierr .NE. NF90_NOERR) THEN
+         write(*,*)"dynetat0: PROBLEM IN ucov"
+         STOP
+      ENDIF
+!================================================================
+	  ierr = NF90_INQ_VARID(nid, "ulat", nvarid)
+      IF (ierr .NE. NF90_NOERR) THEN
+        write(*,*)"dynetat0: ulat is not available in start.nc"
+	   write(*,*)' ierr = ', ierr
+         STOP       
+      ENDIF
+
+      ierr = NF90_GET_VAR(nid, nvarid,icov)
+      IF (ierr .NE. NF90_NOERR) THEN
+         write(*,*)"dynetat0: PROBLEM IN vlat"
+         STOP
+      ENDIF
+!================================================================
+			iu = 0.0 ; iv = 0.0 
+             DO j=d%jj_begin,d%jj_end
+              DO i=d%ii_begin,d%ii_end
+                  k=d%iim*(j-1)+i
+                 IF (d%assign_domain(i,j)==ind ) THEN
+                     ncell=ncell+1
+				 phis(k)= icophis(ncell) 
+	                ps(k)= icops(ncell) 
+	                theta(k,:) = icotheta(ncell,:) 
+	                q(k,:,:)= icoq(ncell,:,:) 
+	                iu(k,:) = icou(ncell,:)
+				 iv(k,:) = icov(ncell,:)
+			   ENDIF
+                ENDDO
+              ENDDO
+
+    DO    l    = 1, llm+1
+      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
+
+   DO l = 1, llm
+     DO j=jj_begin,jj_end
+       DO i=ii_begin,ii_end
+         ij=(j-1)*iim+i
+         mass(ij,l) = ( p(ij,l) - p(ij,l+1) )*Ai(ij)/g
+         rhodz(ij,l) = mass(ij,l) / Ai(ij)
+       ENDDO
+     ENDDO
+   ENDDO
+
+    DO    l    = 1, llm
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          theta_rhodz(ij,l) = theta(ij,l)*rhodz(ij,l)
+        ENDDO
+      ENDDO
+    ENDDO
+
+	  DEALLOCATE(icops)
+	  DEALLOCATE(icophis)
+	  DEALLOCATE(icou)
+	  DEALLOCATE(icov)
+	  DEALLOCATE(icotheta)
+	  DEALLOCATE(p)
+	  DEALLOCATE(theta) 
+       DEALLOCATE(mass)   ! mass    
+       DEALLOCATE(rhodz)   !
+	END SUBROUTINE compute_dynetat0
+!==================================================================
+	  SUBROUTINE compute_dynetatu(iu,iv,iue,ive,u) 
+   use icosa
+   use wind_mod 
+    	IMPLICIT NONE 
+   CHARACTER*20::dimname 
+   REAL(rstd),INTENT(OUT):: u(3*iim*jjm,llm)
+   REAL(rstd) :: iu(iim*jjm,llm),iv(iim*jjm,llm)
+   REAL(rstd) :: iue(3*iim*jjm,llm),ive(3*iim*jjm,llm)
+   INTEGER::halo_size,i,j,k,ij,l
+
+
+  Do l = 1, llm
+   DO j=jj_begin-1,jj_end+1
+      DO i=ii_begin-1,ii_end+1
+         k=iim*(j-1)+i
+	   iue(k+u_right,l)=0.5*(iu(k,l)+iu(k+t_right,l))
+        iue(k+u_lup,l)  =0.5*(iu(k,l)+iu(k+t_lup,l))
+        iue(k+u_ldown,l)=0.5*(iu(k,l)+iu(k+t_ldown,l)) 
+!------------------------------------------------------ 
+	   ive(k+u_right,l)=0.5*(iv(k,l)+iv(k+t_right,l))
+        ive(k+u_lup,l)  =0.5*(iv(k,l)+iv(k+t_lup,l))
+        ive(k+u_ldown,l)=0.5*(iv(k,l)+iv(k+t_ldown,l)) 
+	 END DO
+    END DO 
+  END DO 
+	CALL compute_wind_perp_from_lonlat_compound(iue,ive,u) 
+
+	END SUBROUTINE compute_dynetatu
+
+
+ END MODULE dynetat0_hz_mod 
Index: codes/icosagcm/trunk/src/etat0.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0.f90	(revision 148)
+++ codes/icosagcm/trunk/src/etat0.f90	(revision 149)
@@ -1,3 +1,4 @@
 MODULE etat0_mod
+    CHARACTER(len=255),SAVE :: etat0_type
 
 CONTAINS
@@ -12,4 +13,8 @@
     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 
+    USE dynetat0_hz_mod,  ONLY : dynetat0_hz=>etat0 
+
     IMPLICIT NONE
     TYPE(t_field),POINTER :: f_ps(:)
@@ -19,5 +24,4 @@
     TYPE(t_field),POINTER :: f_q(:)
     
-    CHARACTER(len=255) :: etat0_type
     etat0_type='jablonowsky06'
     CALL getin("etat0",etat0_type)
@@ -28,4 +32,7 @@
     CASE ('academic')
        CALL etat0_academic(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
+    CASE ('heldsz')
+	print*,"heldsz test case"
+       CALL etat0_heldsz(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
     CASE ('dcmip1')
        CALL etat0_dcmip1(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
@@ -38,4 +45,10 @@
      CASE ('dcmip5')
        CALL etat0_dcmip5(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
+     CASE ('readnf_start') 
+	  print*,"readnf_start used"	
+       CALL dynetat0_start(f_ps,f_phis,f_theta_rhodz,f_u,f_q) 
+	CASE ('readnf_hz') 
+	  print*,"readnf_hz used"
+       CALL dynetat0_hz(f_ps,f_phis,f_theta_rhodz,f_u,f_q) 
    CASE DEFAULT
        PRINT*, 'Bad selector for variable etat0 <',etat0_type, &
Index: codes/icosagcm/trunk/src/etat0_heldsz.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0_heldsz.f90	(revision 149)
+++ codes/icosagcm/trunk/src/etat0_heldsz.f90	(revision 149)
@@ -0,0 +1,250 @@
+    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(:)
+
+CONTAINS
+  
+  SUBROUTINE test_etat0_heldsz
+  USE icosa
+  USE kinetic_mod
+  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(:)
+    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)
+    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_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)
+
+    CALL writefield('ps',f_ps)
+    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
+    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
+    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)) 
+
+    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)
+    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
+
+    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
+             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)
+
+  END SUBROUTINE compute_etat0_heldsz
+
+
+	SUBROUTINE held_saurez(f_ps,f_theta_rhodz,f_u) 
+	USE icosa
+	    IMPLICIT NONE
+    TYPE(t_field),POINTER :: f_theta_rhodz(:)
+    TYPE(t_field),POINTER :: f_u(:)
+    TYPE(t_field),POINTER :: f_ps(:)
+    REAL(rstd),POINTER :: theta_rhodz(:,:)
+    REAL(rstd),POINTER :: u(:,:)
+    REAL(rstd),POINTER :: ps(:)
+    REAL(rstd),POINTER :: tetarappel(:,:)
+    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)
+    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
+
+END MODULE etat0_heldsz_mod
Index: codes/icosagcm/trunk/src/icosa_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/icosa_gcm.f90	(revision 148)
+++ codes/icosagcm/trunk/src/icosa_gcm.f90	(revision 149)
@@ -25,4 +25,5 @@
   CALL init_transfert
   CALL init_writefield
+
 !  CALL allocate_field(sum_ne,field_T,type_real)
 !  CALL allocate_field_glo(sum_ne_glo,field_T,type_real)
Index: codes/icosagcm/trunk/src/maxico.f90
===================================================================
--- codes/icosagcm/trunk/src/maxico.f90	(revision 149)
+++ codes/icosagcm/trunk/src/maxico.f90	(revision 149)
@@ -0,0 +1,218 @@
+ MODULE maxicosa
+contains
+    SUBROUTINE maxis(f_fl,maxff,minff) 
+     USE icosa
+     IMPLICIT NONE
+ INTEGER::i,j,ij,l,ind
+ TYPE(t_field),POINTER::f_fl(:) 
+ REAL(rstd),POINTER::fl(:)  
+ REAL(rstd)::maxf,minf 
+ REAL(rstd),INTENT(OUT)::maxff,minff 
+ 
+    maxff=-1e50
+    minff=1e50 
+    DO ind=1,ndomain 
+     CALL swap_dimensions(ind)
+     CALL swap_geometry(ind) 
+     fl=f_fl(ind) 
+     CALL compute_maxis(fl,maxf,minf)  
+     maxff=max(maxff,maxf) 
+     minff=min(minff,minf) 
+    ENDDO 
+   END SUBROUTINE maxis
+!------------------------------------------------------------
+	    SUBROUTINE maxiv(f_fu,maxuu,minuu) 
+     USE icosa
+     IMPLICIT NONE
+ INTEGER::i,j,ij,l,ind
+ TYPE(t_field),POINTER::f_fu(:) 
+ REAL(rstd),POINTER::fu(:,:)  
+ REAL(rstd)::maxu,minu 
+ REAL(rstd),INTENT(OUT)::maxuu,minuu 
+ 
+    maxuu=-1e50
+    minuu=1e50 
+    DO ind=1,ndomain 
+     CALL swap_dimensions(ind)
+     CALL swap_geometry(ind) 
+     fu=f_fu(ind) 
+     CALL compute_maxiv(fu,maxu,minu)  
+	maxuu=max(maxuu,maxu) 
+	minuu=min(minuu,minu) 
+    ENDDO 
+	END SUBROUTINE maxiv
+
+    SUBROUTINE errornorms(f_ref,f_sim,norml1,norml2,normlf) 
+     USE icosa
+     IMPLICIT NONE
+ INTEGER::i,j,ij,l,ind
+ TYPE(t_field),POINTER::f_ref(:),f_sim(:) 
+ REAL(rstd),POINTER::ref(:),sim(:)  
+ REAL(rstd),INTENT(OUT)::norml1,norml2,normlf 
+ REAL(rstd)::maxf,minf
+ REAL(rstd)::maxff,minff
+ 
+    norml1 = 0.0 ; norml2 = 0.0 ; normlf = 0.0
+    DO ind=1,ndomain 
+     CALL swap_dimensions(ind)
+     CALL swap_geometry(ind) 
+      ref=f_ref(ind) 
+      sim=f_sim(ind) 
+     CALL compute_maxis(sim,maxf,minf)  
+     maxff=max(maxff,maxf) 
+     minff=min(minff,minf) 
+    ENDDO 
+	END SUBROUTINE errornorms
+
+	SUBROUTINE compute_maxis(fl,maxf,minf) 
+	USE icosa
+	IMPLICIT NONE 
+	REAL(rstd),INTENT(IN)::fl(iim*jjm)
+	REAL(rstd),INTENT(OUT)::maxf,minf 
+	INTEGER::i,j,ij 
+
+	maxf=-1e50
+	minf= 1e50 
+     DO j=jj_begin,jj_end
+       DO i=ii_begin,ii_end
+           ij=(j-1)*iim+i
+          IF ( maxf .LT. fl(ij) ) Then 
+	      maxf = fl(ij) 
+	     ENDIF 
+	    IF ( minf .GT. fl(ij) ) Then 
+	     minf = fl(ij) 
+	    ENDIF 
+	   ENDDO
+      ENDDO
+	END SUBROUTINE compute_maxis
+!------------------------------------------------------
+		SUBROUTINE compute_maxiv(fu,maxu,minu) 
+	USE icosa
+	IMPLICIT NONE 
+	REAL(rstd),INTENT(IN)::fu(3*iim*jjm,llm) 
+	REAL(rstd),INTENT(OUT)::maxu,minu  
+	INTEGER::i,j,ij 
+
+     maxu=-1e50
+     minu=1e50 
+	DO j=jj_begin-1,jj_end+1
+         DO i=ii_begin-1,ii_end+1
+           ij=(j-1)*iim+i
+          IF ( maxu .LT. fu(ij+u_right,llm) ) Then 
+	      maxu = fu(ij+u_right,llm) 
+	     ENDIF 
+		IF ( maxu .LT. fu(ij+u_lup,llm) ) Then 
+	      maxu = fu(ij+u_lup,llm) 
+	     ENDIF
+	     IF ( maxu .LT. fu(ij+u_ldown,llm) ) Then 
+	      maxu = fu(ij+u_ldown,llm) 
+	     ENDIF
+!------------------------------
+		IF ( minu .GT. fu(ij+u_right,llm) ) Then 
+	      minu = fu(ij+u_right,llm) 
+	     ENDIF 
+		IF ( minu .GT. fu(ij+u_lup,llm) ) Then 
+	      minu = fu(ij+u_lup,llm) 
+	     ENDIF
+	     IF ( maxu .GT. fu(ij+u_ldown,llm) ) Then 
+	      minu = fu(ij+u_ldown,llm) 
+	     ENDIF 
+	   ENDDO
+      ENDDO
+	END SUBROUTINE compute_maxiv
+!--------------------------------------------------
+	Subroutine compute_l1(reference,simulated,nrml1)
+	use icosa
+	IMPLICIT NONE 
+	REAL(rstd),INTENT(IN)::reference(iim*jjm) 
+	REAL(rstd),INTENT(IN)::simulated(iim*jjm) 
+	REAL(rstd),INTENT(OUT)::nrml1
+	REAL(rstd) :: temp1,temp2
+	INTEGER::i,j,ij 
+
+	temp1 = 0.0 ; temp2 = 0.0 
+	nrml1 = 0.0 
+	DO j=jj_begin,jj_end
+         DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+	  temp1 = temp1 + abs(simulated(ij) - reference(ij))*Ai(ij) 
+	  temp2 = temp2 + abs(reference(ij))*Ai(ij) 
+	 END DO
+	END DO
+	  nrml1 = temp1/temp2 
+	END SUBROUTINE compute_l1 
+
+	Subroutine compute_l2(reference,simulated,nrml2)
+	use icosa
+	IMPLICIT NONE 
+	REAL(rstd),INTENT(IN)::reference(iim*jjm) 
+	REAL(rstd),INTENT(IN)::simulated(iim*jjm) 
+	REAL(rstd),INTENT(OUT)::nrml2
+	REAL(rstd) :: temp1,temp2
+	INTEGER::i,j,ij 
+
+	temp1 = 0.0 ; temp2 = 0.0 
+	nrml2 = 0.0 
+	DO j=jj_begin,jj_end
+         DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+	  temp1 = temp1 + (simulated(ij) - reference(ij))*(simulated(ij) - reference(ij))*Ai(ij) 
+	  temp2 = temp2 + reference(ij)*reference(ij)*Ai(ij) 
+	 END DO
+	END DO
+	  temp1 = sqrt(temp1) 
+	  temp2 = sqrt(temp2)  
+	  nrml2 = temp1/temp2 
+	END SUBROUTINE compute_l2 
+
+	Subroutine compute_lf(reference,simulated,nrmlf)
+	use icosa
+	IMPLICIT NONE 
+	REAL(rstd),INTENT(IN)::reference(iim*jjm) 
+	REAL(rstd),INTENT(IN)::simulated(iim*jjm) 
+	REAL(rstd),INTENT(OUT)::nrmlf
+	REAL(rstd):: difference(iim*jjm) 
+	REAL(rstd) :: temp1,temp2
+	INTEGER::i,j,ij 
+
+	temp1 = 0.0 ; temp2 = 0.0 
+	nrmlf = 0.0 
+	
+	temp1 = maxval(ABS(difference))
+	temp2 = maxval(ABS(reference)) 
+	  nrmlf = temp1/temp2 
+	END SUBROUTINE compute_lf 
+!================================================== initial0 
+	 SUBROUTINE initial0(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
+      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.0
+      ps=0.0
+	 phis=0.0
+      theta_rhodz = 0.0
+	 u = 0.0
+    ENDDO
+  END SUBROUTINE initial0
+
+	END MODULE maxicosa
Index: codes/icosagcm/trunk/src/phyparam.F
===================================================================
--- codes/icosagcm/trunk/src/phyparam.F	(revision 149)
+++ codes/icosagcm/trunk/src/phyparam.F	(revision 149)
@@ -0,0 +1,601 @@
+	MODULE PHY
+         USE dimphys_mod
+         LOGICAL:: firstcall,lastcall
+	contains      
+
+	SUBROUTINE phyparam_lmd(it,ngrid,nlayer,nq,
+     s            ptimestep,lati,long,rjourvrai,gmtime,
+     s            pplev,pplay,pphi,pphis,
+     s            pu,pv,pt,pq,
+     s            pdu,pdv,pdt,pdq,pdpsrf)
+
+	USE ICOSA
+	USE dimphys_mod
+	USE RADIATION
+	USE SURFACE_PROCESS
+c
+      IMPLICIT NONE
+c=======================================================================
+c
+c   subject:
+c   --------
+c
+c   Organisation of the physical parametrisations of the LMD 
+c   20 parameters GCM for planetary atmospheres.
+c   It includes:
+c   raditive transfer (long and shortwave) for CO2 and dust.
+c   vertical turbulent mixing
+c   convective adjsutment
+c
+c   author: Frederic Hourdin 15 / 10 /93
+c   -------
+c
+c   arguments:
+c   ----------
+c
+c   input:
+c   ------
+c
+c    ngrid                 Size of the horizontal grid.
+c                          All internal loops are performed on that grid.
+c    nlayer                Number of vertical layers.
+c    nq                    Number of advected fields
+c    firstcall             True at the first call
+c    lastcall              True at the last call
+c    rjourvrai                  Number of days counted from the North. Spring
+c                          equinoxe.
+c    gmtime                 hour (s)
+c    ptimestep             timestep (s)
+c    pplay(ngrid,nlayer+1) Pressure at the middle of the layers (Pa)
+c    pplev(ngrid,nlayer+1) intermediate pressure levels (pa)
+c    pphi(ngrid,nlayer)    Geopotential at the middle of the layers (m2s-2)
+c    pu(ngrid,nlayer)      u component of the wind (ms-1)
+c    pv(ngrid,nlayer)      v component of the wind (ms-1)
+c    pt(ngrid,nlayer)      Temperature (K)
+c    pq(ngrid,nlayer,nq)   Advected fields
+c    pudyn(ngrid,nlayer)    \ 
+c    pvdyn(ngrid,nlayer)     \ Dynamical temporal derivative for the
+c    ptdyn(ngrid,nlayer)     / corresponding variables
+c    pqdyn(ngrid,nlayer,nq) /
+c    pw(ngrid,?)           vertical velocity
+c
+c   output:
+c   -------
+c
+c    pdu(ngrid,nlayer)        \
+c    pdv(ngrid,nlayer)         \  Temporal derivative of the corresponding
+c    pdt(ngrid,nlayer)         /  variables due to physical processes.
+c    pdq(ngrid,nlayer)      /
+c    pdpsrf(ngrid)        /
+c
+c=======================================================================
+c
+!c-----------------------------------------------------------------------
+!c
+!c    0.  Declarations :
+!c    ------------------
+!c    Arguments :
+!c    -----------
+
+!c    inputs:
+!c    -------
+	INTEGER ngrid,nlayer,nq,it,ij,i
+      REAL ptimestep
+      real zdtime
+      REAL pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
+      REAL pphi(ngrid,nlayer)
+      REAL pphis(ngrid)
+      REAL pu(ngrid,nlayer),pv(ngrid,nlayer)
+      REAL pt(ngrid,nlayer),pq(ngrid,nlayer,nq)
+      REAL pdu(ngrid,nlayer),pdv(ngrid,nlayer)
+
+!c   dynamial tendencies
+      REAL pdtdyn(ngrid,nlayer),pdqdyn(ngrid,nlayer,nq)
+      REAL pdudyn(ngrid,nlayer),pdvdyn(ngrid,nlayer)
+      REAL pw(ngrid,nlayer)
+
+!c   Time
+      real rjourvrai
+      REAL gmtime
+
+!c     outputs:
+!c     --------
+
+!c   physical tendencies
+      REAL pdt(ngrid,nlayer),pdq(ngrid,nlayer,nq)
+      REAL pdpsrf(ngrid)
+
+
+!c    Local variables :
+!c    -----------------
+
+      INTEGER j,l,ig,ierr,aslun,nlevel,igout,it1,it2,unit,isoil
+      REAL zps_av
+      REAL zday
+      REAL zh(ngrid,nlayer),z1,z2
+      REAL zzlev(ngrid,nlayer+1),zzlay(ngrid,nlayer)
+      REAL zdvfr(ngrid,nlayer),zdufr(ngrid,nlayer)
+      REAL zdhfr(ngrid,nlayer),zdtsrf(ngrid),zdtsrfr(ngrid)
+      REAL zflubid(ngrid),zpmer(ngrid)
+      REAL zplanck(ngrid),zpopsk(ngrid,nlayer)
+      REAL zdum1(ngrid,nlayer)
+      REAL zdum2(ngrid,nlayer)
+      REAL zdum3(ngrid,nlayer)
+      REAL ztim1,ztim2,ztim3
+      REAL zls,zinsol
+      REAL zdtlw(ngrid,nlayer),zdtsw(ngrid,nlayer)
+      REAL zfluxsw(ngrid),zfluxlw(ngrid)
+      character*2 str2
+
+!c   Local saved variables:
+!c   ----------------------
+	REAL(rstd)::long(ngrid),lati(ngrid) 
+	REAL(rstd)::mu0(ngrid),fract(ngrid),coslat(ngrid)
+        REAL(rstd)::sinlon(ngrid),coslon(ngrid),sinlat(ngrid)
+	REAL(rstd)::dist_sol,declin
+        REAL::totarea !sarvesh  
+!!!!!!!!sarvesh !!!!!!! CHECK SAVE ATTRIBUTE 
+      INTEGER:: icount
+      real:: zday_last
+      REAL:: solarcst
+
+      SAVE icount,zday_last
+      SAVE solarcst
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+      REAL stephan
+      SAVE stephan
+      DATA stephan/5.67e-08/
+      DATA solarcst/1370./
+      REAL presnivs(nlayer)
+	INTEGER:: nn1,nn2
+c-----------------------------------------------------------------------
+c    1. Initialisations :
+c    --------------------
+c     call initial0(ngrid*nlayer*nqmx,pdq)
+
+!	nn1=(jj_begin -1)*iim+ii_begin
+!	nn2=(jj_end -1)*iim+ii_end
+
+      nlevel=nlayer+1
+      igout=ngrid/2+1
+
+	 DO j=jj_begin-offset,jj_end+offset
+         DO i=ii_begin-offset,ii_end+offset
+           ig=(j-1)*iim+i
+	   sinlat(ig) = sin(lati(ig))
+	   coslat(ig) = cos(lati(ig))
+	   sinlon(ig) = sin(long(ig))
+	   coslon(ig) = cos(long(ig))
+	  END DO
+	 ENDDO 
+        zday=rjourvrai+gmtime
+	IF ( it == 0 ) then 
+	 firstcall=.TRUE.
+	ELSE
+	 firstcall=.FALSE. 
+	ENDIF 
+	IF ( it == ndays*day_step ) Then
+	 lastcall = .True.
+        END IF 
+	 
+        IF(firstcall) THEN
+        PRINT*,'FIRSTCALL  ',ngridmx,nlayermx,nsoilmx
+         zday_last=rjourvrai
+         inertie=2000
+         albedo=0.2
+         emissiv=1.
+         z0=0.1
+         rnatur=1.
+         q2=1.e-10
+         q2l=1.e-10
+         tsurf(:)=300.
+         tsoil(:,:)=300.
+!         print*,tsoil(ngrid/2+1,nsoilmx/2+2)
+!         print*,'TS ',tsurf(igout),tsoil(igout,5)
+ 
+          IF (.not.callrad) then
+	    DO j=jj_begin-offset,jj_end+offset
+             DO i=ii_begin-offset,ii_end+offset
+              ig=(j-1)*iim+i
+              fluxrad(ig)=0.
+             enddo
+            enddo
+           ENDIF
+
+         IF(callsoil) THEN
+             CALL soil(ngrid,nsoilmx,firstcall,inertie,
+     s          ptimestep,tsurf,tsoil,capcal,fluxgrd)
+          ELSE
+            PRINT*,'WARNING!!! Thermal conduction in the soil
+     s      turned off'
+	    DO j=jj_begin-offset,jj_end+offset
+             DO i=ii_begin-offset,ii_end+offset
+               ig=(j-1)*iim+i
+               capcal(ig)=1.e5
+               fluxgrd(ig)=0.
+             ENDDO
+           ENDDO
+          ENDIF
+            icount=0
+         ENDIF
+
+        IF (ngrid.NE.ngrid) THEN
+         PRINT*,'STOP in inifis'
+         PRINT*,'Probleme de dimenesions :'
+         PRINT*,'ngrid     = ',ngrid
+         PRINT*,'ngrid   = ',ngrid
+         STOP
+        ENDIF
+
+       DO l=1,nlayer
+         DO j=jj_begin-offset,jj_end+offset
+         DO i=ii_begin-offset,ii_end+offset
+            ig=(j-1)*iim+i
+            pdv(ig,l)=0.
+            pdu(ig,l)=0.
+            pdt(ig,l)=0.
+          ENDDO
+         ENDDO
+        ENDDO
+
+         DO j=jj_begin-offset,jj_end+offset
+         DO i=ii_begin-offset,ii_end+offset
+          ig=(j-1)*iim+i
+          pdpsrf(ig)=0.
+          zflubid(ig)=0.
+          zdtsrf(ig)=0.
+         ENDDO
+	ENDDO
+
+        zps_av=0.0
+         DO j=jj_begin-offset,jj_end+offset
+         DO i=ii_begin-offset,ii_end+offset
+           ig=(j-1)*iim+i
+           zps_av=zps_av+pplev(ig,1)*Ai(ig)
+           totarea=totarea+Ai(ig)
+          END DO
+         END DO
+	zps_av=zps_av/totarea
+
+
+	!print*,"maxplev",maxval(pplev(:,1)),minval(pplev(:,1))
+c-----------------------------------------------------------------------
+c   calcul du geopotentiel aux niveaux intercouches
+c   ponderation des altitudes au niveau des couches en dp/p
+
+       DO l=1,nlayer
+         DO j=jj_begin-offset,jj_end+offset
+         DO i=ii_begin-offset,ii_end+offset
+          ig=(j-1)*iim+i
+          zzlay(ig,l)=pphi(ig,l)/g
+         ENDDO
+        ENDDO
+       ENDDO
+
+	!print*,"zzlay",maxval(zzlay(:,1)),minval(zzlay(:,1))
+
+
+        DO j=jj_begin-offset,jj_end+offset
+         DO i=ii_begin-offset,ii_end+offset
+          ig=(j-1)*iim+i
+          zzlev(ig,1)=0.
+         ENDDO
+        ENDDO
+
+      DO l=2,nlayer
+         DO j=jj_begin-offset,jj_end+offset
+         DO i=ii_begin-offset,ii_end+offset
+          ig=(j-1)*iim+i
+          z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
+          z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
+          zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
+         ENDDO
+        ENDDO
+       ENDDO
+       !print*,"zzlev",maxval(zzlev(:,1)),minval(zzlev(:,1))
+c-----------------------------------------------------------------------
+c   Transformation de la temperature en temperature potentielle
+        DO l=1,nlayer
+          DO j=jj_begin-offset,jj_end+offset
+          DO i=ii_begin-offset,ii_end+offset
+           ig=(j-1)*iim  +i
+           zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**kappa
+          ENDDO
+         ENDDO
+        ENDDO
+
+        DO l=1,nlayer
+          DO j=jj_begin-offset,jj_end+offset
+          DO i=ii_begin-offset,ii_end+offset
+           ig=(j-1)*iim+i
+           zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
+          ENDDO
+         ENDDO
+	ENDDO
+	!print*,"ph pot",maxval(zh(:,1)),minval(zh(:,1))
+!	go to 101
+c-----------------------------------------------------------------------
+c    2. Calcul of the radiative tendencies :
+c    ---------------------------------------
+!      print*,'callrad0'
+      IF(callrad) THEN
+!         print*,'callrad'
+!   WARNING !!! on calcule le ray a chaque appel
+!        IF( MOD(icount,iradia).EQ.0) THEN
+
+           CALL solarlong(zday,zls)
+           CALL orbite(zls,dist_sol,declin)
+!   
+!      print*,'ATTENTIOn : pdeclin = 0','  L_s=',zls
+!      print*,'diurnal=',diurnal
+       IF(diurnal) THEN
+         if ( 1.eq.1 ) then
+               ztim1=SIN(declin)
+               ztim2=COS(declin)*COS(2.*pi*(zday-.5))
+               ztim3=-COS(declin)*SIN(2.*pi*(zday-.5))
+
+               CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
+     s         ztim1,ztim2,ztim3,
+     s         mu0,fract)
+          else
+               zdtime=ptimestep*float(iradia)
+               CALL zenang(zls,gmtime,zdtime,lati,long,mu0,fract)
+              !print*,'ZENANG '
+         endif
+
+          IF(lverbose) THEN
+                  PRINT*,'day, declin, sinlon,coslon,sinlat,coslat'
+                  PRINT*,zday, declin,
+     s            sinlon(igout),coslon(igout),
+     s            sinlat(igout),coslat(igout)
+           ENDIF
+        ELSE
+             !print*,'declin,ngrid,radius',declin,ngrid,radius
+            CALL mucorr(ngrid,declin,lati,mu0,fract,10000.,radius)
+!	      open(100,file="mu0.txt")
+!	      write(100,*)(mu0(ij),ij=1,ngrid) 
+        ENDIF
+!	print*,"iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii"
+!c    2.2 Calcul of the radiative tendencies and fluxes:
+!c    --------------------------------------------------
+!c  2.1.2 levels
+
+            zinsol=solarcst/(dist_sol*dist_sol)
+!            print*,'iim,jjm,llm,ngrid,nlayer,ngrid,nlayer'
+!            print*,iim,jjm,llm,ngrid,nlayer,ngrid,nlayer
+!		print*,"zinsol sol_dist",zinsol,dist_sol
+!	STOP 
+
+            CALL sw(ngrid,nlayer,diurnal,coefvis,albedo,
+     $              pplev,zps_av,
+     $              mu0,fract,zinsol,
+     $              zfluxsw,zdtsw,
+     $              lverbose)
+
+	!print*,"sw",maxval(zfluxsw),minval(zfluxsw),
+!	$            maxval(zdtsw),minval(zdtsw), "   it",it
+
+!	print*,"lllllllllllllllllllllllllllllllllllllllll"
+!	print*,"pplev",maxval(pplev),minval(pplev)
+!	print*,"zps,tsurf",zps_av,maxval(tsurf),minval(tsurf)
+!	print*,"pt",maxval(pt),minval(pt) 
+
+
+            CALL lw(ngrid,nlayer,coefir,emissiv,
+     $             pplev,zps_av,tsurf,pt,
+     $             zfluxlw,zdtlw,
+     $             lverbose)
+
+	!print*,"lw",maxval(zfluxlw),minval(zfluxlw),
+!	$            maxval(zdtlw),minval(zdtlw)
+
+!	print*,"lw",maxval(
+
+!	print*,"after lw"
+c    2.4 total flux and tendencies:
+c    ------------------------------
+
+c    2.4.1 fluxes
+
+          DO j=jj_begin-offset,jj_end+offset
+           DO i=ii_begin-offset,ii_end+offset
+               ig=(j-1)*iim+i
+               fluxrad(ig)=emissiv(ig)*zfluxlw(ig)
+     $         +zfluxsw(ig)*(1.-albedo(ig))
+
+               zplanck(ig)=tsurf(ig)*tsurf(ig)
+
+               zplanck(ig)=emissiv(ig)*
+     $         stephan*zplanck(ig)*zplanck(ig)
+
+               fluxrad(ig)=fluxrad(ig)-zplanck(ig)
+            ENDDO
+          ENDDO
+c    2.4.2 temperature tendencies
+
+            DO l=1,nlayer
+	       DO j=jj_begin-offset,jj_end+offset
+                DO i=ii_begin-offset,ii_end+offset
+                  ig=(j-1)*iim+i
+                  dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
+                 ENDDO
+               ENDDO
+             ENDDO
+
+
+c    2.5 Transformation of the radiative tendencies:
+c    -----------------------------------------------
+
+         DO l=1,nlayer
+	    DO j=jj_begin-offset,jj_end+offset
+            DO i=ii_begin-offset,ii_end+offset
+              ig=(j-1)*iim+i
+              pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
+            ENDDO
+           ENDDO
+          ENDDO
+
+         IF(lverbose) THEN
+            PRINT*,'Diagnotique for the radaition'
+            PRINT*,'albedo, emissiv, mu0,fract,Frad,Planck'
+            PRINT*,albedo(igout),emissiv(igout),mu0(igout),
+     s           fract(igout),
+     s           fluxrad(igout),zplanck(igout)
+            PRINT*,'Tlay Play Plev dT/dt SW dT/dt LW (K/day)'
+!            PRINT*,'unjours',unjours
+            DO l=1,nlayer
+               WRITE(*,'(3f15.5,2e15.2)') pt(igout,l),
+     s         pplay(igout,l),pplev(igout,l),
+     s         zdtsw(igout,l),zdtlw(igout,l)
+            ENDDO
+          ENDIF
+        ENDIF   !( CALL RADIATION ) 
+!	print*,"eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee"
+
+c-----------------------------------------------------------------------
+c    3. Vertical diffusion (turbulent mixing):
+c    -----------------------------------------
+c
+       IF(calldifv) THEN
+
+         DO ig=1,ngrid
+            zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
+         ENDDO
+
+         CALL zerophys(ngrid*nlayer,zdum1)
+         CALL zerophys(ngrid*nlayer,zdum2)
+c        CALL zerophys(ngrid*nlayer,zdum3)
+         do l=1,nlayer
+            do ig=1,ngrid
+               zdum3(ig,l)=pdt(ig,l)/zpopsk(ig,l)
+            enddo
+         enddo
+
+         CALL vdif(ngrid,nlayer,zday,
+     $        ptimestep,capcal,z0,
+     $        pplay,pplev,zzlay,zzlev,
+     $        pu,pv,zh,tsurf,emissiv,
+     $        zdum1,zdum2,zdum3,zflubid,
+     $        zdufr,zdvfr,zdhfr,zdtsrfr,q2,q2l,
+     $        lverbose)
+c        igout=ngrid/2+1
+c        PRINT*,'zdufr zdvfr zdhfr'
+c        DO l=1,nlayer
+c           PRINT*,zdufr(igout,1),zdvfr(igout,l),zdhfr(igout,l)
+c        ENDDO
+c        CALL difv  (ngrid,nlayer,
+c    $                  area,lati,capcal,
+c    $                  pplev,pphi,
+c    $                  pu,pv,zh,tsurf,emissiv,
+c    $                  zdum1,zdum2,zdum3,zflubid,
+c    $                  zdufr,zdvfr,zdhfr,zdtsrf,
+c    $                  .true.)
+c        PRINT*,'zdufr zdvfr zdhfr'
+c        DO l=1,nlayer
+c           PRINT*,zdufr(igout,1),zdvfr(igout,l),zdhfr(igout,l)
+c        ENDDO
+c        STOP
+
+         DO l=1,nlayer
+            DO ig=1,ngrid
+               pdv(ig,l)=pdv(ig,l)+zdvfr(ig,l)
+               pdu(ig,l)=pdu(ig,l)+zdufr(ig,l)
+               pdt(ig,l)=pdt(ig,l)+zdhfr(ig,l)*zpopsk(ig,l)
+            ENDDO
+         ENDDO
+
+         DO ig=1,ngrid
+            zdtsrf(ig)=zdtsrf(ig)+zdtsrfr(ig)
+         ENDDO
+
+       ELSE
+
+         DO j=jj_begin-offset,jj_end+offset
+          DO i=ii_begin-offset,ii_end+offset
+            ig=(j-1)*iim+i
+            zdtsrf(ig)=zdtsrf(ig)+
+     s      (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
+          ENDDO
+         ENDDO
+
+!       write(66,*)"tsrf",maxval(zdtsrf(:)),minval(zdtsrf(:))
+!       write(66,*)"frd",maxval(fluxrad(:)),minval(fluxrad(:))
+!       write(66,*)"fgd",maxval(fluxgrd(:)),minval(fluxgrd(:))
+        ENDIF
+c-----------------------------------------------------------------------
+c   4. Dry convective adjustment:
+c   -----------------------------
+
+      IF(calladj) THEN
+
+         DO l=1,nlayer
+            DO ig=1,ngrid
+               zdum1(ig,l)=pdt(ig,l)/zpopsk(ig,l)
+            ENDDO
+         ENDDO
+         CALL zerophys(ngrid*nlayer,zdufr)
+         CALL zerophys(ngrid*nlayer,zdvfr)
+         CALL zerophys(ngrid*nlayer,zdhfr)
+         CALL convadj(ngrid,nlayer,ptimestep,
+     $                pplay,pplev,zpopsk,
+     $                pu,pv,zh,
+     $                pdu,pdv,zdum1,
+     $                zdufr,zdvfr,zdhfr)
+
+         DO l=1,nlayer
+            DO ig=1,ngrid
+               pdu(ig,l)=pdu(ig,l)+zdufr(ig,l)
+               pdv(ig,l)=pdv(ig,l)+zdvfr(ig,l)
+               pdt(ig,l)=pdt(ig,l)+zdhfr(ig,l)*zpopsk(ig,l)
+            ENDDO
+         ENDDO
+
+        ENDIF
+
+!101		continue
+c-----------------------------------------------------------------------
+c   On incremente les tendances physiques de la temperature du sol:
+c   ---------------------------------------------------------------
+!       WRITE(55,*)"tsurf",maxval(tsurf(:)),minval(tsurf(:)),it
+
+	  DO j=jj_begin-offset,jj_end+offset
+           DO i=ii_begin-offset,ii_end+offset
+             ig=(j-1)*iim+i
+             tsurf(ig)=tsurf(ig)+ptimestep*zdtsrf(ig) 
+           ENDDO
+         ENDDO
+c-----------------------------------------------------------------------
+c   soil temperatures:
+c   --------------------
+
+       IF (callsoil) THEN
+         CALL soil(ngrid,nsoilmx,.false.,inertie,
+     s          ptimestep,tsurf,tsoil,capcal,fluxgrd)
+         IF(lverbose) THEN
+            PRINT*,'Surface Heat capacity,conduction Flux, Ts,
+     s          dTs, dt'
+            PRINT*,capcal(igout),fluxgrd(igout),tsurf(igout),
+     s          zdtsrf(igout),ptimestep
+         ENDIF
+       ENDIF
+
+c-----------------------------------------------------------------------
+c   sorties:
+c   --------
+      if(zday.GT.zday_last+period_sort) then
+       zday_last=zday
+c  Ecriture/extension de la coordonnee temps
+
+         do ig=1,ngrid
+            zpmer(ig)=pplev(ig,1)*exp(pphi(ig,1)/(kappa*cpp*285.))
+         enddo
+       endif
+c-----------------------------------------------------------------------
+      IF(lastcall) THEN
+       PRINT*,'Ecriture du fichier de reinitialiastion de la physique'
+      ENDIF
+
+      icount=icount+1
+!      RETURN
+      END SUBROUTINE phyparam_lmd
+	END MODULE PHY
Index: codes/icosagcm/trunk/src/physics.f90
===================================================================
--- codes/icosagcm/trunk/src/physics.f90	(revision 148)
+++ codes/icosagcm/trunk/src/physics.f90	(revision 149)
@@ -8,5 +8,6 @@
   SUBROUTINE init_physics
   USE icosa
-  USE physics_dcmip_mod, init_physics_dcmip=>init_physics
+  USE physics_dcmip_mod,init_physics_dcmip=>init_physics
+  USE physics_dry_mod
   IMPLICIT NONE
     
@@ -18,18 +19,25 @@
       CASE ('dcmip')
         CALL init_physics_dcmip
+
+      CASE ('lmd')
+	CALL init_physics_dry
       
       CASE DEFAULT
-         PRINT*, 'Bad selector for variable physics <',physics_type, &
+         PRINT*, 'Bad selector for variable physics init <',physics_type, &
               '> options are <none>, <dcmip>,'
-         STOP
+
     END SELECT
     
   END SUBROUTINE init_physics
 
-  SUBROUTINE physics(it,f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
+  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
     INTEGER, INTENT(IN)   :: it
+    REAL(rstd),INTENT(IN)::jD_cur,jH_cur
     TYPE(t_field),POINTER :: f_phis(:)
     TYPE(t_field),POINTER :: f_ps(:)
@@ -37,15 +45,29 @@
     TYPE(t_field),POINTER :: f_ue(:)
     TYPE(t_field),POINTER :: f_q(:)
+    LOGICAL:: firstcall,lastcall
     
     SELECT CASE(TRIM(physics_type))
       CASE ('none')
+
+	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)
+
+      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
+  STOP
     END SELECT
     
Index: codes/icosagcm/trunk/src/physics_dry.f90
===================================================================
--- codes/icosagcm/trunk/src/physics_dry.f90	(revision 149)
+++ codes/icosagcm/trunk/src/physics_dry.f90	(revision 149)
@@ -0,0 +1,307 @@
+	MODULE physics_dry_mod
+  USE ICOSA
+  PUBLIC init_physics_dry, physics_dry
+
+CONTAINS
+
+   SUBROUTINE init_physics_dry
+      USE ICOSA !sarvesh
+      USE time_mod !sarvesh 
+      USE dimphys_mod
+      USE RADIATION
+
+      IMPLICIT NONE
+      INTEGER::i,j,ij
+!-------------------------------------- ORBITAL PARAMETER----
+    periheli=150.
+    CALL getin('periheli', periheli)
+    aphelie=150.
+    CALL getin('aphelie',aphelie) 
+    coefir=0.08
+    CALL getin('coefir',coefir)
+    coefvis=0.99
+    CALL getin('coefvis',coefvis)
+    obliquit=0.0
+    CALL getin('obliquit',obliquit) 
+    peri_day=0.
+    CALL getin('peri_day',peri_day)
+    year_day=360.
+    CALL getin('year_day',year_day) 
+    callrad=.true.
+    CALL getin('callrad', callrad)
+    calldifv=.true.
+    CALL getin('calldifv', calldifv)
+    calladj=.true.
+    CALL getin('calladj', calladj)
+    callcond=.true.
+    callsoil=.true.
+    CALL getin('callsoil',callsoil)
+    season=.true.
+    CALL getin('season',season)
+    diurnal=.true.
+    CALL getin('diurnal',diurnal)
+    lverbose=.false.
+    CALL getin('lverbose',lverbose)
+    period_sort=1.
+   CALL getin('period_sort',period_sort)
+!    ptimestep=dt
+!   CALL getin('ptimestep',ptimestep) 
+
+      print*,'Activation de la physique:'
+      print*,' Rayonnement ',callrad
+      print*,' Diffusion verticale turbulente ', calldifv
+      print*,' Ajustement convectif ',calladj
+      print*,' Sol ',callsoil
+      print*,' Cycle diurne ',diurnal
+!      choice of the frequency of the computation of radiations
+      IF(diurnal) THEN
+         iradia=NINT(daysec/(20.*dt))
+      ELSE
+         iradia=NINT(daysec/(4.*dt))
+      ENDIF
+      iradia=1
+
+      ngridmx=iim*jjm ; nlayermx=llm
+	offset=halo
+
+      ALLOCATE(albedo(ngridmx));ALLOCATE(emissiv(ngridmx))
+      ALLOCATE(inertie(ngridmx));ALLOCATE(z0(ngridmx))
+      ALLOCATE(rnatur(ngridmx));ALLOCATE(tsurf(ngridmx))
+      ALLOCATE(tsoil(ngridmx,nlayermx));ALLOCATE(fluxgrd(ngridmx))
+      ALLOCATE(fluxrad(ngridmx));ALLOCATE(dtrad(ngridmx,llm+1))
+      ALLOCATE(q2(ngridmx,llm+1));ALLOCATE(q2l(ngridmx,llm+1))
+      ALLOCATE(capcal(ngridmx))
+
+      CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
+
+      PRINT*,'unjours',daysec
+      PRINT*,'The radiative transfer is computed each ', iradia,' physical time-step or each ', &
+              iradia*dt,' seconds'
+  END SUBROUTINE init_physics_dry
+
+
+  SUBROUTINE physics_dry( it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
+  USE icosa
+  IMPLICIT NONE
+    INTEGER,INTENT(IN)    :: it
+    REAL(rstd),INTENT(IN) :: jD_cur,jH_cur
+    TYPE(t_field),POINTER :: f_phis(:)
+    TYPE(t_field),POINTER :: f_ps(:)
+    TYPE(t_field),POINTER :: f_theta_rhodz(:)
+    TYPE(t_field),POINTER :: f_ue(:)
+    TYPE(t_field),POINTER :: f_q(:)
+  
+    REAL(rstd),POINTER :: phis(:)
+    REAL(rstd),POINTER :: ps(:)
+    REAL(rstd),POINTER :: theta_rhodz(:,:)
+    REAL(rstd),POINTER :: ue(:,:)
+    REAL(rstd),POINTER :: q(:,:,:)
+!    REAL(rstd),POINTER :: precl(:)
+    INTEGER :: ind
+!    LOGICAL:: firstcall,lastcall
+
+    CALL transfert_request(f_ue,req_e1_vect)
+    CALL transfert_request(f_theta_rhodz,req_i1)
+
+    DO ind=1,ndomain
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      phis=f_phis(ind)
+      ps=f_ps(ind)
+      theta_rhodz=f_theta_rhodz(ind)
+      ue=f_ue(ind)
+      q=f_q(ind)
+!      out_i=f_out_i(ind)
+!      precl=f_precl(ind) 
+!     print*,"====================================ind",ind,"----------it",it
+      CALL compute_physics_dry(it,jD_cur,jH_cur,phis, ps, theta_rhodz, ue, q(:,:,1))
+     ENDDO	
+
+!    CALL writefield("out_i",f_out_i)
+    
+!    IF (mod(it,itau_out)==0 ) THEN
+!      CALL writefield("precl",f_precl)
+!    ENDIF
+
+  END SUBROUTINE physics_dry
+  
+  SUBROUTINE compute_physics_dry(it,jD_cur,jH_cur,phis, ps, theta_rhodz, ue, q)
+  USE icosa
+  USE pression_mod
+  USE exner_mod
+  USE theta2theta_rhodz_mod
+  USE geopotential_mod
+  USE wind_mod
+  USE PHY
+  
+  IMPLICIT NONE
+    INTEGER::it
+    REAL(rstd) :: jD_cur
+    REAL(rstd) :: jH_cur
+    REAL(rstd) :: phis(iim*jjm)
+    REAL(rstd) :: ps(iim*jjm)
+    REAL(rstd) :: theta_rhodz(iim*jjm,llm)
+    REAL(rstd) :: ue(3*iim*jjm,llm)
+    REAL(rstd) :: q(iim*jjm,llm)
+!    REAL(rstd) :: precl(iim*jjm)
+
+    REAL(rstd) :: p(iim*jjm,llm+1)
+    REAL(rstd) :: pks(iim*jjm)
+    REAL(rstd) :: pk(iim*jjm,llm)
+    REAL(rstd) :: phi(iim*jjm,llm)
+    REAL(rstd) :: T(iim*jjm,llm)
+    REAL(rstd) :: Tfi(iim*jjm,llm)
+    REAL(rstd) :: theta(iim*jjm,llm)
+
+   REAL(rstd) :: uc(iim*jjm,3,llm)
+   REAL(rstd) :: u(iim*jjm,llm)
+   REAL(rstd) :: v(iim*jjm,llm)
+    REAL(rstd) :: ufi(iim*jjm,llm)
+    REAL(rstd) :: vfi(iim*jjm,llm)
+    REAL(rstd) :: qfi(iim*jjm,llm)
+    REAL(rstd) :: utemp(iim*jjm,llm)
+    REAL(rstd) :: vtemp(iim*jjm,llm)
+    REAL(rstd) :: lat(iim*jjm)
+    REAL(rstd) :: lon(iim*jjm)
+    REAL(rstd) :: pmid(iim*jjm,llm)
+    REAL(rstd) :: pint(iim*jjm,llm+1)
+    REAL(rstd) :: pdel(iim*jjm,llm)
+    REAL(rstd) :: plev(iim*jjm,llm+1),play(iim*jjm,llm) 
+    REAL(rstd) :: pkbycp
+    INTEGER :: i,j,l,ij,ig
+	
+!-------------------
+!    LOGICAL:: firstcall,lastcall
+    REAL(rstd) :: dufi(iim*jjm,llm)
+    REAL(rstd) :: dvfi(iim*jjm,llm)
+    REAL(rstd) :: dTfi(iim*jjm,llm)
+    REAL(rstd) :: dpsfi(iim*jjm)
+    REAL(rstd) :: dqfi(iim*jjm,llm)
+!    PRINT *,'Entering in LMD SIMPLE physics'    
+
+
+	offset=halo 
+    CALL compute_pression(ps,p,halo)
+    CALL compute_exner(ps,p,pks,pk,halo)
+    CALL compute_theta_rhodz2theta(ps,theta_rhodz,theta,halo)
+    CALL compute_geopotential(phis,pks,pk,theta,phi,halo)
+    CALL compute_theta_rhodz2temperature(ps,theta_rhodz,T,halo)
+    CALL compute_wind_centered(ue,uc)
+    CALL compute_wind_centered_lonlat_compound(uc, u, v)
+
+     DO j=jj_begin-offset,jj_end+offset
+      DO i=ii_begin-offset,ii_end+offset
+        ij=(j-1)*iim+i
+        CALL xyz2lonlat(xyz_i(ij,:),lon(ij),lat(ij)) 
+      ENDDO
+     ENDDO	
+   
+    DO l=1,llm
+      DO j=jj_begin-offset,jj_end+offset
+        DO i=ii_begin-offset,ii_end+offset
+          ij=(j-1)*iim+i
+!          Tfi(ij,l)=T(ij,l)
+!          ufi(ij,l)=u(ij,l)
+!          vfi(ij,l)=v(ij,l)
+!          qfi(ij,l)=q(ij,l)
+          dTfi(ij,l)=0.0
+          dufi(ij,l)=0.0
+          dvfi(ij,l)=0.0
+          dqfi(ij,l)=0.0
+        ENDDO
+      ENDDO
+    ENDDO    
+	plev(:,:) = p(:,:)  
+	dpsfi=0.0  
+
+   DO l=1,llm
+     DO j=jj_begin-offset,jj_end+offset
+       DO i=ii_begin-offset,ii_end+offset
+          ij=(j-1)*iim+i
+          pkbycp=pk(ij,l)/cpp
+          play(ij,l)=preff*pkbycp**(1./kappa) 
+        ENDDO
+      ENDDO
+     ENDDO 
+
+
+      CALL phyparam_lmd(it,iim*jjm,llm,1,dt,lat,lon,jD_cur,jH_cur, &
+      plev,play,phi,phis,u,v,T,q,dufi,dvfi,dTfi,dqfi,dpsfi) 
+
+      CALL ADDFI(u,v,T,q,ps,dufi,dvfi,dTfi,dqfi,dpsfi)
+
+!	CALL SARCHECKF(llm)
+!	print*,"plev",(maxval(plev(:,l)),l=1,llm+1) 
+
+!      CALL phyparam_lmd(it,iim*jjm,llm,1,dt,lat,lon,jD_cur,jH_cur, &
+!      plev,play,phi,phis,ufi,vfi,Tfi,qfi,dufi,dvfi,dTfi,dqfi,dpsfi) 
+
+!		Print*,"going ADD FI",it
+!      CALL ADDFI(ufi,vfi,Tfi,qfi,ps,dufi,dvfi,dTfi,dqfi,dpsfi)
+
+!       WRITE(11,*)"ducovfi",maxval(dufi),minval(dufi),it
+!       WRITE(11,*)"ucovfi",maxval(ufi),minval(ufi)
+!       WRITE(11,*)"dtetafi",maxval(dTfi),minval(dTfi)
+
+!=============================================
+!	go to 1234
+    DO l=1,llm
+      DO j=jj_begin-offset,jj_end+offset
+        DO i=ii_begin-offset,ii_end+offset
+          ij=(j-1)*iim+i
+          uc(ij,:,l)=(dufi(ij,l)*elon_i(ij,:)+dvfi(ij,l)*elat_i(ij,:))*dt
+        ENDDO
+      ENDDO
+    ENDDO
+
+    
+    DO l=1,llm
+      DO j=jj_begin-offset,jj_end+offset
+        DO i=ii_begin-offset,ii_end+offset
+          ij=(j-1)*iim+i
+          ue(ij+u_right,l)=ue(ij+u_right,l)+sum( 0.5*(uc(ij,:,l) + uc(ij+t_right,:,l))*ep_e(ij+u_right,:) )
+          ue(ij+u_lup,l)=ue(ij+u_lup,l)+sum( 0.5*(uc(ij,:,l) + uc(ij+t_lup,:,l))*ep_e(ij+u_lup,:) )
+          ue(ij+u_ldown,l)=ue(ij+u_ldown,l)+sum( 0.5*(uc(ij,:,l) + uc(ij+t_ldown,:,l))*ep_e(ij+u_ldown,:) )
+        ENDDO
+      ENDDO
+    ENDDO
+!1234	 continue
+   
+!   CALL compute_temperature2theta_rhodz(ps,Tfi,theta_rhodz,halo)
+   CALL compute_temperature2theta_rhodz(ps,T,theta_rhodz,halo)
+
+!       WRITE(13,*)"tetafi",maxval(theta_rhodz),minval(theta_rhodz)
+    RETURN   
+  END SUBROUTINE compute_physics_dry
+
+    SUBROUTINE addfi(ufi,vfi,Tfi,qfi,ps,dufi,dvfi,dTfi,dqfi,dpsfi)
+	USE ICOSA
+	IMPLICIT NONE
+	  REAL(rstd) :: dufi(iim*jjm,llm)
+	  REAL(rstd) :: dvfi(iim*jjm,llm)
+        REAL(rstd) :: dTfi(iim*jjm,llm)
+        REAL(rstd) :: dpsfi(iim*jjm)
+        REAL(rstd) :: dqfi(iim*jjm,llm)
+        REAL(rstd) :: ufi(iim*jjm,llm)
+        REAL(rstd) :: vfi(iim*jjm,llm)
+        REAL(rstd) :: qfi(iim*jjm,llm)
+	  REAL(rstd) :: ps(iim*jjm)
+	  REAL(rstd) :: Tfi(iim*jjm,llm)
+	  INTEGER::i,j,l,ij,offset
+		offset=halo
+	
+     DO l=1,llm
+      DO j=jj_begin-offset,jj_end+offset
+        DO i=ii_begin-offset,ii_end+offset
+          ij=(j-1)*iim+i
+          Tfi(ij,l)=Tfi(ij,l)+dTfi(ij,l)*dt
+          ufi(ij,l)=ufi(ij,l)+dufi(ij,l)*dt 
+          vfi(ij,l)=vfi(ij,l)+dvfi(ij,l)*dt 
+          qfi(ij,l)=qfi(ij,l)+dqfi(ij,l)*dt  
+        END DO
+      END DO
+     END DO  
+	 ps(:)=ps(:) + dpsfi(:)*dt
+	END SUBROUTINE addfi 
+   
+END MODULE physics_dry_mod
Index: codes/icosagcm/trunk/src/radiation_mod.F
===================================================================
--- codes/icosagcm/trunk/src/radiation_mod.F	(revision 149)
+++ codes/icosagcm/trunk/src/radiation_mod.F	(revision 149)
@@ -0,0 +1,1033 @@
+	MODULE RADIATION
+	 USE ICOSA
+	 USE dimphys_mod
+!	 USE PHY
+	contains
+
+!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	 SUBROUTINE zerophys(n,x)
+      IMPLICIT NONE
+      INTEGER n
+      REAL x(n)
+      INTEGER i
+
+      DO i=1,n
+         x(i)=0.
+      ENDDO
+      RETURN
+      END subroutine zerophys
+!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+	 SUBROUTINE solarlong(pday,psollong)
+      IMPLICIT NONE
+c=======================================================================
+c
+c   Objet:
+c   ------
+c
+c      Calcul de la distance soleil-planete et de la declinaison
+c   en fonction du jour de l'annee.
+c
+c
+c   Methode:
+c   --------
+c
+c      Calcul complet de l'elipse
+c
+c   Interface:
+c   ----------
+c
+c      Uncommon comprenant les parametres orbitaux.
+c
+c   Arguments:
+c   ----------
+c
+c   Input:
+c   ------
+c   pday          jour de l'annee (le jour 0 correspondant a l'equinoxe)
+c   lwrite        clef logique pour sorties de controle
+c
+c   Output:
+c   -------
+c   pdist_sol     distance entre le soleil et la planete
+c                 ( en unite astronomique pour utiliser la constante 
+c                  solaire terrestre 1370 Wm-2 )
+c   pdecli        declinaison ( en radians )
+c
+c=======================================================================
+c-----------------------------------------------------------------------
+c   Declarations:
+c   -------------
+
+c arguments:
+c ----------
+
+      REAL pday,pdist_sol,pdecli,psollong
+      LOGICAL lwrite
+
+c Local:
+c ------
+
+      REAL zanom,xref,zx0,zdx,zteta,zz
+      INTEGER iter
+
+
+c-----------------------------------------------------------------------
+c calcul de l'angle polaire et de la distance au soleil :
+c -------------------------------------------------------
+
+c  calcul de l'zanomalie moyenne
+
+      zz=(pday-peri_day)/year_day
+      zanom=2.*pi*(zz-nint(zz))
+      xref=abs(zanom)
+
+c  resolution de lequation horaire  zx0 - e * sin (zx0) = xref
+c  methode de Newton
+
+      zx0=xref+e_elips*sin(xref)
+      DO 110 iter=1,10
+         zdx=-(zx0-e_elips*sin(zx0)-xref)/(1.-e_elips*cos(zx0))
+         if(abs(zdx).le.(1.e-7)) goto 120
+         zx0=zx0+zdx
+110   continue
+120   continue
+      zx0=zx0+zdx
+      if(zanom.lt.0.) zx0=-zx0
+
+c zteta est la longitude solaire
+
+      zteta=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))
+
+      psollong=zteta-timeperi
+
+      IF(psollong.LT.0.) psollong=psollong+2.*pi
+      IF(psollong.GT.2.*pi) psollong=psollong-2.*pi
+
+      RETURN
+      END SUBROUTINE solarlong
+
+!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+        SUBROUTINE orbite(pls,pdist_sol,pdecli)
+      IMPLICIT NONE
+c====================================================================
+c
+c   Objet:
+c   ------
+c
+c   Distance from sun and declimation as a function of the solar
+c   longitude Ls
+c
+c   Interface:
+c   ----------
+c   Arguments:
+c   ----------
+c
+c   Input:
+c   ------
+c   pls          Ls
+c
+c   Output:
+c   -------
+c   pdist_sol     Distance Sun-Planet in UA
+c   pdecli        declinaison ( en radians )
+c
+c====================================================================
+c-----------------------------------------------------------------------
+c   Declarations:
+c   -------------
+c arguments:
+c ----------
+
+      REAL pday,pdist_sol,pdecli,pls
+
+c-----------------------------------------------------------------------
+
+c Distance Sun-Planet
+
+      pdist_sol=p_elips/(1.+e_elips*cos(pls+timeperi))
+
+c Solar declination
+
+      pdecli= asin (sin(pls)*sin(obliquit*pi/180.))
+
+c-----------------------------------------------------------------------
+c   sorties eventuelles:
+c   ---------------------
+
+      END SUBROUTINE orbite
+!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+        SUBROUTINE iniorbit
+     $     (paphelie,pperiheli,pyear_day,pperi_day,pobliq)
+      IMPLICIT NONE
+c=====================================================================
+c
+c   Auteur:
+c   -------
+c     Frederic Hourdin      22 Fevrier 1991
+c
+c   Objet:
+c   ------
+c    Initialisation du sous programme orbite qui calcule
+c    a une date donnee de l'annee de duree year_day commencant
+c    a l'equinoxe de printemps et dont le perihelie se situe
+c    a la date peri_day, la distance au soleil et la declinaison.
+c
+c   Interface:
+c   ----------
+c   - Doit etre appele avant d'utiliser orbite.
+c   - initialise le common comorbit
+c
+c   Arguments:
+c   ----------
+c
+c   Input:
+c   ------
+c   aphelie       \   aphelie et perihelie de l'orbite
+c   periheli      /   en millions de kilometres.
+c
+c=====================================================================
+c   Declarations:
+c   -------------
+c   Arguments:
+c   ----------
+
+      REAL paphelie,pperiheli,pyear_day,pperi_day,pobliq
+
+c   Local:
+c   ------
+
+      REAL zxref,zanom,zz,zx0,zdx
+      INTEGER iter
+
+c'-----------------------------------------------------------------------
+
+      aphelie =paphelie
+      periheli=pperiheli
+      year_day=pyear_day
+      obliquit=pobliq
+      peri_day=pperi_day
+
+      PRINT*,'Perihelie en Mkm  ',periheli
+      PRINT*,'Aphelise  en Mkm  ',aphelie 
+      PRINT*,'obliquite en degres  :',obliquit
+      unitastr=149.597927
+      e_elips=(aphelie-periheli)/(periheli+aphelie)
+      p_elips=0.5*(periheli+aphelie)*(1-e_elips*e_elips)/unitastr
+
+      print*,'e_elips',e_elips
+      print*,'p_elips',p_elips
+
+c-----------------------------------------------------------------------
+c calcul de l'angle polaire et de la distance au soleil :
+c -------------------------------------------------------
+
+c  calcul de l'zanomalie moyenne
+
+      zz=(year_day-pperi_day)/year_day
+      zanom=2.*pi*(zz-nint(zz))
+      zxref=abs(zanom)
+      PRINT*,'zanom  ',zanom
+
+c  resolution de lequation horaire  zx0 - e * sin (zx0) = zxref
+c  methode de Newton
+
+      zx0=zxref+e_elips*sin(zxref)
+      DO 110 iter=1,100
+         zdx=-(zx0-e_elips*sin(zx0)-zxref)/(1.-e_elips*cos(zx0))
+         if(abs(zdx).le.(1.e-12)) goto 120
+         zx0=zx0+zdx
+110   continue
+120   continue
+      zx0=zx0+zdx
+      if(zanom.lt.0.) zx0=-zx0
+      PRINT*,'zx0   ',zx0
+
+c zteta est la longitude solaire
+
+      timeperi=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))
+      PRINT*,'longitude solaire du perihelie timeperi = ',timeperi
+
+      RETURN
+      END SUBROUTINE iniorbit
+!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+      SUBROUTINE mucorr(npts,pdeclin, plat, pmu, pfract,phaut,prad)
+      IMPLICIT NONE
+
+c=======================================================================
+c
+c   Calcul of equivalent solar angle and and fraction of day whithout 
+c   diurnal cycle.
+c
+c   parmeters :
+c   -----------
+c
+c      Input :
+c      -------
+c         npts             number of points
+c         pdeclin          solar declinaison
+c         plat(npts)        latitude
+c         phaut            hauteur typique de l'atmosphere
+c         prad             rayon planetaire '
+c
+c      Output :
+c      --------
+c         pmu(npts)          equivalent cosinus of the solar angle
+c         pfract(npts)       fractionnal day
+c
+c=======================================================================
+
+c-----------------------------------------------------------------------
+c
+c    0. Declarations :
+c    -----------------
+
+c     Arguments :
+c     -----------
+      INTEGER npts
+      REAL plat(npts),pmu(npts), pfract(npts)
+      REAL phaut,prad,pdeclin
+c
+c     Local variables :
+c     -----------------
+      INTEGER j,i,ij,ig
+      REAL pi
+      REAL z,cz,sz,tz,phi,cphi,sphi,tphi
+      REAL ap,a,t,b
+      REAL alph
+
+c-----------------------------------------------------------------------
+
+      !print*,'npts,pdeclin'
+      !print*,npts,pdeclin
+      pi = 4. * atan(1.0)
+      !print*,'PI=',pi
+      pi=2.*asin(1.)
+      z = pdeclin
+      cz = cos (z)
+      sz = sin (z)
+      !print*,'cz,sz',cz,sz
+
+!	DO j=jj_begin-offset,jj_end+offset
+!         DO i=ii_begin-offset,ii_end+offset
+!             ig=(j-1)*iim+i
+
+	  DO ig=1,npts
+             phi = plat(ig)
+             cphi = cos(phi)
+           if (cphi.le.1.e-9) cphi=1.e-9
+            sphi = sin(phi)
+            tphi = sphi / cphi
+               b = cphi * cz
+               t = -tphi * sz / cz
+               a = 1.0 - t*t
+              ap = a
+             IF(t.eq.0.) then
+               t=0.5*pi
+             ELSE
+                IF (a.lt.0.) a = 0.
+                  t = sqrt(a) / t
+                IF (t.lt.0.) then
+                  t = -atan (-t) + pi
+                 ELSE
+                  t = atan(t)
+                ENDIF
+             ENDIF
+   
+            pmu(ig) = (sphi*sz*t) / pi + b*sin(t)/pi
+            pfract(ig) = t / pi
+            IF (ap .lt.0.) then
+             pmu(ig) = sphi * sz
+             pfract(ig) = 1.0
+            ENDIF
+            IF (pmu(ig).le.0.0) pmu(ig) = 0.0
+             pmu(ig) = pmu(ig) / pfract(ig)
+            IF (pmu(ig).eq.0.) pfract(ig) = 0.
+	  ENDDO
+!	ENDDO
+c-----------------------------------------------------------------------
+c   correction de rotondite:
+c   ------------------------
+
+	      alph=phaut/prad
+!	DO j=jj_begin-offset,jj_end+offset
+!         DO i=ii_begin-offset,ii_end+offset
+!		           ig=(j-1)*iim+i
+	   DO ig = 1,npts
+            pmu(ig)=sqrt(1224.*pmu(ig)*pmu(ig)+1.)/35.
+
+c    $    (sqrt(alph*alph*pmu(ig)*pmu(ig)+2.*alph+1.)-alph*pmu(ig))
+	 ENDDO
+!	ENDDO
+
+      RETURN
+      END SUBROUTINE mucorr
+!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	
+	      SUBROUTINE sw(ngrid,nlayer,ldiurn,
+     $              coefvis,albedo,
+     $              plevel,ps_rad,pmu,pfract,psolarf0,
+     $              fsrfvis,dtsw,
+     $              lwrite)
+      IMPLICIT NONE
+c=======================================================================
+c
+c   Rayonnement solaire en atmosphere non diffusante avec un 
+c   coefficient d'absoprption gris.
+c'
+c=======================================================================
+c
+c   declarations:
+c   -------------
+c   arguments:
+c   ----------
+c
+      INTEGER ngrid,nlayer,i,j,ij
+      REAL albedo(ngrid),coefvis
+      REAL pmu(ngrid),pfract(ngrid)
+      REAL plevel(ngrid,nlayer+1),ps_rad
+      REAL psolarf0
+      REAL fsrfvis(ngrid),dtsw(ngrid,nlayer)
+      LOGICAL lwrite,ldiurn
+c
+c   variables locales:
+c   ------------------
+c
+
+      REAL zalb(ngrid),zmu(ngrid),zfract(ngrid)
+      REAL zplev(ngrid,nlayer+1)
+      REAL zflux(ngrid),zdtsw(ngrid,nlayer)
+
+      INTEGER ig,l,nlevel,ncount,igout
+	INTEGER,DIMENSION(ngrid)::index
+      REAL ztrdir(ngrid,nlayer+1),ztrref(ngrid,nlayer+1)
+      REAL zfsrfref(ngrid)
+      REAL z1(ngrid)
+      REAL zu(ngrid,nlayer+1)
+      REAL tau0
+
+      EXTERNAL SSUM
+      EXTERNAL ismax
+      REAL ismax
+
+      LOGICAL firstcall
+      SAVE firstcall
+      DATA firstcall/.true./
+
+c-----------------------------------------------------------------------
+c   1. initialisations:
+c   -------------------
+
+
+ 
+      IF (firstcall) THEN
+ 
+      IF (ngrid.NE.ngrid) THEN
+         PRINT*,'STOP in inifis'
+         PRINT*,'Probleme de dimenesions :'
+         PRINT*,'ngrid     = ',ngrid
+         PRINT*,'ngrid   = ',ngrid
+         STOP
+      ENDIF
+ 
+ 
+      IF (nlayer.NE.llm) THEN
+         PRINT*,'STOP in inifis'
+         PRINT*,'Probleme de dimenesions :'
+         PRINT*,'nlayer     = ',nlayer
+         PRINT*,'llm   = ',llm
+         STOP
+      ENDIF
+ 
+      ENDIF
+
+      nlevel=nlayer+1
+c-----------------------------------------------------------------------
+c   Definitions des tableaux locaux pour les points ensoleilles:
+c   ------------------------------------------------------------
+
+      IF (ldiurn) THEN
+         ncount=0
+        DO j=jj_begin-offset,jj_end+offset
+           DO i=ii_begin-offset,ii_end+offset
+            ig=(j-1)*iim+i
+            index(ig)=0
+           ENDDO
+	 ENDDO
+         DO j=jj_begin-offset,jj_end+offset
+          DO i=ii_begin-offset,ii_end+offset
+             ig=(j-1)*iim+i
+           IF(pfract(ig).GT.1.e-6) THEN
+             ncount=ncount+1
+             index(ncount)=ig
+           ENDIF
+          ENDDO
+         ENDDO
+!	SARVESH CHANGED NCOUNT TO NGRID TO BE CONSISTENT ???
+
+         CALL monGATHER(ncount,zfract,pfract,index)
+         CALL monGATHER(ncount,zmu,pmu,index)
+         CALL monGATHER(ncount,zalb,albedo,index)
+         DO l=1,nlevel
+            CALL monGATHER(ncount,zplev(1,l),plevel(1,l),index)
+         ENDDO
+      ELSE
+         ncount=ngrid
+         zfract(:)=pfract(:)
+         zmu(:)=pmu(:)
+         zalb(:)=albedo(:)
+         zplev(:,:)=plevel(:,:)
+      ENDIF
+
+c-----------------------------------------------------------------------
+c   calcul des profondeurs optiques integres depuis p=0:
+c   ----------------------------------------------------
+
+      tau0=-.5*log(coefvis)
+
+c calcul de la partie homogene de l'opacite
+c'
+      tau0=tau0/ps_rad
+
+
+      DO l=1,nlayer+1
+        DO j=jj_begin-offset,jj_end+offset
+          DO i=ii_begin-offset,ii_end+offset
+            ig=(j-1)*iim+i
+            zu(ig,l)=tau0*zplev(ig,l)
+           ENDDO
+         ENDDO
+        ENDDO
+
+c-----------------------------------------------------------------------
+c   2. calcul de la transmission depuis le sommet de l'atmosphere:
+c'   -----------------------------------------------------------
+
+       DO l=1,nlevel
+         DO j=jj_begin-offset,jj_end+offset
+           DO i=ii_begin-offset,ii_end+offset
+            ig=(j-1)*iim+i
+            ztrdir(ig,l)=exp(-zu(ig,l)/zmu(ig))
+           ENDDO
+          ENDDO
+        ENDDO
+
+      IF (lwrite) THEN
+         igout=ncount/2+1
+         PRINT*
+         PRINT*,'Diagnostique des transmission dans le spectre solaire'
+         PRINT*,'zfract, zmu, zalb'
+         PRINT*,zfract(igout), zmu(igout), zalb(igout)
+         PRINT*,'Pression, quantite d abs, transmission'
+         DO l=1,nlayer+1
+            PRINT*,zplev(igout,l),zu(igout,l),ztrdir(igout,l)
+         ENDDO
+      ENDIF
+
+c-----------------------------------------------------------------------
+c   3. taux de chauffage, ray. solaire direct:
+c   ------------------------------------------
+
+      DO l=1,nlayer
+        DO j=jj_begin-offset,jj_end+offset
+          DO i=ii_begin-offset,ii_end+offset
+            ig=(j-1)*iim+i
+            zdtsw(ig,l)=g*psolarf0*zfract(ig)*zmu(ig)*
+     $                     (ztrdir(ig,l+1)-ztrdir(ig,l))/
+     $                     (cpp*(zplev(ig,l)-zplev(ig,l+1)))
+           ENDDO
+          ENDDO
+       ENDDO
+      IF (lwrite) THEN
+         PRINT*
+         PRINT*,'Diagnostique des taux de chauffage solaires:'
+         PRINT*,' 1 taux de chauffage lie au ray. solaire  direct'
+         DO l=1,nlayer
+            PRINT*,zdtsw(igout,l)
+         ENDDO
+      ENDIF
+
+
+c-----------------------------------------------------------------------
+c   4. calcul du flux solaire arrivant sur le sol:
+c   ----------------------------------------------
+
+       DO j=jj_begin-offset,jj_end+offset
+        DO i=ii_begin-offset,ii_end+offset
+          ig=(j-1)*iim+i
+          z1(ig)=zfract(ig)*zmu(ig)*psolarf0*ztrdir(ig,1)
+          zflux(ig)=(1.-zalb(ig))*z1(ig)
+          zfsrfref(ig)=    zalb(ig)*z1(ig)
+        ENDDO
+       ENDDO
+      IF (lwrite) THEN
+         PRINT*
+         PRINT*,'Diagnostique des taux de chauffage solaires:'
+         PRINT*,' 2 flux solaire net incident sur le sol'
+         PRINT*,zflux(igout)
+      ENDIF
+
+
+c-----------------------------------------------------------------------
+c   5.calcul des traansmissions depuis le sol, cas diffus:
+c   ------------------------------------------------------
+
+      DO l=1,nlevel
+       DO j=jj_begin-offset,jj_end+offset
+        DO i=ii_begin-offset,ii_end+offset
+          ig=(j-1)*iim+i
+            ztrref(ig,l)=exp(-(zu(ig,1)-zu(ig,l))*1.66)
+         ENDDO
+        ENDDO
+       ENDDO
+
+      IF (lwrite) THEN
+         PRINT*
+         PRINT*,'Diagnostique des taux de chauffage solaires'
+         PRINT*,' 3 transmission avec les sol'
+         PRINT*,'niveau     transmission'
+         DO l=1,nlevel
+            PRINT*,l,ztrref(igout,l)
+         ENDDO
+      ENDIF
+
+c-----------------------------------------------------------------------
+c   6.ajout a l'echauffement de la contribution du ray. sol. reflechit: 
+c' -------------------------------------------------------------------
+
+      DO l=1,nlayer
+       DO j=jj_begin-offset,jj_end+offset
+        DO i=ii_begin-offset,ii_end+offset
+          ig=(j-1)*iim+i
+            zdtsw(ig,l)=zdtsw(ig,l)+
+     $      g*zfsrfref(ig)*(ztrref(ig,l+1)-ztrref(ig,l))/
+     $      (cpp*(zplev(ig,l+1)-zplev(ig,l)))
+         ENDDO
+        ENDDO
+       ENDDO
+
+c-----------------------------------------------------------------------
+c   10. sorites eventuelles:
+c   ------------------------
+
+      IF (lwrite) THEN
+         PRINT*
+         PRINT*,'Diagnostique des taux de chauffage solaires:'
+         PRINT*,' 3 taux de chauffage total'
+         DO l=1,nlayer
+            PRINT*,zdtsw(igout,l)
+         ENDDO
+      ENDIF
+
+      IF (ldiurn) THEN
+         CALL zerophys(ngrid,fsrfvis)
+         CALL monscatter(ncount,fsrfvis,index,zflux)
+         CALL zerophys(ngrid*nlayer,dtsw)
+         DO l=1,nlayer
+            CALL monscatter(ncount,dtsw(1,l),index,zdtsw(1,l))
+         ENDDO
+      ELSE
+         !print*,'NOT DIURNE'
+         fsrfvis(:)=zflux(:)
+         dtsw(:,:)=zdtsw(:,:)
+      ENDIF
+
+      RETURN
+      END  SUBROUTINE sw
+!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+      SUBROUTINE lw(ngrid,nlayer,coefir,emissiv,
+     $             pp,ps_rad,ptsurf,pt,
+     $             pfluxir,pdtlw,
+     $             lwrite)
+
+      IMPLICIT NONE
+c=======================================================================
+c
+c   calcul de l'evolution de la temperature sous l'effet du rayonnement
+c   infra-rouge.
+c   Pour simplifier, les transmissions sont precalculees et ne
+c   dependent que de l'altitude.
+c
+c   arguments:
+c   ----------
+c'
+c   entree:
+c   -------
+c      ngrid             nombres de points de la grille horizontale
+c      nlayer              nombre de couches
+c      ptsurf(ngrid)     temperature de la surface
+c      pt(ngrid,nlayer)    temperature des couches
+c      pp(ngrid,nlayer+1)  pression entre les couches
+c      lwrite            variable logique pour sorties
+c
+c   sortie:
+c   -------
+c      pdtlw(ngrid,nlayer) taux de refroidissement
+c      pfluxir(ngrid)    flux infrarouge sur le sol
+c
+c=======================================================================
+
+!c   declarations:
+!c   -------------
+!c   arguments:
+!c'   ----------
+
+      INTEGER ngrid,nlayer
+      REAL coefir,emissiv(ngrid),ps_rad
+      REAL ptsurf(ngrid),pt(ngrid,nlayer),pp(ngrid,nlayer+1)
+      REAL pdtlw(ngrid,nlayer),pfluxir(ngrid)
+      LOGICAL lwrite
+
+c   variables locales:
+c   ------------------
+
+      INTEGER nlevel,ilev,ig,i,il
+      REAL zplanck(ngridmx,nlayermx+1),zcoef
+      REAL zfluxup(ngridmx,nlayermx+1),zfluxdn(ngridmx,nlayermx+1)
+      REAL zflux(ngridmx,nlayermx+1)
+      REAL zlwtr1(ngridmx),zlwtr2(ngridmx)
+      REAL zup(ngridmx,nlayermx+1),zdup(ngridmx)
+      REAL stephan
+
+      LOGICAL lstrong
+      SAVE lstrong,stephan
+      DATA lstrong/.true./
+c-----------------------------------------------------------------------
+c   initialisations:
+c   ----------------
+
+      nlevel=nlayer+1
+      stephan=5.67e-08
+
+
+	pfluxir=0.0
+	pdtlw=0.0
+	!print*,"ngr,nlay,coefi",ngrid,nlayer,coefir
+c-----------------------------------------------------------------------
+c   2. calcul des quantites d'absorbants:
+c'   -------------------------------------
+
+c   absorption forte
+      IF(lstrong) THEN
+         DO ilev=1,nlevel
+            DO ig=1,ngrid
+               zup(ig,ilev)=pp(ig,ilev)*pp(ig,ilev)/(2.*g)
+            ENDDO
+         ENDDO
+         IF(lwrite) THEN
+            DO ilev=1,nlayer
+             PRINT*,' up(',ilev,')  =  ',zup(ngrid/2+1,ilev)
+            ENDDO
+         ENDIF
+         zcoef=-log(coefir)/sqrt(ps_rad*ps_rad/(2.*g))
+
+c   absorption faible
+      ELSE
+         DO ilev=1,nlevel
+            DO ig=1,ngrid
+               zup(ig,ilev)=pp(ig,ilev)
+            ENDDO
+         ENDDO
+         zcoef=-log(coefir)/ps_rad
+      ENDIF
+
+
+c-----------------------------------------------------------------------
+c   2. calcul de la fonction de corps noir:
+c   ---------------------------------------
+
+      DO ilev=1,nlayer
+         DO ig=1,ngrid
+            zplanck(ig,ilev)=pt(ig,ilev)*pt(ig,ilev)
+            zplanck(ig,ilev)=stephan*
+     $      zplanck(ig,ilev)*zplanck(ig,ilev)
+	 ENDDO
+      ENDDO
+
+c-----------------------------------------------------------------------
+c   4. flux descendants:
+c   --------------------
+
+      DO ilev=1,nlayer
+         DO ig=1,ngrid
+            zfluxdn(ig,ilev)=0.
+		 ENDDO
+         DO ig=1,ngrid
+            zdup(ig)=zup(ig,ilev)-zup(ig,nlevel)
+         ENDDO
+         CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)
+
+         DO il=nlayer,ilev,-1
+            zlwtr2(:)=zlwtr1(:)
+            DO ig=1,ngrid
+               zdup(ig)=zup(ig,ilev)-zup(ig,il)
+            ENDDO
+            CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)
+            DO ig=1,ngrid
+               zfluxdn(ig,ilev)=zfluxdn(ig,ilev)+
+     $         zplanck(ig,il)*(zlwtr1(ig)-zlwtr2(ig))
+			ENDDO
+		 ENDDO
+      ENDDO
+
+      DO ig=1,ngrid
+         zfluxdn(ig,nlevel)=0.
+         pfluxir(ig)=emissiv(ig)*zfluxdn(ig,1)
+      ENDDO
+
+      DO ig=1,ngrid
+         zfluxup(ig,1)=ptsurf(ig)*ptsurf(ig)
+         zfluxup(ig,1)=emissiv(ig)*stephan*zfluxup(ig,1)*zfluxup(ig,1)
+     $   +(1.-emissiv(ig))*zfluxdn(ig,1)
+      ENDDO
+
+c-----------------------------------------------------------------------
+c   3. flux montants:
+c   ------------------
+
+      DO ilev=1,nlayer
+         DO ig=1,ngrid
+            zdup(ig)=zup(ig,1)-zup(ig,ilev+1)
+         ENDDO
+         CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)
+         DO ig=1,ngrid
+            zfluxup(ig,ilev+1)=zfluxup(ig,1)*zlwtr1(ig)
+		 ENDDO
+         DO il=1,ilev
+            zlwtr2(:)=zlwtr1(:)
+            DO ig=1,ngrid
+               zdup(ig)=zup(ig,il+1)-zup(ig,ilev+1)
+            ENDDO
+            CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)
+            DO ig=1,ngrid
+               zfluxup(ig,ilev+1)=zfluxup(ig,ilev+1)+
+     $         zplanck(ig,il)*(zlwtr1(ig)-zlwtr2(ig))
+			ENDDO
+         ENDDO
+
+      ENDDO
+
+c-----------------------------------------------------------------------
+c   5. calcul des flux nets:
+c   ------------------------
+
+      DO ilev=1,nlevel
+         DO ig=1,ngrid
+            zflux(ig,ilev)=zfluxup(ig,ilev)-zfluxdn(ig,ilev)
+		 ENDDO
+      ENDDO
+
+c-----------------------------------------------------------------------
+c   6. Calcul des taux de refroidissement:
+c   --------------------------------------
+   
+      DO ilev=1,nlayer
+         DO ig=1,ngrid
+            pdtlw(ig,ilev)=(zflux(ig,ilev+1)-zflux(ig,ilev))*
+     $           g/(cpp*(pp(ig,ilev+1)-pp(ig,ilev)))
+		 ENDDO
+      ENDDO
+
+c-----------------------------------------------------------------------
+c   10. sorties eventuelles:
+c   ------------------------
+
+      IF (lwrite) THEN
+
+         PRINT*
+         PRINT*,'Diagnostique rayonnement thermique'
+         PRINT*
+         PRINT*,'temperature     ',
+     $   'flux montant    flux desc.     taux de refroid.'
+         i=ngrid/2+1
+         WRITE(6,9000) ptsurf(i)
+         DO ilev=1,nlayer
+            WRITE(6,'(i4,4e18.4)') ilev,pt(i,ilev),
+     $      zfluxup(i,ilev),zfluxdn(i,ilev),pdtlw(i,ilev)
+		 ENDDO
+         WRITE(6,9000) zfluxup(i,nlevel),zfluxdn(i,nlevel)
+
+      ENDIF
+
+c-----------------------------------------------------------------------
+
+      RETURN
+9000  FORMAT(4e18.4)
+      END SUBROUTINE lw 
+!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+	      subroutine solang ( kgrid,psilon,pcolon,psilat,pcolat,
+     &                    ptim1,ptim2,ptim3,pmu0,pfract )
+      IMPLICIT NONE
+
+C
+C**** *LW*   - ORGANIZES THE LONGWAVE CALCULATIONS
+C
+C     PURPOSE.
+C     --------
+C          CALCULATES THE SOLAR ANGLE FOR ALL THE POINTS OF THE GRID
+C     ==== INPUTS  ===
+C
+C PSILON(KGRID)   : SINUS OF THE LONGITUDE
+C PCOLON(KGRID)   : COSINUS OF THE LONGITUDE
+C PSILAT(KGRID)   : SINUS OF THE LATITUDE
+C PCOLAT(KGRID)   : COSINUS OF THE LATITUDE
+C PTIM1           : SIN(DECLI)
+C PTIM2           : COS(DECLI)*COS(TIME)
+C PTIM3           : SIN(DECLI)*SIN(TIME)
+C
+C     ==== OUTPUTS ===
+C
+C PMU0 (KGRID)    : SOLAR ANGLE
+C PFRACT(KGRID)   : DAY FRACTION OF THE TIME INTERVAL
+C
+C        IMPLICIT ARGUMENTS :   NONE
+C        --------------------
+C
+C     METHOD.
+C     -------
+C
+C     EXTERNALS.
+C     ----------
+C
+C         NONE
+C
+C     REFERENCE.
+C     ----------
+C
+C         RADIATIVE PROCESSES IN METEOROLOGIE AND CLIMATOLOGIE
+C         PALTRIDGE AND PLATT
+C
+C     AUTHOR.
+C     -------
+C        FREDERIC HOURDIN
+C
+C     MODIFICATIONS.
+C     --------------
+C        ORIGINAL :90-01-14
+C                  92-02-14 CALCULATIONS DONE THE ENTIER GRID (J.Polcher)
+C-----------------------------------------------------------------------
+C
+C     ------------------------------------------------------------------
+
+C-----------------------------------------------------------------------
+C
+C*      0.1   ARGUMENTS
+C             ---------
+C
+      INTEGER kgrid
+      REAL ptim1,ptim2,ptim3
+      REAL psilon(kgrid),pcolon(kgrid),pmu0(kgrid),pfract(kgrid)
+      REAL psilat(kgrid), pcolat(kgrid)
+C
+      INTEGER jl,i,j
+      REAL ztim1,ztim2,ztim3
+C------------------------------------------------------------------------
+C------------------------------------------------------------------------
+C---------------------------------------------------------------------
+C
+C--------------------------------------------------------------------
+C
+C*     1.     INITIALISATION
+C             --------------
+C
+!!!!!! 100  CONTINUE
+C
+      DO j=jj_begin-offset,jj_end+offset
+       DO i=ii_begin-offset,ii_end+offset
+         jl=(j-1)*iim+i
+        pmu0(jl)=0.
+        pfract(jl)=0.
+       ENDDO
+      ENDDO
+!C
+!C*     1.1     COMPUTATION OF THE SOLAR ANGLE
+!C              ------------------------------
+!C
+       DO j=jj_begin-offset,jj_end+offset
+         DO i=ii_begin-offset,ii_end+offset
+          jl=(j-1)*iim+i
+         ztim1=psilat(jl)*ptim1
+         ztim2=pcolat(jl)*ptim2
+         ztim3=pcolat(jl)*ptim3
+         pmu0(jl)=ztim1+ztim2*pcolon(jl)+ztim3*psilon(jl)
+        ENDDO
+      ENDDO
+!C
+!C*     1.2      DISTINCTION BETWEEN DAY AND NIGHT
+!C               ---------------------------------
+!C
+       DO j=jj_begin-offset,jj_end+offset
+        DO i=ii_begin-offset,ii_end+offset
+          jl=(j-1)*iim+i      
+          IF (pmu0(jl).gt.0.) THEN
+            pfract(jl)=1.
+	   ELSE
+            pmu0(jl)=0.
+	    pfract(jl)=0.
+          ENDIF
+         ENDDO
+        ENDDO
+      RETURN
+      END SUBROUTINE solang
+!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+	      SUBROUTINE monGATHER(n,a,b,index)
+c
+      IMPLICIT NONE
+C
+      INTEGER n,ng,index(n),i,j,ij
+      REAL a(n),b(n)
+c
+      DO 100 i=1,n
+        a(i)=b(index(i))
+100   CONTINUE
+
+!       DO j=jj_begin-offset,jj_end+offset
+!         DO i=ii_begin-offset,ii_end+offset
+!          ij=(j-1)*iim+i
+!          a(ij)=b(index(ij))
+!c
+      RETURN
+      END  SUBROUTINE monGATHER
+!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+       subroutine monscatter(n,a,index,b)
+C
+       implicit none
+       integer N,INDEX(n),I
+       real A(n),B(n)
+c
+c
+       DO 100 I=1,N
+          A(INDEX(I))=B(I)
+100    CONTINUE
+c
+       return
+       end SUBROUTINE monscatter
+!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+	      SUBROUTINE lwtr(ngrid,coef,lstrong,dup,transm)
+      IMPLICIT NONE
+      INTEGER ngrid
+      REAL coef
+      LOGICAL lstrong
+      REAL dup(ngrid),transm(ngrid)
+      INTEGER ig
+
+      IF(lstrong) THEN
+         DO ig=1,ngrid
+            transm(ig)=exp(-coef*sqrt(dup(ig)))
+         ENDDO
+      ELSE
+         DO ig=1,ngrid
+            transm(ig)=exp(-coef*dup(ig))
+         ENDDO
+      ENDIF
+      RETURN
+      END subroutine lwtr 
+!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+       END MODULE RADIATION 
Index: codes/icosagcm/trunk/src/surface_mod.F
===================================================================
--- codes/icosagcm/trunk/src/surface_mod.F	(revision 149)
+++ codes/icosagcm/trunk/src/surface_mod.F	(revision 149)
@@ -0,0 +1,935 @@
+	MODULE SURFACE_PROCESS
+	
+	USE ICOSA
+	USE dimphys_mod
+	USE RADIATION
+	DATA lmixmin,emin_turb,karman/100.,1.e-8,.4/
+
+	contains
+
+!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+	      subroutiNE vdif(ngrid,nlay,ptime,
+     $                ptimestep,pcapcal,pz0,
+     $                pplay,pplev,pzlay,pzlev,
+     $                pu,pv,ph,ptsrf,pemis,
+     $                pdufi,pdvfi,pdhfi,pfluxsrf,
+     $                pdudif,pdvdif,pdhdif,pdtsrf,pq2,pq2l,
+     $                lwrite)
+      IMPLICIT NONE
+
+c=======================================================================
+c
+c   Diffusion verticale
+c   Shema implicite
+c   On commence par rajouter au variables x la tendance physique
+c   et on resoult en fait:
+c      x(t+1) =  x(t) + dt * (dx/dt)phys(t)  +  dt * (dx/dt)difv(t+1)
+c
+c   !!! attention :
+c   pour utilisation sur une machine sans allocation dynamique de 
+c   memoires (sur SUN par exemple) il faut que ngrid soit egal
+c   a ngrid.
+c
+c   arguments:
+c   ----------
+c
+c   entree:
+c   -------
+c
+c
+c=======================================================================
+
+c-----------------------------------------------------------------------
+c   declarations:
+c   -------------
+c
+c   arguments:
+c   ----------
+
+      INTEGER ngrid,nlay
+      REAL ptime,ptimestep
+      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
+      REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
+      REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
+      REAL ptsrf(ngrid),pemis(ngrid)
+      REAL pdufi(ngrid,nlay),pdvfi(ngrid,nlay),pdhfi(ngrid,nlay)
+      REAL pfluxsrf(ngrid)
+      REAL pdudif(ngrid,nlay),pdvdif(ngrid,nlay),pdhdif(ngrid,nlay)
+      REAL pdtsrf(ngrid),pcapcal(ngrid),pz0(ngrid)
+      REAL pq2(ngrid,nlay+1),pq2l(ngrid,nlay+1)
+      LOGICAL lwrite
+c
+c   local:
+c   ------
+
+      INTEGER ilev,ig,ilay,nlev
+      INTEGER unit,ierr,it1,it2,icount
+      SAVE icount
+      INTEGER cluvdb,putdat,putvdim,setname,setvdim
+      REAL z4st,zdplanck(ngrid),zu2
+      REAL zkv(ngrid,nlayermx+1),zkh(ngrid,nlayermx+1)
+      REAL zcdv(ngrid),zcdh(ngrid)
+      REAL zu(ngrid,nlayermx),zv(ngrid,nlayermx)
+      REAL zh(ngrid,nlayermx)
+      REAL ztsrf2(ngrid)
+      REAL z1(ngrid),z2(ngrid)
+      REAL za(ngrid,nlayermx),zb(ngrid,nlayermx)
+      REAL zb0(ngrid,nlayermx)
+      REAL zc(ngrid,nlayermx),zd(ngrid,nlayermx)
+      REAL zout_dyn(iim+1,jjm+1,nlayermx+1),zout_fi(ngrid,nlayermx+1)
+      REAL zcst1
+      REAL karman
+
+      EXTERNAL coefdifv
+      EXTERNAL SSUM
+      REAL SSUM
+      SAVE karman
+
+      DATA karman/0.4/
+      DATA icount/0/
+c
+c-----------------------------------------------------------------------
+c   initialisations:
+c   ----------------
+
+      nlev=nlay+1
+
+      IF(ngrid.NE.ngrid) THEN
+         PRINT*,'STOP dans coefdifv'
+         PRINT*,'probleme de dimensions :'
+         PRINT*,'ngrid  =',ngrid
+         PRINT*,'ngrid  =',ngrid
+         STOP
+      ENDIF
+
+c   computation of rho*dz and dt*rho/dz=dt*rho**2 g/dp:
+c   with rho=p/RT=p/ (R Theta) (p/ps)**kappa
+c   ---------------------------------
+
+      DO ilay=1,nlay
+         DO ig=1,ngrid
+            za(ig,ilay)=
+     s       (pplev(ig,ilay)-pplev(ig,ilay+1))/g
+         ENDDO
+      ENDDO
+
+      zcst1=4.*g*ptimestep/(kappa*cpp)**2 
+      DO ilev=2,nlev-1
+         DO ig=1,ngrid
+            zb0(ig,ilev)=pplev(ig,ilev)*
+     s      (pplev(ig,1)/pplev(ig,ilev))**kappa /
+     s      (ph(ig,ilev-1)+ph(ig,ilev))
+            zb0(ig,ilev)=zcst1*zb0(ig,ilev)*zb0(ig,ilev)/
+     s      (pplay(ig,ilev-1)-pplay(ig,ilev))
+         ENDDO
+      ENDDO
+      DO ig=1,ngrid
+         zb0(ig,1)=ptimestep*pplev(ig,1)/(kappa*cpp*ptsrf(ig))
+      ENDDO
+      IF(lwrite) THEN
+         ig=ngrid/2+1
+         PRINT*,'Pression (mbar) ,altitude (km),u,v,theta, rho dz'
+         DO ilay=1,nlay
+            WRITE(*,*) .01*pplay(ig,ilay),.001*pzlay(ig,ilay),
+     s      pu(ig,ilay),pv(ig,ilay),ph(ig,ilay),za(ig,ilay)
+         ENDDO
+         PRINT*,'Pression (mbar) ,altitude (km),zb'
+         DO ilev=1,nlay
+            WRITE(*,*) .01*pplev(ig,ilev),.001*pzlev(ig,ilev),
+     s      zb0(ig,ilev)
+         ENDDO
+      ENDIF
+
+c-----------------------------------------------------------------------
+c   2. ajout des tendances physiques:
+c   ------------------------------
+
+      DO ilev=1,nlay
+         DO ig=1,ngrid
+            zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep
+            zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep
+            zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
+         ENDDO
+      ENDDO
+
+c-----------------------------------------------------------------------
+c   3. calcul de  cd :
+c   ----------------
+c
+      CALL vdif_cd( ngrid,pz0,g,pzlay,pu,pv,ptsrf,ph,zcdv,zcdh)
+
+c     CALL my_25(ptimestep,g,pzlev,pzlay,pu,pv,ph,zcdv,
+c    a          pq2,pq2l,zkv,zkh)
+
+        CALL vdif_k(ngrid,nlay,
+     s   ptimestep,g,pzlev,pzlay,pz0,pu,pv,ph,zcdv,zkv,zkh,pq2,pq2l)
+
+      DO ig=1,ngrid
+         zu2=pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1)
+         zcdv(ig)=zcdv(ig)*sqrt(zu2)
+         zcdh(ig)=zcdh(ig)*sqrt(zu2)
+      ENDDO
+
+      IF(lwrite) THEN
+         PRINT*
+         PRINT*,'Diagnostique diffusion verticale'
+         PRINT*,'coefficients Cd pour v et h'
+         PRINT*,zcdv(ngrid/2+1),zcdh(ngrid/2+1)
+         PRINT*,'coefficients K pour v et h'
+         DO ilev=1,nlay
+            PRINT*,zkv(ngrid/2+1,ilev),zkh(ngrid/2+1,ilev)
+         ENDDO
+      ENDIF
+
+c-----------------------------------------------------------------------
+c   integration verticale pour u:
+c   -----------------------------
+c
+      CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2))
+      CALL multipl(ngrid,zcdv,zb0,zb)
+      DO ig=1,ngrid
+         z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
+         zc(ig,nlay)=za(ig,nlay)*zu(ig,nlay)*z1(ig)
+         zd(ig,nlay)=zb(ig,nlay)*z1(ig)
+      ENDDO
+
+      DO ilay=nlay-1,1,-1
+         DO ig=1,ngrid
+            z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
+     $         zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
+            zc(ig,ilay)=(za(ig,ilay)*zu(ig,ilay)+
+     $         zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
+            zd(ig,ilay)=zb(ig,ilay)*z1(ig)
+         ENDDO
+      ENDDO
+
+      DO ig=1,ngrid
+         zu(ig,1)=zc(ig,1)
+      ENDDO
+      DO ilay=2,nlay
+         DO ig=1,ngrid
+            zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1)
+         ENDDO
+      ENDDO
+
+c-----------------------------------------------------------------------
+c   integration verticale pour v:
+c   -----------------------------
+c
+      DO ig=1,ngrid
+         z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
+         zc(ig,nlay)=za(ig,nlay)*zv(ig,nlay)*z1(ig)
+         zd(ig,nlay)=zb(ig,nlay)*z1(ig)
+      ENDDO
+
+      DO ilay=nlay-1,1,-1
+         DO ig=1,ngrid
+            z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
+     $         zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
+            zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay)+
+     $         zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
+            zd(ig,ilay)=zb(ig,ilay)*z1(ig)
+         ENDDO
+      ENDDO
+
+      DO ig=1,ngrid
+         zv(ig,1)=zc(ig,1)
+      ENDDO
+      DO ilay=2,nlay
+         DO ig=1,ngrid
+            zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1)
+         ENDDO
+      ENDDO
+
+c-----------------------------------------------------------------------
+c   integration verticale pour h:
+c   -----------------------------
+c
+      CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
+      CALL multipl(ngrid,zcdh,zb0,zb)
+      DO ig=1,ngrid
+         z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
+         zc(ig,nlay)=za(ig,nlay)*zh(ig,nlay)*z1(ig)
+         zd(ig,nlay)=zb(ig,nlay)*z1(ig)
+      ENDDO
+
+      DO ilay=nlay-1,1,-1
+         DO ig=1,ngrid
+            z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
+     $         zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
+            zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+
+     $         zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
+            zd(ig,ilay)=zb(ig,ilay)*z1(ig)
+         ENDDO
+      ENDDO
+
+c-----------------------------------------------------------------------
+c   rajout eventuel de planck dans le shema implicite:
+c   --------------------------------------------------
+
+      z4st=4.*5.67e-8*ptimestep
+c     z4st=0.
+      DO ig=1,ngrid
+         zdplanck(ig)=z4st*pemis(ig)*ptsrf(ig)*ptsrf(ig)*ptsrf(ig)
+      ENDDO
+
+c-----------------------------------------------------------------------
+c   calcul le l'evolution de la temperature du sol':
+c   -----------------------------------------------
+
+      DO ig=1,ngrid
+         z1(ig)=pcapcal(ig)*ptsrf(ig)+cpp*zb(ig,1)*zc(ig,1)
+     s     +zdplanck(ig)*ptsrf(ig)+ pfluxsrf(ig)*ptimestep
+         z2(ig)= pcapcal(ig)+cpp*zb(ig,1)*(1.-zd(ig,1))+zdplanck(ig)
+         ztsrf2(ig)=z1(ig)/z2(ig)
+         zh(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig)
+         pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep
+      ENDDO
+
+c-----------------------------------------------------------------------
+c   integration verticale finale:
+c   -----------------------------
+
+      DO ilay=2,nlay
+         DO ig=1,ngrid
+            zh(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zh(ig,ilay-1)
+         ENDDO
+      ENDDO
+
+c-----------------------------------------------------------------------
+c   calcul final des tendances de la diffusion verticale:
+c   -----------------------------------------------------
+
+      DO ilev = 1, nlay
+         DO ig=1,ngrid
+            pdudif(ig,ilev)=(    zu(ig,ilev)-
+     $      (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep)    )/ptimestep
+            pdvdif(ig,ilev)=(    zv(ig,ilev)-
+     $      (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep)    )/ptimestep
+            pdhdif(ig,ilev)=(    zh(ig,ilev)-
+     $      (ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep)    )/ptimestep
+         ENDDO
+      ENDDO
+
+      IF(lwrite) THEN
+         PRINT*
+         PRINT*,'Diagnostique de la diffusion verticale'
+         PRINT*,'h avant et apres diffusion verticale'
+         PRINT*,ptsrf(ngrid/2+1),ztsrf2(ngrid/2+1)
+         DO 3110 ilev=1,nlay
+            PRINT*,ph(ngrid/2+1,ilev),zh(ngrid/2+1,ilev)
+3110     CONTINUE
+      ENDIF
+c---------------------------------------------------------------------
+      RETURN
+      END SUBROUTINE vdif
+!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+	      SUBROUTINE convadj(ngrid,nlay,ptimestep,
+     S                   pplay,pplev,ppopsk,
+     $                   pu,pv,ph,
+     $                   pdufi,pdvfi,pdhfi,
+     $                   pduadj,pdvadj,pdhadj)
+      IMPLICIT NONE
+
+c=======================================================================
+c
+c   ajustement convectif sec
+c   on peut ajouter les tendances pdhfi au profil pdh avant l'ajustement
+c'
+c=======================================================================
+
+c-----------------------------------------------------------------------
+c   declarations:
+c   -------------
+c   arguments:
+c   ----------
+
+      INTEGER ngrid,nlay
+      REAL ptimestep
+      REAL ph(ngrid,nlay),pdhfi(ngrid,nlay),pdhadj(ngrid,nlay)
+      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1),ppopsk(ngrid,nlay)
+      REAL pu(ngrid,nlay),pdufi(ngrid,nlay),pduadj(ngrid,nlay)
+      REAL pv(ngrid,nlay),pdvfi(ngrid,nlay),pdvadj(ngrid,nlay)
+
+c   local:
+c   ------
+
+      INTEGER ig,i,l,l1,l2,jj
+      INTEGER jcnt, jadrs(ngrid)
+
+      REAL*8 sig(nlayermx+1),sdsig(nlayermx),dsig(nlayermx)
+      REAL*8 zu(ngrid,nlayermx),zv(ngrid,nlayermx)
+      REAL*8 zh(ngrid,nlayermx)
+      REAL*8 zu2(ngrid,nlayermx),zv2(ngrid,nlayermx)
+      REAL*8 zh2(ngrid,nlayermx)
+      REAL*8 zhm,zsm,zum,zvm,zalpha
+
+      LOGICAL vtest(ngrid),down
+
+c
+c-----------------------------------------------------------------------
+c   initialisation:
+c   ---------------
+c
+      IF(ngrid.NE.ngrid) THEN
+         PRINT*
+         PRINT*,'STOP dans convadj'
+         PRINT*,'ngrid    =',ngrid
+         PRINT*,'ngrid  =',ngrid
+      ENDIF
+c
+c-----------------------------------------------------------------------
+c   detection des profils a modifier:
+c   ---------------------------------
+c   si le profil est a modifier
+c   (i.e. ph(niv_sup) < ph(niv_inf) )
+c   alors le tableau "vtest" est mis a .TRUE. ;
+c   sinon, il reste a sa valeur initiale (.FALSE.)
+c   cette operation est vectorisable
+c   On en profite pour copier la valeur initiale de "ph"
+c   dans le champ de travail "zh"
+
+
+      DO 1010 l=1,nlay
+         DO 1015 ig=1,ngrid
+            zh(ig,l)=ph(ig,l)+pdhfi(ig,l)*ptimestep
+            zu(ig,l)=pu(ig,l)+pdufi(ig,l)*ptimestep
+            zv(ig,l)=pv(ig,l)+pdvfi(ig,l)*ptimestep
+1015     CONTINUE
+1010  CONTINUE
+
+      zu2(:,:)=zu(:,:)
+      zv2(:,:)=zv(:,:)
+      zh2(:,:)=zh(:,:)
+
+      DO 1020 ig=1,ngrid
+        vtest(ig)=.FALSE.
+ 1020 CONTINUE
+c
+      DO 1040 l=2,nlay
+        DO 1060 ig=1,ngrid
+CRAY      vtest(ig)=CVMGM(.TRUE. , vtest(ig),
+CRAY .                      zh2(ig,l)-zh2(ig,l-1))
+          IF(zh2(ig,l).LT.zh2(ig,l-1)) vtest(ig)=.TRUE.
+ 1060   CONTINUE
+ 1040 CONTINUE
+c
+CRAY  CALL WHENNE(ngrid, vtest, 1, 0, jadrs, jcnt)
+      jcnt=0
+      DO 1070 ig=1,ngrid
+         IF(vtest(ig)) THEN
+            jcnt=jcnt+1
+            jadrs(jcnt)=ig
+         ENDIF
+1070  CONTINUE
+
+
+c-----------------------------------------------------------------------
+c  Ajustement des "jcnt" profils instables indices par "jadrs":
+c  ------------------------------------------------------------
+c
+      DO 1080 jj = 1, jcnt
+c
+          i = jadrs(jj)
+c
+c   Calcul des niveaux sigma sur cette colonne
+          DO l=1,nlay+1
+            sig(l)=pplev(i,l)/pplev(i,1)
+         ENDDO
+         DO l=1,nlay
+            dsig(l)=sig(l)-sig(l+1)
+            sdsig(l)=ppopsk(i,l)*dsig(l)
+         ENDDO
+          l2 = 1
+c
+c      -- boucle de sondage vers le haut
+c
+cins$     Loop
+ 8000     CONTINUE
+c
+            l2 = l2 + 1
+c
+cins$       Exit
+            IF (l2 .GT. nlay) Goto 8001
+c
+            IF (zh2(i, l2) .LT. zh2(i, l2-1)) THEN
+c
+c          -- l2 est le niveau le plus haut de la colonne instable
+c
+              l1 = l2 - 1
+              l  = l1
+              zsm = sdsig(l2)
+              zhm = zh2(i, l2)
+c
+c          -- boucle de sondage vers le bas
+c
+cins$         Loop
+ 8020         CONTINUE
+c
+                zsm = zsm + sdsig(l)
+                zhm = zhm + sdsig(l) * (zh2(i, l) - zhm) / zsm
+c
+c            -- doit on etendre la colonne vers le bas ?
+c
+c_EC (M1875) 20/6/87 : AND -> AND THEN
+c
+                down = .FALSE.
+                IF (l1 .NE. 1) THEN    !--  and then
+                  IF (zhm .LT. zh2(i, l1-1)) THEN
+                    down = .TRUE.
+                  END IF
+                END IF
+c
+                IF (down) THEN
+c
+                  l1 = l1 - 1
+                  l  = l1
+c
+                ELSE
+c
+c              -- peut on etendre la colonne vers le haut ?
+c
+cins$             Exit
+                  IF (l2 .EQ. nlay) Goto 8021
+c
+cins$             Exit
+                  IF (zh2(i, l2+1) .GE. zhm) Goto 8021
+c
+                  l2 = l2 + 1
+                  l  = l2
+c
+                END IF
+c
+cins$         End Loop
+              GO TO 8020
+ 8021         CONTINUE
+c
+c          -- nouveau profil : constant (valeur moyenne)
+c
+              zalpha=0.
+              zum=0.
+              zvm=0.
+              DO 1100 l = l1, l2
+                zalpha=zalpha+ABS(zh2(i,l)-zhm)*dsig(l)
+                zh2(i, l) = zhm
+                zum=zum+dsig(l)*zu(i,l)
+                zvm=zvm+dsig(l)*zv(i,l)
+ 1100         CONTINUE
+              zalpha=zalpha/(zhm*(sig(l1)-sig(l2+1)))
+              zum=zum/(sig(l1)-sig(l2+1))
+              zvm=zvm/(sig(l1)-sig(l2+1))
+              IF(zalpha.GT.1.) THEN
+                 PRINT*,'WARNING dans convadj zalpha=',zalpha
+                 if(ig.eq.1) then
+                    print*,'Au pole nord'
+                 elseif (ig.eq.ngrid) then
+                    print*,'Au pole sud'
+                 else
+                    print*,'Point i=',
+     .              ig-((ig-1)/iim)*iim,'j=',(ig-1)/iim+1
+                 endif
+!                 STOP !problem with icosa pole 
+                 zalpha=1.
+              ELSE
+c                IF(zalpha.LT.0.) STOP'zalpha=0'
+                 IF(zalpha.LT.1.e-5) zalpha=1.e-5
+              ENDIF
+              DO l=l1,l2
+                 zu2(i,l)=zu2(i,l)+zalpha*(zum-zu2(i,l))
+                 zv2(i,l)=zv2(i,l)+zalpha*(zvm-zv2(i,l))
+              ENDDO
+ 
+              l2 = l2 + 1
+c
+            END IF
+c
+cins$     End Loop
+          GO TO 8000
+ 8001     CONTINUE
+c
+ 1080 CONTINUE
+c
+      DO 4000 l=1,nlay
+        DO 4020 ig=1,ngrid
+          pdhadj(ig,l)=(zh2(ig,l)-zh(ig,l))/ptimestep
+          pduadj(ig,l)=(zu2(ig,l)-zu(ig,l))/ptimestep
+          pdvadj(ig,l)=(zv2(ig,l)-zv(ig,l))/ptimestep
+ 4020   CONTINUE
+ 4000 CONTINUE
+c
+      RETURN
+      END SUBROUTINE convadj 
+!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+	      SUBROUTINE soil(ngrid,nsoil,firstcall,ptherm_i,
+     s          ptimestep,ptsrf,ptsoil,
+     s          pcapcal,pfluxgrd)
+      IMPLICIT NONE
+
+c=======================================================================
+c
+c   Auteur:  Frederic Hourdin     30/01/92
+c   -------
+c
+c   objet:  computation of : the soil temperature evolution
+c   ------                   the surfacic heat capacity "Capcal"
+c                            the surface conduction flux pcapcal
+c
+c
+c   Method: implicit time integration
+c   -------
+c   Consecutive ground temperatures are related by:
+c           T(k+1) = C(k) + D(k)*T(k)  (1)
+c   the coefficients C and D are computed at the t-dt time-step.
+c   Routine structure:
+c   1)new temperatures are computed  using (1)
+c   2)C and D coefficients are computed from the new temperature
+c     profile for the t+dt time-step
+c   3)the coefficients A and B are computed where the diffusive
+c     fluxes at the t+dt time-step is given by
+c            Fdiff = A + B Ts(t+dt)
+c     or     Fdiff = F0 + Capcal (Ts(t+dt)-Ts(t))/dt
+c            with F0 = A + B (Ts(t))
+c                 Capcal = B*dt
+c           
+c   Interface:
+c   ----------
+c
+c   Arguments:
+c   ----------
+c   ngird               number of grid-points
+c   ptimestep              physical timestep (s)
+c   pto(ngrid,nsoil)     temperature at time-step t (K)
+c   ptn(ngrid,nsoil)     temperature at time step t+dt (K)
+c   pcapcal(ngrid)      specific heat (W*m-2*s*K-1)
+c   pfluxgrd(ngrid)      surface diffusive flux from ground (Wm-2)
+c   
+c=======================================================================
+c   declarations:
+c   -------------
+c-----------------------------------------------------------------------
+c  arguments
+c  ---------
+
+      INTEGER ngrid,nsoil
+      REAL ptimestep
+      REAL ptsrf(ngrid),ptsoil(ngrid,nsoilmx),ptherm_i(ngrid)
+      REAL pcapcal(ngrid),pfluxgrd(ngrid)
+      LOGICAL firstcall
+
+c-----------------------------------------------------------------------
+c  local arrays
+c  ------------
+
+      INTEGER ig,jk
+      REAL za(ngrid),zb(ngrid)
+      REAL zdz2(nsoilmx),z1(ngrid)
+      REAL min_period,dalph_soil
+
+c   local saved variables:
+c   ----------------------
+      REAL dz1(nsoilmx),dz2(nsoilmx)
+      REAL zc(ngrid,nsoilmx),zd(ngrid,nsoilmx)
+      REAL lambda
+
+!!!!!!!! SARVESH !!!!!!! SAVE ATTRIBUTE 
+!!      SAVE dz1,dz2,zc,zd,lambda
+
+c-----------------------------------------------------------------------
+c   Depthts:
+c   --------
+
+      REAL fz,rk,fz1,rk1,rk2
+      fz(rk)=fz1*(dalph_soil**rk-1.)/(dalph_soil-1.)
+
+      IF (firstcall) THEN
+
+c-----------------------------------------------------------------------
+c   ground levels 
+c   grnd=z/l where l is the skin depth of the diurnal cycle:
+c   --------------------------------------------------------
+
+         min_period=20000.
+         dalph_soil=2.
+
+         OPEN(99,file='soil.def',status='old',form='formatted',err=9999)
+         READ(99,*) min_period
+         READ(99,*) dalph_soil
+         PRINT*,'Discretization for the soil model'
+         PRINT*,'First level e-folding depth',min_period,
+     s   '   dalph',dalph_soil
+         CLOSE(99)
+9999     CONTINUE
+
+c   la premiere couche represente un dixieme de cycle diurne
+         fz1=sqrt(min_period/3.14)
+
+         DO jk=1,nsoil
+            rk1=jk
+            rk2=jk-1
+            dz2(jk)=fz(rk1)-fz(rk2)
+         ENDDO
+         DO jk=1,nsoil-1
+            rk1=jk+.5
+            rk2=jk-.5
+            dz1(jk)=1./(fz(rk1)-fz(rk2))
+         ENDDO
+         lambda=fz(.5)*dz1(1)
+         PRINT*,'full layers, intermediate layers (secoonds)'
+         DO jk=1,nsoil
+            rk=jk
+            rk1=jk+.5
+            rk2=jk-.5
+            PRINT*,fz(rk1)*fz(rk2)*3.14,
+     s      fz(rk)*fz(rk)*3.14
+         ENDDO
+
+c   Initialisations:
+c   ----------------
+
+      ELSE
+c-----------------------------------------------------------------------
+c   Computation of the soil temperatures using the Cgrd and Dgrd
+c  coefficient computed at the previous time-step:
+c  -----------------------------------------------
+
+c    surface temperature
+         DO ig=1,ngrid
+            ptsoil(ig,1)=(lambda*zc(ig,1)+ptsrf(ig))/
+     s      (lambda*(1.-zd(ig,1))+1.)
+         ENDDO
+
+c   other temperatures
+         DO jk=1,nsoil-1
+            DO ig=1,ngrid
+               ptsoil(ig,jk+1)=zc(ig,jk)+zd(ig,jk)*ptsoil(ig,jk)
+            ENDDO
+         ENDDO
+
+      ENDIF
+c-----------------------------------------------------------------------
+c   Computation of the Cgrd and Dgrd coefficient for the next step:
+c   ---------------------------------------------------------------
+
+      DO jk=1,nsoil
+         zdz2(jk)=dz2(jk)/ptimestep
+      ENDDO
+
+      DO ig=1,ngrid
+         z1(ig)=zdz2(nsoil)+dz1(nsoil-1)
+         zc(ig,nsoil-1)=zdz2(nsoil)*ptsoil(ig,nsoil)/z1(ig)
+         zd(ig,nsoil-1)=dz1(nsoil-1)/z1(ig)
+      ENDDO
+
+      DO jk=nsoil-1,2,-1
+         DO ig=1,ngrid
+            z1(ig)=1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk)))
+            zc(ig,jk-1)=
+     s      (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk))*z1(ig)
+            zd(ig,jk-1)=dz1(jk-1)*z1(ig)
+         ENDDO
+      ENDDO
+
+c-----------------------------------------------------------------------
+c   computation of the surface diffusive flux from ground and
+c   calorific capacity of the ground:
+c   ---------------------------------
+
+      DO ig=1,ngrid
+         pfluxgrd(ig)=ptherm_i(ig)*dz1(1)*
+     s   (zc(ig,1)+(zd(ig,1)-1.)*ptsoil(ig,1))
+         pcapcal(ig)=ptherm_i(ig)*
+     s   (dz2(1)+ptimestep*(1.-zd(ig,1))*dz1(1))
+         z1(ig)=lambda*(1.-zd(ig,1))+1.
+         pcapcal(ig)=pcapcal(ig)/z1(ig)
+         pfluxgrd(ig)=pfluxgrd(ig)
+     s   +pcapcal(ig)*(ptsoil(ig,1)*z1(ig)-lambda*zc(ig,1)-ptsrf(ig))
+     s   /ptimestep
+      ENDDO
+
+      RETURN
+      END SUBROUTINE SOIL
+!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+	      SUBROUTINE vdif_cd( ngrid,pz0,pg,pz,pu,pv,pts,ph,pcdv,pcdh)
+      IMPLICIT NONE
+c=======================================================================
+c
+c   Subject: computation of the surface drag coefficient using the
+c   -------  approch developed by Loui for ECMWF.
+c
+c   Author: Frederic Hourdin  15 /10 /93
+c   -------
+c
+c   Arguments:
+c   ----------
+c
+c   inputs:
+c   ------
+c     ngrid            size of the horizontal grid
+c     pg               gravity (m s -2)
+c     pz(ngrid)        height of the first atmospheric layer
+c     pu(ngrid)        u component of the wind in that layer
+c     pv(ngrid)        v component of the wind in that layer
+c     pts(ngrid)       surfacte temperature
+c     ph(ngrid)        potential temperature T*(p/ps)^kappa
+c
+c   outputs:
+c   --------
+c     pcdv(ngrid)      Cd for the wind
+c     pcdh(ngrid)      Cd for potential temperature
+c
+c=======================================================================
+c
+c-----------------------------------------------------------------------
+c   Declarations:
+c   -------------
+
+c   Arguments:
+c   ----------
+
+      INTEGER ngrid,nlay
+      REAL pz0(ngrid)
+      REAL pg,pz(ngrid)
+      REAL pu(ngrid),pv(ngrid)
+      REAL pts(ngrid),ph(ngrid)
+      REAL pcdv(ngrid),pcdh(ngrid)
+
+c   Local:
+c   ------
+
+      INTEGER ig
+
+      REAL zu2,z1,zri,zcd0,zz
+
+      REAL karman,b,c,d,c2b,c3bc,c3b,z0,umin2
+      LOGICAL firstcal
+      DATA karman,b,c,d,umin2/.4,5.,5.,5.,1.e-12/
+      DATA firstcal/.true./
+      SAVE b,c,d,karman,c2b,c3bc,c3b,firstcal,umin2
+
+c-----------------------------------------------------------------------
+c   couche de surface:
+c   ------------------
+
+c     DO ig=1,ngrid
+c        zu2=pu(ig)*pu(ig)+pv(ig)*pv(ig)+umin2
+c        pcdv(ig)=pz0(ig)*(1.+sqrt(zu2))
+c        pcdh(ig)=pcdv(ig)
+c     ENDDO
+c     RETURN
+
+      IF (firstcal) THEN
+         c2b=2.*b
+         c3bc=3.*b*c
+         c3b=3.*b
+         firstcal=.false.
+      ENDIF
+
+c!!!! WARNING, verifier la formule originale de Louis!
+      DO ig=1,ngrid
+         zu2=pu(ig)*pu(ig)+pv(ig)*pv(ig)+umin2
+         zri=pg*pz(ig)*(ph(ig)-pts(ig))/(ph(ig)*zu2)
+         z1=1.+pz(ig)/pz0(ig)
+         zcd0=karman/log(z1)
+         zcd0=zcd0*zcd0*sqrt(zu2)
+         IF(zri.LT.0.) THEN
+            z1=b*zri/(1.+c3bc*zcd0*sqrt(-z1*zri))
+            pcdv(ig)=zcd0*(1.-2.*z1)
+            pcdh(ig)=zcd0*(1.-3.*z1)
+         ELSE
+            zz=sqrt(1.+d*zri)
+            pcdv(ig)=zcd0/(1.+c2b*zri/zz)
+            pcdh(ig)=zcd0/(1.+c3b*zri*zz)
+         ENDIF
+      ENDDO
+
+c-----------------------------------------------------------------------
+
+      RETURN
+      END SUBROUTINE vdif_cd
+!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+        SUBROUTINE vdif_k(ngrid,nlay,
+     s   ptimestep,pg,pzlev,pzlay,pz0,pu,pv,ph,pcdv,pkv,pkh,pq2,pq2l)
+
+      IMPLICIT NONE
+
+      INTEGER ngrid,nlay
+
+      REAL ptimestep
+      REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
+      REAL pz0(ngrid)
+      REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
+      REAL pg,pcdv(ngrid)
+      REAL pkv(ngrid,nlay+1),pkh(ngrid,nlay+1)
+	REAL pq2(ngrid,nlay+1),pq2l(ngrid,nlay+1) !!!! SARVESH ADDED to 
+
+      INTEGER ig,il
+      REAL zdu,zdv,zri,zdvodz2,zdz,z1,lmix
+
+      REAL karman
+      SAVE karman
+!      DATA lmixmin,emin_turb,karman/100.,1.e-8,.4/
+!!!!! SARVESH !!!!!!
+!Error: Host associated variable 'lmixmin' may not be in the DATA statement 
+
+!      print*,'LMIXMIN',lmixmin
+      DO ig=1,ngrid
+         pkv(ig,1)=0.
+         pkh(ig,1)=0.
+         pkv(ig,nlay+1)=0.
+         pkh(ig,nlay+1)=0.
+      ENDDO
+c    s     ' zdu,zdv,zdz,zdovdz2,ph(ig,il)+ph(ig,il-1)'
+      DO il=2,nlay
+         DO ig=1,ngrid
+            z1=pzlev(ig,il)+pz0(ig)
+            lmix=karman*z1/(1.+karman*z1/lmixmin)
+c           lmix=lmixmin
+c WARNING test lmix=lmixmin
+            zdu=pu(ig,il)-pu(ig,il-1)
+            zdv=pv(ig,il)-pv(ig,il-1)
+            zdz=pzlay(ig,il)-pzlay(ig,il-1)
+            zdvodz2=(zdu*zdu+zdv*zdv)/(zdz*zdz)
+            IF(zdvodz2.LT.1.e-5) THEN
+                pkv(ig,il)=lmix*sqrt(emin_turb)
+            ELSE
+               zri=2.*pg*(ph(ig,il)-ph(ig,il-1))
+     s         / (zdz* (ph(ig,il)+ph(ig,il-1)) *zdvodz2  )
+               pkv(ig,il)=
+     s         lmix*sqrt(MAX(lmix*lmix*zdvodz2*(1-zri/.4),emin_turb))
+            ENDIF
+            pkh(ig,il)=pkv(ig,il)
+c           IF(ig.EQ.ngrid/2+1) PRINT*,il,lmix,pkv(ig,il),
+c    s      zdu,zdv,zdz,zdvodz2,ph(ig,il)+ph(ig,il-1),
+c    s      lmix*lmix*zdvodz2*(1-zri/.4),emin_turb,zri,ph(ig,il)-ph(ig,il-1),
+c    s      ph(ig,il),ph(ig,il-1)
+         ENDDO
+      ENDDO
+
+      RETURN
+      END SUBROUTINE vdif_k
+!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+	      SUBROUTINE multipl(n,x1,x2,y)
+      IMPLICIT NONE
+c====================================================================
+c
+c   multiplication de deux vecteurs
+c
+c=======================================================================
+c
+      INTEGER n,i
+      REAL x1(n),x2(n),y(n)
+c
+      DO 10 i=1,n
+	 y(i)=x1(i)*x2(i)
+10    CONTINUE
+c
+      RETURN
+      END SUBROUTINE multipl
+!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        END MODULE SURFACE_PROCESS
Index: codes/icosagcm/trunk/src/time.f90
===================================================================
--- codes/icosagcm/trunk/src/time.f90	(revision 148)
+++ codes/icosagcm/trunk/src/time.f90	(revision 149)
@@ -11,6 +11,18 @@
   INTEGER,SAVE    :: itau_out, itau_adv, itau_dissip, itau_physics, itaumax
   
+  INTEGER,SAVE :: day_step,ndays
+  REAL(rstd),SAVE :: jD_ref,jH_ref
+  INTEGER,SAVE :: day_ini,day_end,annee_ref,day_ref
+  REAL(rstd),SAVE::start_time
+  CHARACTER(LEN=255) :: time_style
+  INTEGER,SAVE:: an, mois, jour
+  REAL(rstd),SAVE:: heure
+  CHARACTER (LEN=10):: calend
+
   PUBLIC create_time_counter_header, update_time_counter, close_time_counter, init_time,  &
-         dt, write_period, itau_out, itau_adv, itau_dissip, itau_physics, itaumax
+         dt, write_period, itau_out, itau_adv, itau_dissip, itau_physics, itaumax, & 
+day_step,ndays,jD_ref,jH_ref,day_ini,day_end,annee_ref,day_ref,an, mois, jour,heure, &
+            calend,time_style
+
 
 
@@ -24,5 +36,10 @@
   IMPLICIT NONE
   REAL(rstd) :: run_length
-  
+
+
+     time_style='dcmip'
+   CALL getin('time_style',time_style)
+
+   IF (TRIM(time_style)=='dcmip')  Then
     dt=90.
     CALL getin('dt',dt)
@@ -30,13 +47,4 @@
     itaumax=100
     CALL getin('itaumax',itaumax)
-
-    itau_adv=1
-    CALL getin('itau_adv',itau_adv)
-
-    itau_dissip=1
-    CALL getin('itau_dissip',itau_dissip)
-
-    itau_physics=1
-    CALL getin('itau_physics',itau_physics)
 
     run_length=dt*itaumax
@@ -54,4 +62,15 @@
     itau_out=FLOOR(.5+write_period/dt)
     IF (is_mpi_root) PRINT *, 'Output frequency (scaled) set to ',write_period, ' : itau_out = ',itau_out
+    ENDIF 
+
+    itau_adv=1
+    CALL getin('itau_adv',itau_adv)
+
+    itau_dissip=1
+    CALL getin('itau_dissip',itau_dissip)
+
+    itau_physics=1
+    CALL getin('itau_physics',itau_physics)
+
     
     CALL create_time_counter_header
Index: codes/icosagcm/trunk/src/timeloop_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 148)
+++ codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 149)
@@ -19,4 +19,7 @@
   USE physics_mod
   USE mpipara
+  USE IOIPSL
+  USE maxicosa
+  USE check_conserve_mod
   USE trace
   USE transfert_mod
@@ -43,5 +46,4 @@
   REAL(rstd),POINTER :: dtheta_rhodz(:,:),dtheta_rhodzm1(:,:),dtheta_rhodzm2(:,:)
   REAL(rstd),POINTER :: hflux(:,:),wflux(:,:),hfluxt(:,:),wfluxt(:,:)
-
 !  REAL(rstd) :: dt, run_length
   INTEGER :: ind
@@ -49,4 +51,7 @@
   CHARACTER(len=255) :: scheme_name
   LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time
+  CHARACTER(len=7) :: first 
+  REAL(rstd),SAVE :: jD_cur, jH_cur
+  REAL(rstd),SAVE :: start_time
   LOGICAL, PARAMETER :: check=.FALSE.
 !  INTEGER :: itaumax
@@ -88,4 +93,56 @@
   CALL allocate_field(f_wfluxt,field_t,type_real,llm+1)
 
+!----------------------------------------------------
+  IF (TRIM(time_style)=='lmd')  Then
+
+   day_step=180
+   CALL getin('day_step',day_step)
+
+   ndays=1
+   CALL getin('ndays',ndays)
+
+   dt = daysec/REAL(day_step)
+   itaumax = ndays*day_step
+
+   calend = 'earth_360d'
+   CALL getin('calend', calend)
+
+   day_ini = 0
+   CALL getin('day_ini',day_ini)
+
+   day_end = 0
+   CALL getin('day_end',day_end)
+
+   annee_ref = 1998
+   CALL getin('annee_ref',annee_ref)
+
+   start_time = 0
+   CALL getin('start_time',start_time) 
+
+   write_period=0
+   CALL getin('write_period',write_period)
+      
+   write_period=write_period/scale_factor
+   itau_out=FLOOR(write_period/dt)
+   
+   PRINT *, 'Output frequency (scaled) set to ',write_period, ' : itau_out = ',itau_out 
+
+  mois = 1 ; heure = 0.
+  call ymds2ju(annee_ref, mois, day_ref, heure, jD_ref)
+  jH_ref = jD_ref - int(jD_ref) 
+  jD_ref = int(jD_ref) 
+  
+  CALL ioconf_startdate(INT(jD_ref),jH_ref) 
+  write(*,*)'annee_ref, mois, day_ref, heure, jD_ref'
+  write(*,*)annee_ref, mois, day_ref, heure, jD_ref
+  write(*,*)"ndays,day_step,itaumax,dt======>"
+  write(*,*)ndays,day_step,itaumax,dt
+  call ju2ymds(jD_ref+jH_ref,an, mois, jour, heure)
+  write(*,*)'jD_ref+jH_ref,an, mois, jour, heure'
+  write(*,*)jD_ref+jH_ref,an, mois, jour, heure
+  day_end = day_ini + ndays 
+ END IF 
+!----------------------------------------------------
+
   scheme_name='runge_kutta'
   CALL getin('scheme',scheme_name)
@@ -126,5 +183,4 @@
 !  CALL allocate_field(f_dtheta_rhodzm1,field_t,type_real,llm)
 !  CALL allocate_field(f_dtheta_rhodzm2,field_t,type_real,llm)
-
 !  CALL allocate_field(f_theta,field_t,type_real,llm)
 !  CALL allocate_field(f_dtheta,field_t,type_real,llm)
@@ -135,7 +191,7 @@
   CALL init_advect_tracer
   CALL init_physics
-
-   
-  CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
+!========================================= INITIALIZATION 
+! CALL dynetat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q)
+  CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q)
   CALL writefield("phis",f_phis,once=.TRUE.)
   CALL transfert_request(f_q,req_i1) 
@@ -158,6 +214,7 @@
      CALL compute_rhodz(.FALSE., ps, rhodz)   
   END DO
-  
+
   CALL transfert_request(f_phis,req_i0) 
+  CALL transfert_request(f_phis,req_i1) 
   CALL transfert_request(f_phis,req_i1) 
 
@@ -170,11 +227,11 @@
     ENDIF
     
-    IF (is_mpi_root) PRINT *,"It No :",It,"   t :",dt*It
     IF (mod(it,itau_out)==0 ) THEN
-      CALL writefield("q",f_q)
+!      IF (is_mpi_root) PRINT *,"It No :",It,"   t :",dt*It
       CALL update_time_counter(dt*it)
+      CALL compute_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it)  
     ENDIF
-    
-    CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q)
+
+      CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q)
 
     DO stage=1,nb_stage
@@ -189,11 +246,10 @@
        CASE (mlf)
           CALL  leapfrog_matsuno_scheme(stage)
-          
           !      CASE ('leapfrog')
-          !        CALL leapfrog_scheme
+          !      CALL leapfrog_scheme
           !
           !      CASE ('adam_bashforth')
-          !        CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz)
-          !        CALL adam_bashforth_scheme
+          !      CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz)
+          !      CALL adam_bashforth_scheme
        CASE DEFAULT
           STOP
@@ -208,5 +264,5 @@
     IF(MOD(it+1,itau_adv)==0) THEN
 !       CALL transfert_request(f_wfluxt,req_i1) ! FIXME
-!       CALL transfert_request(f_hfluxt,req_e1) ! FIXME
+!      CALL transfert_request(f_hfluxt,req_e1) ! FIXME
 
        CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz)  ! update q and rhodz after RK step
@@ -227,7 +283,13 @@
          END DO
        ENDIF
-
     END IF
-
+!----------------------------------------------------
+  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)
+!	print*,"Just b4 phys"
+    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)
     ENDDO
@@ -493,5 +555,5 @@
           DO i=ii_begin-dd,ii_end+dd
              ij=(j-1)*iim+i
-             m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g 
+             m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g
              IF(comp) THEN
                 rhodz(ij,l) = m
@@ -508,5 +570,5 @@
           STOP
        ELSE
-!          PRINT *, 'No discrepancy between ps and rhodz detected'
+          PRINT *, 'No discrepancy between ps and rhodz detected'
        END IF
     END IF
