Index: codes/icosagcm/trunk/src/etat0_dcmip1.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0_dcmip1.f90	(revision 55)
+++ codes/icosagcm/trunk/src/etat0_dcmip1.f90	(revision 55)
@@ -0,0 +1,259 @@
+MODULE etat0_dcmip1_mod
+  USE icosa
+  PRIVATE
+
+  REAL(rstd), SAVE  :: h0=1.
+  REAL(rstd), SAVE  :: lon0=3*pi/2
+  REAL(rstd), SAVE  :: lat0=0.0
+  REAL(rstd), SAVE  :: alpha=0.0
+  REAL(rstd), SAVE  :: R0 
+  REAL(rstd), SAVE  :: lat1=0.
+  REAL(rstd), SAVE  :: lat2=0.
+  REAL(rstd), SAVE  :: lon1=pi/6
+  REAL(rstd), SAVE  :: lon2=-pi/6
+  REAL(rstd), SAVE  :: latc1=0.
+  REAL(rstd), SAVE  :: latc2=0.
+  REAL(rstd), SAVE  :: lonc1=5*pi/6
+  REAL(rstd), SAVE  :: lonc2=7*pi/6
+  REAL(rstd), SAVE  :: zt=1000.0
+  REAL(rstd), SAVE  :: rt
+  REAL(rstd), SAVE  :: zc=5000.0
+
+  PUBLIC etat0
+  
+CONTAINS
+  
+  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q)
+  USE icosa
+  IMPLICIT NONE
+    TYPE(t_field),POINTER :: f_ps(:)
+    TYPE(t_field),POINTER :: f_phis(:)
+    TYPE(t_field),POINTER :: f_theta_rhodz(:)
+    TYPE(t_field),POINTER :: f_u(:)
+    TYPE(t_field),POINTER :: f_q(:)
+  
+    REAL(rstd),POINTER :: ps(:)
+    REAL(rstd),POINTER :: phis(:)
+    REAL(rstd),POINTER :: theta_rhodz(:,:)
+    REAL(rstd),POINTER :: u(:,:)
+    REAL(rstd),POINTER :: q(:,:,:)
+    INTEGER :: ind
+    
+    R0=radius*0.5
+    rt=radius*0.5
+    
+    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)
+      CALL compute_etat0_ncar(ps, phis, theta_rhodz, u, q(:,:,1))
+    ENDDO
+
+  END SUBROUTINE etat0
+  
+  SUBROUTINE compute_etat0_ncar(ps, phis, theta_rhodz, u, q)
+  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) :: q(iim*jjm,llm)  
+
+
+  REAL(rstd) :: qxt1(iim*jjm,llm)
+  REAL(rstd) :: lon, lat
+  REAL(rstd) ::dd1,dd2,dd1t1,dd1t2,dd2t1
+  REAL(rstd) :: pr, zr(llm+1), zrl(llm)
+  REAL(rstd) :: rr1,rr2,bb,cc,aa,hmx
+  REAL(rstd) :: X2(3),X1(3)
+  INTEGER :: i,j,n,l
+  CHARACTER(len=255) :: ncar_adv_shape
+
+  u = 0.0 ; phis = 0 ; theta_rhodz = 0 ; ps = ncar_p0
+  
+  DO l=1, llm+1
+     pr = ap(l) + bp(l)*ncar_p0
+     zr(l) = -kappa*cpp*ncar_T0/g*log(pr/ncar_p0) 
+  ENDDO
+  
+  DO l=1, llm 
+     zrl(l) = 0.5*(zr(l) + zr(l+1)) 
+  END DO
+  
+  ncar_adv_shape='cos_bell'
+  CALL getin('ncar_adv_shape',ncar_adv_shape)
+  
+  SELECT CASE(TRIM(ncar_adv_shape))
+     !--------------------------------------------- SINGLE COSINE BELL
+  CASE('const')
+     q=1
+  CASE('cos_bell')
+     CALL cosine_bell_1(q) 
+     
+  CASE('slotted_cyl')
+     CALL slotted_cylinders(q) 	
+     
+  CASE('dbl_cos_bell_q1')
+     CALL cosine_bell_2(q) 
+     
+  CASE('dbl_cos_bell_q2')
+     CALL cosine_bell_2(q) 
+     DO l=1,llm
+        q(:,l)= 0.9 - 0.8*q(:,l)*q(:,l)
+     END DO
+     
+  CASE('complement')
+     ! tracer such that, in combination with the other tracer fields 
+     ! with weight (3/10), the sum is equal to one
+     CALL cosine_bell_2(qxt1)
+     DO l = 1,llm
+        q(:,l) = 0.9 - 0.8*qxt1(:,l)*qxt1(:,l)  
+     END DO
+     q = q + qxt1 
+     CALL slotted_cylinders(qxt1)
+     q = q + qxt1 
+     q = 1. - q*0.3  
+     
+  CASE('hadley')  ! hadley like meridional circulation 
+     CALL hadleyq(q) 
+     
+  CASE DEFAULT
+     PRINT *, 'Bad selector for variable ncar_adv_shape : <', TRIM(ncar_adv_shape),   &
+          '> options are <const>, <slotted_cyl>, <cos_bell>, <dbl_cos_bell_q1>', &
+          '<dbl_cos_bell_q2>, <complement>, <hadley>'
+     STOP
+     
+  END SELECT
+      
+CONTAINS
+
+!======================================================================
+
+    SUBROUTINE cosine_bell_1(hx)
+    IMPLICIT NONE 
+    REAL(rstd) :: hx(iim*jjm,llm)
+    
+      DO l=1,llm 
+        DO j=jj_begin-1,jj_end+1
+          DO i=ii_begin-1,ii_end+1
+            n=(j-1)*iim+i
+            CALL xyz2lonlat(xyz_i(n,:),lon,lat)
+            CALL dist_lonlat(lon0,lat0,lon,lat,rr1)  ! GC distance from center 
+            rr1 = radius*rr1
+            IF ( rr1 .LT. R0 ) then 
+              hx(n,l)= 0.5*h0*(1+cos(pi*rr1/R0)) 
+	    ELSE
+              hx(n,l)=0.0
+	    END IF
+	  END DO
+        END DO
+      END DO
+   
+    END SUBROUTINE cosine_bell_1
+
+!==============================================================================
+
+    SUBROUTINE   cosine_bell_2(hx) 
+    IMPLICIT NONE 	
+    REAL(rstd) :: hx(iim*jjm,llm)
+   
+      DO l=1,llm
+        DO j=jj_begin-1,jj_end+1
+          DO i=ii_begin-1,ii_end+1
+            n=(j-1)*iim+i
+            CALL xyz2lonlat(xyz_i(n,:),lon,lat)
+            CALL dist_lonlat(lonc1,latc1,lon,lat,rr1)  ! GC distance from center 
+            rr1 = radius*rr1
+            CALL dist_lonlat(lonc2,latc2,lon,lat,rr2)  ! GC distance from center 
+            rr2 = radius*rr2
+            dd1t1 = rr1/rt
+            dd1t1 = dd1t1*dd1t1
+            dd1t2 = (zrl(l) - zc)/zt
+            dd1t2 = dd1t2*dd1t2 
+            dd1 = dd1t1 + dd1t2 
+            dd1 = Min(1.0,dd1)  
+            dd2t1 = rr2/rt
+            dd2t1 = dd2t1*dd2t1
+            dd2 = dd2t1 + dd1t2 
+            dd2 = Min(1.0,dd2) 
+            hx(n,l)= 0.5*(1. + cos(pi*dd1))+0.5*(1.+cos(pi*dd2))
+          END DO
+        END DO
+      END DO
+
+    END SUBROUTINE cosine_bell_2
+
+!=============================================================================
+
+    SUBROUTINE slotted_cylinders(hx) 
+    IMPLICIT NONE
+    REAL(rstd) :: hx(iim*jjm,llm)
+    
+      DO l=1,llm
+        DO j=jj_begin-1,jj_end+1
+          DO i=ii_begin-1,ii_end+1
+            n=(j-1)*iim+i
+            CALL xyz2lonlat(xyz_i(n,:),lon,lat)
+            CALL dist_lonlat(lonc1,latc1,lon,lat,rr1)  ! GC distance from center 
+            rr1 = radius*rr1
+            CALL dist_lonlat(lonc2,latc2,lon,lat,rr2)  ! GC distance from center 
+            rr2 = radius*rr2
+            dd1t1 = rr1/rt
+            dd1t1 = dd1t1*dd1t1
+            dd1t2 = (zrl(l) - zc)/zt
+            dd1t2 = dd1t2*dd1t2 
+            dd1 = dd1t1 + dd1t2  
+            dd2t1 = rr2/rt
+            dd2t1 = dd2t1*dd2t1
+            dd2 = dd2t1 + dd1t2
+            
+            IF ( dd1 .LT. 0.5 ) Then 
+              hx(n,l) = 1.0
+            ELSEIF ( dd2 .LT. 0.5 ) Then 
+              hx(n,l) = 1.0
+            ELSE
+	      hx(n,l) = 0.1 
+            END IF  
+    
+            IF ( zrl(l) .GT. zc ) Then 
+              IF ( ABS(latc1 - lat) .LT. 0.125 ) Then 
+                hx(n,l)= 0.0
+ 	      ENDIF 
+            ENDIF  
+  
+          ENDDO
+       END DO
+      END DO
+
+    END SUBROUTINE slotted_cylinders 
+
+!==============================================================================
+
+    SUBROUTINE hadleyq(hx) 
+    IMPLICIT NONE 
+    REAL(rstd)::hx(iim*jjm,llm) 
+    REAL(rstd),PARAMETER:: zz1=3500.,zz2=6500.,zz0=0.5*(zz1+zz2)
+   
+      DO l=1,llm
+        IF ( ( zz1 .LT. zrl(l) ) .and. ( zrl(l) .LT. zz2 ) )  THEN 
+          hx(:,l) = 0.5*(1. + cos(0.002*pi*(zrl(l)-zz0)/3.))  
+        ELSE
+	  hx(:,l) = 0.0 
+        END IF
+      END DO
+    END SUBROUTINE hadleyq 
+
+  END SUBROUTINE compute_etat0_ncar
+  
+
+END MODULE etat0_dcmip1_mod
Index: codes/icosagcm/trunk/src/etat0_dcmip3.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0_dcmip3.f90	(revision 54)
+++ codes/icosagcm/trunk/src/etat0_dcmip3.f90	(revision 55)
@@ -10,21 +10,7 @@
 
   PRIVATE
-! some of the following should be obtained via getin
-!  REAL(rstd),PARAMETER :: ztop=10000.   ! Height position of the model top (m)
-!  REAL(rstd),PARAMETER :: ptop=27391.9  ! Pressure at the model top (Pa)
-!  REAL(rstd),PARAMETER :: Xscale=1.     ! Small-planet scaling factor (1=Earth)
-!  REAL(rstd),PARAMETER :: Omega =0      ! Planetary rotation rate (s-1)
-  REAL(rstd),PARAMETER :: u0=20.         ! Maximum amplitude of the zonal wind (m.s-1)
-  REAL(rstd),PARAMETER :: zs=0.          ! Surface elevation (m)
-  REAL(rstd),PARAMETER :: N=0.01         ! Brunt-Vaisala frequency (s-1)
-  REAL(rstd),PARAMETER :: Teq=300.       ! Surface temperature at the equator (K)
-  REAL(rstd),PARAMETER :: Peq=1e5        ! Reference surface pressure at the equator (hPa)
-  REAL(rstd),PARAMETER :: d=5000.        ! Witdth parameter for theta
-  REAL(rstd),PARAMETER :: lonc=2*pi/3    ! Longitudinal centerpoint of theta
-  REAL(rstd),PARAMETER :: latc=0         ! Longitudinal centerpoint of theta
-  REAL(rstd),PARAMETER :: dtheta=1.      ! Maximum amplitude of theta (K)
-  REAL(rstd),PARAMETER :: Lz=20000.      ! Vertical wave lenght of the theta perturbation
 
   PUBLIC  etat0
+
 CONTAINS
 
@@ -68,4 +54,14 @@
   USE wind_mod
   IMPLICIT NONE
+  REAL(rstd),PARAMETER :: u0=20.         ! Maximum amplitude of the zonal wind (m.s-1)
+  REAL(rstd),PARAMETER :: N=0.01         ! Brunt-Vaisala frequency (s-1)
+  REAL(rstd),PARAMETER :: Teq=300.       ! Surface temperature at the equator (K)
+  REAL(rstd),PARAMETER :: Peq=1e5        ! Reference surface pressure at the equator (hPa)
+  REAL(rstd),PARAMETER :: d=5000.        ! Witdth parameter for theta
+  REAL(rstd),PARAMETER :: lonc=2*pi/3    ! Longitudinal centerpoint of theta
+  REAL(rstd),PARAMETER :: latc=0         ! Longitudinal centerpoint of theta
+  REAL(rstd),PARAMETER :: dtheta=1.      ! Maximum amplitude of theta (K)
+  REAL(rstd),PARAMETER :: Lz=20000.      ! Vertical wave lenght of the theta perturbation
+
   REAL(rstd), INTENT(OUT) :: ps(iim*jjm)
   REAL(rstd), INTENT(OUT) :: phis(iim*jjm)
@@ -85,5 +81,5 @@
   INTEGER :: i,j,l,ij
   REAL(rstd) :: Rd        ! gas constant of dry air, P=rho.Rd.T
-  REAL(rstd) :: alpha, rm, zs
+  REAL(rstd) :: alpha, rm
   REAL(rstd) :: lon,lat, C0, C1, GG
   REAL(rstd) :: p0psk, pspsk,r,zz, thetab, thetap
@@ -156,5 +152,4 @@
      CALL compute_temperature2theta_rhodz(ps,T,theta_rhodz,0)
   ELSE
-!     CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,0)
      CALL compute_temperature2theta_rhodz(ps,T,theta_rhodz,0)
   END IF
Index: codes/icosagcm/trunk/src/etat0_ncar.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0_ncar.f90	(revision 54)
+++ 	(revision )
@@ -1,259 +1,0 @@
-MODULE etat0_ncar_mod
-  USE icosa
-  PRIVATE
-
-  REAL(rstd), SAVE  :: h0=1.
-  REAL(rstd), SAVE  :: lon0=3*pi/2
-  REAL(rstd), SAVE  :: lat0=0.0
-  REAL(rstd), SAVE  :: alpha=0.0
-  REAL(rstd), SAVE  :: R0 
-  REAL(rstd), SAVE  :: lat1=0.
-  REAL(rstd), SAVE  :: lat2=0.
-  REAL(rstd), SAVE  :: lon1=pi/6
-  REAL(rstd), SAVE  :: lon2=-pi/6
-  REAL(rstd), SAVE  :: latc1=0.
-  REAL(rstd), SAVE  :: latc2=0.
-  REAL(rstd), SAVE  :: lonc1=5*pi/6
-  REAL(rstd), SAVE  :: lonc2=7*pi/6
-  REAL(rstd), SAVE  :: zt=1000.0
-  REAL(rstd), SAVE  :: rt
-  REAL(rstd), SAVE  :: zc=5000.0
-
-  PUBLIC etat0
-  
-CONTAINS
-  
-  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q)
-  USE icosa
-  IMPLICIT NONE
-    TYPE(t_field),POINTER :: f_ps(:)
-    TYPE(t_field),POINTER :: f_phis(:)
-    TYPE(t_field),POINTER :: f_theta_rhodz(:)
-    TYPE(t_field),POINTER :: f_u(:)
-    TYPE(t_field),POINTER :: f_q(:)
-  
-    REAL(rstd),POINTER :: ps(:)
-    REAL(rstd),POINTER :: phis(:)
-    REAL(rstd),POINTER :: theta_rhodz(:,:)
-    REAL(rstd),POINTER :: u(:,:)
-    REAL(rstd),POINTER :: q(:,:,:)
-    INTEGER :: ind
-    
-    R0=radius*0.5
-    rt=radius*0.5
-    
-    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)
-      CALL compute_etat0_ncar(ps, phis, theta_rhodz, u, q(:,:,1))
-    ENDDO
-
-  END SUBROUTINE etat0
-  
-  SUBROUTINE compute_etat0_ncar(ps, phis, theta_rhodz, u, q)
-  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) :: q(iim*jjm,llm)  
-
-
-  REAL(rstd) :: qxt1(iim*jjm,llm)
-  REAL(rstd) :: lon, lat
-  REAL(rstd) ::dd1,dd2,dd1t1,dd1t2,dd2t1
-  REAL(rstd) :: pr, zr(llm+1), zrl(llm)
-  REAL(rstd) :: rr1,rr2,bb,cc,aa,hmx
-  REAL(rstd) :: X2(3),X1(3)
-  INTEGER :: i,j,n,l
-  CHARACTER(len=255) :: ncar_adv_shape
-
-  u = 0.0 ; phis = 0 ; theta_rhodz = 0 ; ps = ncar_p0
-  
-  DO l=1, llm+1
-     pr = ap(l) + bp(l)*ncar_p0
-     zr(l) = -kappa*cpp*ncar_T0/g*log(pr/ncar_p0) 
-  ENDDO
-  
-  DO l=1, llm 
-     zrl(l) = 0.5*(zr(l) + zr(l+1)) 
-  END DO
-  
-  ncar_adv_shape='cos_bell'
-  CALL getin('ncar_adv_shape',ncar_adv_shape)
-  
-  SELECT CASE(TRIM(ncar_adv_shape))
-     !--------------------------------------------- SINGLE COSINE BELL
-  CASE('const')
-     q=1
-  CASE('cos_bell')
-     CALL cosine_bell_1(q) 
-     
-  CASE('slotted_cyl')
-     CALL slotted_cylinders(q) 	
-     
-  CASE('dbl_cos_bell_q1')
-     CALL cosine_bell_2(q) 
-     
-  CASE('dbl_cos_bell_q2')
-     CALL cosine_bell_2(q) 
-     DO l=1,llm
-        q(:,l)= 0.9 - 0.8*q(:,l)*q(:,l)
-     END DO
-     
-  CASE('complement')
-     ! tracer such that, in combination with the other tracer fields 
-     ! with weight (3/10), the sum is equal to one
-     CALL cosine_bell_2(qxt1)
-     DO l = 1,llm
-        q(:,l) = 0.9 - 0.8*qxt1(:,l)*qxt1(:,l)  
-     END DO
-     q = q + qxt1 
-     CALL slotted_cylinders(qxt1)
-     q = q + qxt1 
-     q = 1. - q*0.3  
-     
-  CASE('hadley')  ! hadley like meridional circulation 
-     CALL hadleyq(q) 
-     
-  CASE DEFAULT
-     PRINT *, 'Bad selector for variable ncar_adv_shape : <', TRIM(ncar_adv_shape),   &
-          '> options are <const>, <slotted_cyl>, <cos_bell>, <dbl_cos_bell_q1>', &
-          '<dbl_cos_bell_q2>, <complement>, <hadley>'
-     STOP
-     
-  END SELECT
-      
-CONTAINS
-
-!======================================================================
-
-    SUBROUTINE cosine_bell_1(hx)
-    IMPLICIT NONE 
-    REAL(rstd) :: hx(iim*jjm,llm)
-    
-      DO l=1,llm 
-        DO j=jj_begin-1,jj_end+1
-          DO i=ii_begin-1,ii_end+1
-            n=(j-1)*iim+i
-            CALL xyz2lonlat(xyz_i(n,:),lon,lat)
-            CALL dist_lonlat(lon0,lat0,lon,lat,rr1)  ! GC distance from center 
-            rr1 = radius*rr1
-            IF ( rr1 .LT. R0 ) then 
-              hx(n,l)= 0.5*h0*(1+cos(pi*rr1/R0)) 
-	    ELSE
-              hx(n,l)=0.0
-	    END IF
-	  END DO
-        END DO
-      END DO
-   
-    END SUBROUTINE cosine_bell_1
-
-!==============================================================================
-
-    SUBROUTINE   cosine_bell_2(hx) 
-    IMPLICIT NONE 	
-    REAL(rstd) :: hx(iim*jjm,llm)
-   
-      DO l=1,llm
-        DO j=jj_begin-1,jj_end+1
-          DO i=ii_begin-1,ii_end+1
-            n=(j-1)*iim+i
-            CALL xyz2lonlat(xyz_i(n,:),lon,lat)
-            CALL dist_lonlat(lonc1,latc1,lon,lat,rr1)  ! GC distance from center 
-            rr1 = radius*rr1
-            CALL dist_lonlat(lonc2,latc2,lon,lat,rr2)  ! GC distance from center 
-            rr2 = radius*rr2
-            dd1t1 = rr1/rt
-            dd1t1 = dd1t1*dd1t1
-            dd1t2 = (zrl(l) - zc)/zt
-            dd1t2 = dd1t2*dd1t2 
-            dd1 = dd1t1 + dd1t2 
-            dd1 = Min(1.0,dd1)  
-            dd2t1 = rr2/rt
-            dd2t1 = dd2t1*dd2t1
-            dd2 = dd2t1 + dd1t2 
-            dd2 = Min(1.0,dd2) 
-            hx(n,l)= 0.5*(1. + cos(pi*dd1))+0.5*(1.+cos(pi*dd2))
-          END DO
-        END DO
-      END DO
-
-    END SUBROUTINE cosine_bell_2
-
-!=============================================================================
-
-    SUBROUTINE slotted_cylinders(hx) 
-    IMPLICIT NONE
-    REAL(rstd) :: hx(iim*jjm,llm)
-    
-      DO l=1,llm
-        DO j=jj_begin-1,jj_end+1
-          DO i=ii_begin-1,ii_end+1
-            n=(j-1)*iim+i
-            CALL xyz2lonlat(xyz_i(n,:),lon,lat)
-            CALL dist_lonlat(lonc1,latc1,lon,lat,rr1)  ! GC distance from center 
-            rr1 = radius*rr1
-            CALL dist_lonlat(lonc2,latc2,lon,lat,rr2)  ! GC distance from center 
-            rr2 = radius*rr2
-            dd1t1 = rr1/rt
-            dd1t1 = dd1t1*dd1t1
-            dd1t2 = (zrl(l) - zc)/zt
-            dd1t2 = dd1t2*dd1t2 
-            dd1 = dd1t1 + dd1t2  
-            dd2t1 = rr2/rt
-            dd2t1 = dd2t1*dd2t1
-            dd2 = dd2t1 + dd1t2
-            
-            IF ( dd1 .LT. 0.5 ) Then 
-              hx(n,l) = 1.0
-            ELSEIF ( dd2 .LT. 0.5 ) Then 
-              hx(n,l) = 1.0
-            ELSE
-	      hx(n,l) = 0.1 
-            END IF  
-    
-            IF ( zrl(l) .GT. zc ) Then 
-              IF ( ABS(latc1 - lat) .LT. 0.125 ) Then 
-                hx(n,l)= 0.0
- 	      ENDIF 
-            ENDIF  
-  
-          ENDDO
-       END DO
-      END DO
-
-    END SUBROUTINE slotted_cylinders 
-
-!==============================================================================
-
-    SUBROUTINE hadleyq(hx) 
-    IMPLICIT NONE 
-    REAL(rstd)::hx(iim*jjm,llm) 
-    REAL(rstd),PARAMETER:: zz1=3500.,zz2=6500.,zz0=0.5*(zz1+zz2)
-   
-      DO l=1,llm
-        IF ( ( zz1 .LT. zrl(l) ) .and. ( zrl(l) .LT. zz2 ) )  THEN 
-          hx(:,l) = 0.5*(1. + cos(0.002*pi*(zrl(l)-zz0)/3.))  
-        ELSE
-	  hx(:,l) = 0.0 
-        END IF
-      END DO
-    END SUBROUTINE hadleyq 
-
-  END SUBROUTINE compute_etat0_ncar
-  
-
-END MODULE etat0_ncar_mod
