Index: /codes/icosagcm/trunk/run_adv.def
===================================================================
--- /codes/icosagcm/trunk/run_adv.def	(revision 24)
+++ /codes/icosagcm/trunk/run_adv.def	(revision 25)
@@ -1,7 +1,17 @@
+#------------------------------- Mesh ---------------------------------
+
 # Number of subdivision on a main triangle (nbp) : integer (default=40)
 nbp=20
 
-# Number of vertical layer (llm) : integer (default=19)
+# optim_it : mesh optimisation : number of iteration : integer (default=0)
+optim_it=0
+
+# Number of vertical layers (llm) : integer (default=19)
 llm=19
+
+# disvert : vertical discretisation : string (default='std') : std, ncar
+disvert=ncar
+
+#---------------------------------- Misc --------------------------------
 
 # number of tracer (nqtot) : integer (default 1)
@@ -11,12 +21,25 @@
 dt = 180.
 
-# scheme type : string ( default='adam_bashforth') euler, leapfrog, leapfrog_matsuno, adam_bashforth)
+# number of timestep (default 100)
+itaumax = 38400
+
+# output field period : integer (default none)
+write_period=3600
+
+# etat0 : initial state : string (default=jablonowsky06) : 
+# jablonowsky06, academic, ncar
+etat0=ncar
+
+# caldyn : computation type for gcm equation : string (default=gcm) : gcm, adv
+caldyn=adv
+
+#------------------------------ Dynamics --------------------------------
+
+# scheme type : string ( default='adam_bashforth')
+# euler, leapfrog, leapfrog_matsuno, adam_bashforth)
 scheme = euler
 
 # matsuno period : integer ( default=5)
 matsuno_period = 5
-
-# number of timestep (default 100)
-itaumax = 38400
 
 # dissipation time graddiv : real (default=5000)
@@ -25,31 +48,33 @@
 
 # dissipation time nxgradrot (default=5000)
-tetagrot = 1800
-tau_gradrot = 1800
-
+tetagrot=1800
+tau_gradrot=1800
 tau_divgrad=1800
-
-# output field period : integer (default none)
-write_period=3600
-
-
-#NCAR testcase : integer ( default=5)
-ncar_test_case=1
-
-# disvert : vertical discretisation : string (default='std') : std, ncar
-disvert=ncar
-
-# optim_it : mesh optimisation : number of iteration : integer (default=0)
-optim_it=0
-
-# initial_state
-
 
 # guided_type : string (default=none) : none, ncar
 guided_type=ncar
 
-# etat0 : initial state : string (default=jablonowsky06) : jablonowsky06, academic, ncar
-etat0=ncar
+#-------------------- parameters for NCAR test cases ------------------------
 
-# caldyn : computation type for gcm equation : string (default=gcm) : gcm, adv
-caldyn=adv
+# NCAR advection test, initial tracer : string ( default='cos_bell')
+# const, slotted_cyl, cos_bell, dbl_cos_bell_q1, dbl_cos_bell_q2, complement, hadley
+ncar_adv_shape=complement
+
+# NCAR advection test, wind field : string (default='deform') : solid, deform, hadley
+ncar_adv_wind=deform
+
+# ncar_dz : model layer thickness in meters: real (default=400) 
+# used if disvert=ncar
+ncar_dz=400
+
+# ncar_T0 : reference temperature for NCAR test cases : real (default=300) 
+# also used by disvert if disvert=ncar
+ncar_T0=300
+
+# ncar_p0 : reference pressure for NCAR test cases : real (default=1e5) 
+# also used by disvert if disvert=ncar
+ncar_p0=1e5
+
+# ncar_disvert_c : exponent for B(eta) : integer (default=1)
+# used by disvert if disvert=ncar
+ncar_disvert_c=1
Index: /codes/icosagcm/trunk/src/advect_tracer.f90
===================================================================
--- /codes/icosagcm/trunk/src/advect_tracer.f90	(revision 24)
+++ /codes/icosagcm/trunk/src/advect_tracer.f90	(revision 25)
@@ -133,4 +133,5 @@
 
     IF ( iadvtr == iapp_tracvl ) THEN 
+       PRINT *, 'Advection scheme'
        bigt = dt*iapp_tracvl
        DO ind=1,ndomain
Index: /codes/icosagcm/trunk/src/caldyn_sw.f90
===================================================================
--- /codes/icosagcm/trunk/src/caldyn_sw.f90	(revision 24)
+++ /codes/icosagcm/trunk/src/caldyn_sw.f90	(revision 25)
@@ -24,5 +24,5 @@
   TYPE(t_request),POINTER :: req(:) 
   TYPE(t_request),POINTER :: req_u(:)   
-  PUBLIC :: allocate_caldyn,swap_caldyn,init_williamson,caldyn,write_caldyn,Compute_PV,Compute_enstrophy
+  PUBLIC :: allocate_caldyn,swap_caldyn,caldyn,write_caldyn,Compute_PV,Compute_enstrophy
 CONTAINS
 
Index: /codes/icosagcm/trunk/src/disvert_ncar.f90
===================================================================
--- /codes/icosagcm/trunk/src/disvert_ncar.f90	(revision 24)
+++ /codes/icosagcm/trunk/src/disvert_ncar.f90	(revision 25)
@@ -6,5 +6,5 @@
   REAL(rstd), SAVE, ALLOCATABLE,TARGET :: presnivs(:)
 
-CONTAINS 
+CONTAINS
 !=========================================================================
 
@@ -21,5 +21,4 @@
   END SUBROUTINE init_disvert  
 
-
   SUBROUTINE disvert(ap,bp,presnivs)
   USE icosa
@@ -28,37 +27,39 @@
   REAL(rstd),INTENT(OUT) :: bp(:)
   REAL(rstd),INTENT(OUT) :: presnivs(:)
+  ! reads from run.def : ncar_dz, ncar_T0, ncar_p0, ncar_disvert_c
+  INTEGER :: l,cindx
+  REAL(rstd) :: H, eta_top, eta
+
+  ncar_dz=400 ; CALL getin('ncar_dz',ncar_dz);
+
+! SELECT CASE( 
+! pressure profile depends on test case
+! coded here for 1-x (transport)
+
+  ncar_T0=300; CALL getin('ncar_T0',ncar_T0)
+  ncar_p0=1e5; CALL getin('ncar_p0',ncar_p0)
+  cindx=1 ; CALL getin('ncar_disvert_c',cindx)
+
+  H = ncar_T0*cpp*kappa/g ! height scale R.T0/g with R=cpp*kappa
+  eta_top = exp(-llm*ncar_dz/H)
+ 
+  do l = 1,llm+1
+     eta = exp(-(l-1)*ncar_dz/H)
+     bp(l) = ((eta - eta_top)/(1 - eta_top))**cindx
+     ap(l) = ncar_p0 * ( eta - bp(l) )
+  ENDDO
+  bp(1)=1.
+  ap(1)=0.
+  bp(llm+1) = 0
   
-  REAL(rstd) :: sig(llm+1)
-  REAL(rstd) :: sigtop
-  REAL(rstd),PARAMETER:: p0=100000.0 
-  INTEGER :: l,cindx
-  REAL(rstd) :: hdz, ehdz 
+  DO l = 1, llm
+     presnivs(l) = 0.5 *( ap(l)+bp(l)*ncar_p0 + ap(l+1)+bp(l+1)*ncar_p0 )
+  ENDDO
 
+  PRINT *, 'Vertical placement of model levels according to DCMIP Appendix E.3'
+  PRINT *, 'Parameters : ncar_dz=', ncar_dz, '  ncar_p0=',ncar_p0, '  ncar_disvert_c=',cindx
+  PRINT *, 'Isothermal amtosphere with ncar_T0=',ncar_T0 
 
-	 hdz = 400.*g/(300.*287.) 
-	 ehdz = exp(-hdz) 
-
-    do l = 1,llm+1
-	sig(l) = ehdz**(l-1) 
-    end do  
-	
-       sigtop = sig(llm+1) 
-        cindx = 1 
-
-    bp(llm+1) =   0.
-    DO l = 1, llm
-      bp(l) = (sig(l) - sigtop)/(1 - sigtop)
-      bp(l) = bp(l)**cindx
-      ap(l) = p0 * ( sig(l) - bp(l) )
-    ENDDO
-    bp(1)=1.
-    ap(1)=0.
-    ap(llm+1) = p0 * ( sig(llm+1) - bp(llm+1) )
-  
-    DO l = 1, llm
-      presnivs(l) = 0.5 *( ap(l)+bp(l)*p0 + ap(l+1)+bp(l+1)*p0 )
-   ENDDO
-
-  END SUBROUTINE disvert
+END SUBROUTINE disvert
 
 END  MODULE disvert_ncar_mod
Index: /codes/icosagcm/trunk/src/etat0_ncar.f90
===================================================================
--- /codes/icosagcm/trunk/src/etat0_ncar.f90	(revision 24)
+++ /codes/icosagcm/trunk/src/etat0_ncar.f90	(revision 25)
@@ -19,7 +19,4 @@
   REAL(rstd), PARAMETER :: rt=radius*0.5
   REAL(rstd), PARAMETER :: zc=5000.0
-  REAL(rstd), PARAMETER :: ps0=100000.0
-  REAL(rstd), PARAMETER :: T0=300
-  REAL(rstd), PARAMETER :: R_d=287.0   
 
   PUBLIC etat0
@@ -78,57 +75,63 @@
   REAL(rstd) :: X2(3),X1(3)
   INTEGER :: i,j,n,l
-  INTEGER :: testcase
-
-    u = 0.0 ; phis = 0 ; theta_rhodz = 0 ; ps = ps0
-	
-    DO l=1, llm+1
-      pr = ap(l) + bp(l)*ps0
-      zr(l) = -R_d*T0/g*log(pr/ps0) 
-    ENDDO    
-
-    DO l=1, llm 
-      zrl(l) = 0.5*(zr(l) + zr(l+1)) 
-    END DO
-
-    testcase=5
-    CALL getin('ncar_test_case',testcase)
- 
-    SELECT CASE(testcase)
-!--------------------------------------------- SINGLE COSINE BELL
-      CASE(0)
-        CALL cosine_bell_1(q) 
-        
-      CASE(1)
-        CALL cosine_bell_2(q) 
-        
-      CASE(2)
-        CALL cosine_bell_2(q) 
-        DO l=1,llm
-          q(:,l)= 0.9 - 0.8*q(:,l)*q(:,l)
-        END DO  
-        
-      CASE(3)
-        CALL slotted_cylinders(q) 	
-
-      CASE(4) 
-        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(5)  ! hadley like meridional circulation 
-        CALL hadleyq(q) 
-	
-      CASE DEFAULT
-        PRINT*,"no such initial profile"
-        STOP
-        
-      END SELECT
- 
-  CONTAINS
+  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
 
 !======================================================================
Index: /codes/icosagcm/trunk/src/guided_mod.f90
===================================================================
--- /codes/icosagcm/trunk/src/guided_mod.f90	(revision 24)
+++ /codes/icosagcm/trunk/src/guided_mod.f90	(revision 25)
@@ -19,5 +19,5 @@
       
       CASE ('ncar')
-        CALL init_guided_ncar(dt)
+        CALL init_guided_ncar
         
       CASE DEFAULT
@@ -30,9 +30,9 @@
 
   
-  SUBROUTINE guided(it, f_ps, f_theta_rhodz, f_u, f_q)
+  SUBROUTINE guided(tt, f_ps, f_theta_rhodz, f_u, f_q)
   USE icosa
   USE guided_ncar_mod, ONLY : guided_ncar => guided
   IMPLICIT NONE
-    INTEGER, INTENT(IN)   :: it
+    REAL(rstd), INTENT(IN):: tt
     TYPE(t_field),POINTER :: f_ps(:)
     TYPE(t_field),POINTER :: f_phis(:)
@@ -45,5 +45,5 @@
       
       CASE ('ncar')
-        CALL guided_ncar(it, f_ps, f_theta_rhodz, f_u, f_q)
+        CALL guided_ncar(tt, f_ps, f_theta_rhodz, f_u, f_q)
     END SELECT 
   
Index: /codes/icosagcm/trunk/src/guided_ncar_mod.f90
===================================================================
--- /codes/icosagcm/trunk/src/guided_ncar_mod.f90	(revision 24)
+++ /codes/icosagcm/trunk/src/guided_ncar_mod.f90	(revision 25)
@@ -3,26 +3,37 @@
   PRIVATE
   
-  REAL(rstd),SAVE :: dt
-  INTEGER,SAVE    :: test
+  INTEGER,SAVE :: case_wind
   
+  REAL(rstd), PARAMETER :: alpha=0.0 ! tilt of solid-body rotation
+  REAL(rstd), PARAMETER :: tau = 12*daysec ! 12 days               ! see p. 16
+  REAL(rstd), PARAMETER :: w0_deform = 23000*pi/tau, b=0.2, ptop=25494.4  ! see p. 16 
+  REAL(rstd), PARAMETER :: u0_hadley=40.,w0_hadley=0.15 ,ztop= 12000. 
+  INTEGER, PARAMETER    :: K_hadley=5
+
   PUBLIC init_guided, guided
 
 CONTAINS
 
-
-  SUBROUTINE init_guided(dt_in)
-  IMPLICIT NONE
-    REAL(rstd),INTENT(IN) :: dt_in
-    
-    dt=dt_in
-    test=2
-    
+  SUBROUTINE init_guided
+    IMPLICIT NONE
+    CHARACTER(LEN=255) :: wind
+    wind='deform'
+    CALL getin('ncar_adv_wind',wind)
+    SELECT CASE(TRIM(wind))
+    CASE('solid')
+       case_wind=0
+    CASE('deform')
+       case_wind=1
+    CASE('hadley')
+       case_wind=2
+    CASE DEFAULT
+       PRINT*,'Bad selector for variable ncar_adv_wind : <', TRIM(wind),'> options are <solid>, <deform>, <hadley>'
+       END SELECT
   END SUBROUTINE init_guided
-
   
-  SUBROUTINE guided(it, f_ps, f_theta_rhodz, f_u, f_q)
+  SUBROUTINE guided(tt, f_ps, f_theta_rhodz, f_u, f_q)
   USE icosa
     IMPLICIT NONE
-    INTEGER, INTENT(IN)   :: it
+    REAL(rstd), INTENT(IN):: tt
     TYPE(t_field),POINTER :: f_ps(:)
     TYPE(t_field),POINTER :: f_phis(:)
@@ -38,5 +49,5 @@
       CALL swap_geometry(ind) 
       ue = f_u(ind)  
-      CALL wind_profile(it,ue)
+      CALL wind_profile(tt,ue)
     END DO    
 
@@ -44,57 +55,64 @@
   
 
-  SUBROUTINE wind_profile(it,ue)
-  USE icosa
-  IMPLICIT NONE
-    INTEGER,INTENT(IN) :: it 
+  SUBROUTINE wind_profile(tt,ue)
+    USE icosa
+    USE disvert_mod
+    IMPLICIT NONE
+    REAL(rstd),INTENT(IN)  :: tt ! current time 
     REAL(rstd),INTENT(OUT) :: ue(iim*3*jjm,llm)
     REAL(rstd) :: lon, lat
     REAL(rstd) :: nx(3),n_norm,Velocity(3,llm)
-    REAL(rstd), PARAMETER :: lon0=3*pi/2,lat0=0.0
     REAL(rstd) :: rr1,rr2,bb,cc,aa,hmx
     REAL(rstd) :: v1(3),v2(3),ny(3)
     INTEGER :: i,j,n,l
-    REAL(rstd) :: zrl(llm)
+    REAL(rstd) :: pitbytau,kk, pr, zr, u0, u1, v0
 
-    CALL compute_zrl(zrl)
+    pitbytau = pi*tt/tau
+    kk = 10*radius/tau
+    u0 = 2*pi*radius/tau ! for solid-body rotation
 !---------------------------------------------------------
     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 compute_velocity(xyz_e(n+u_right,:),l,velocity(:,l))
-          CALL cross_product2(xyz_v(n+z_rdown,:)/radius,xyz_v(n+z_rup,:)/radius,nx)       
-          ue(n+u_right,l)=1e-10
-          n_norm=sqrt(sum(nx(:)**2))
-          IF (n_norm>1e-30) THEN
-            nx=-nx/n_norm*ne(n,right)
-            ue(n+u_right,l)=sum(nx(:)*velocity(:,l))
-            IF (ABS(ue(n+u_right,l))<1e-100) PRINT *,"ue(n+u_right) ==0",i,j,velocity(:,1)
-          ENDIF
+       pr = presnivs(l)
+       zr = -kappa*cpp*ncar_T0/g*log(pr/ncar_p0)  ! reciprocal of (1) p. 13, isothermal atmosphere
+       u1 = w0_deform*radius/b/ptop*cos(2*pitbytau)*(exp((ptop-pr)/b/ptop)-exp((pr-ncar_p0)/b/ptop))
+       v0 = -radius*w0_hadley*pi/(5.0*ztop)*(ncar_p0/pr)*cos(pi*zr/ztop)*cos(pitbytau) ! for Hadley cell
+
+       DO j=jj_begin-1,jj_end+1
+          DO i=ii_begin-1,ii_end+1
+             n=(j-1)*iim+i
+             CALL compute_velocity(xyz_e(n+u_right,:),l,velocity(:,l))
+             CALL cross_product2(xyz_v(n+z_rdown,:)/radius,xyz_v(n+z_rup,:)/radius,nx)       
+             ue(n+u_right,l)=1e-10
+             n_norm=sqrt(sum(nx(:)**2))
+             IF (n_norm>1e-30) THEN
+                nx=-nx/n_norm*ne(n,right)
+                ue(n+u_right,l)=sum(nx(:)*velocity(:,l))
+                IF (ABS(ue(n+u_right,l))<1e-100) PRINT *,"ue(n+u_right)==0",i,j,velocity(:,1)
+             ENDIF
              
-          CALL compute_velocity(xyz_e(n+u_lup,:),l,velocity(:,l))
-          CALL cross_product2(xyz_v(n+z_up,:)/radius,xyz_v(n+z_lup,:)/radius,nx)
-
-          ue(n+u_lup,l)=1e-10
-          n_norm=sqrt(sum(nx(:)**2))
-          IF (n_norm>1e-30) THEN
-            nx=-nx/n_norm*ne(n,lup)
-            ue(n+u_lup,l)=sum(nx(:)*velocity(:,l))
-	  ENDIF
-       
-          CALL compute_velocity(xyz_e(n+u_ldown,:),l,velocity(:,l))
-          CALL cross_product2(xyz_v(n+z_ldown,:)/radius,xyz_v(n+z_down,:)/radius,nx)
-
-          ue(n+u_ldown,l)=1e-10
-          n_norm=sqrt(sum(nx(:)**2))
-          IF (n_norm>1e-30) THEN
-            nx=-nx/n_norm*ne(n,ldown)
-            ue(n+u_ldown,l)=sum(nx(:)*velocity(:,l))
-            IF (ABS(ue(n+u_ldown,l))<1e-100) PRINT *,"ue(n+u_ldown) ==0",i,j
-          ENDIF        
-        ENDDO
-      ENDDO
+             CALL compute_velocity(xyz_e(n+u_lup,:),l,velocity(:,l))
+             CALL cross_product2(xyz_v(n+z_up,:)/radius,xyz_v(n+z_lup,:)/radius,nx)
+             
+             ue(n+u_lup,l)=1e-10
+             n_norm=sqrt(sum(nx(:)**2))
+             IF (n_norm>1e-30) THEN
+                nx=-nx/n_norm*ne(n,lup)
+                ue(n+u_lup,l)=sum(nx(:)*velocity(:,l))
+             ENDIF
+             
+             CALL compute_velocity(xyz_e(n+u_ldown,:),l,velocity(:,l))
+             CALL cross_product2(xyz_v(n+z_ldown,:)/radius,xyz_v(n+z_down,:)/radius,nx)
+             
+             ue(n+u_ldown,l)=1e-10
+             n_norm=sqrt(sum(nx(:)**2))
+             IF (n_norm>1e-30) THEN
+                nx=-nx/n_norm*ne(n,ldown)
+                ue(n+u_ldown,l)=sum(nx(:)*velocity(:,l))
+                IF (ABS(ue(n+u_ldown,l))<1e-100) PRINT *,"ue(n+u_ldown)==0",i,j
+             ENDIF
+          ENDDO
+       ENDDO
     END DO
- 
+    
   CONTAINS
 	   
@@ -104,23 +122,8 @@
       INTEGER,INTENT(IN)::l
       REAL(rstd),INTENT(OUT) :: velocity(3)
-      REAL(rstd) :: u0
-      REAL(rstd)  :: e_lat(3), e_lon(3)
+      REAL(rstd) :: e_lat(3), e_lon(3)
       REAL(rstd) :: lon,lat
-      REAL(rstd),PARAMETER::alpha=0.0
       REAL(rstd) :: u,v
-!--------------------------------- test 1
-      REAL(rstd),parameter:: nday=12,daysec=86400
-      REAL(rstd),PARAMETER :: tau = nday*daysec
-      REAL(rstd) :: tt, pitbytau,kk
-!---------------------------------- hadley 
-      INTEGER,PARAMETER::K_hadly=5
-      REAL(rstd),PARAMETER::u0_hadly=40.,w0_hadly=0.25 ,ztop= 12000. 
-      REAL(rstd)::coff_hadly
-	
-	
-      tt = it*dt 
-      pitbytau = pi*tt/tau 
-      kk = 10*radius/tau 
-      
+	      
       CALL xyz2lonlat(x/radius,lon,lat)
       e_lat(1) = -cos(lon)*sin(lat)
@@ -133,49 +136,25 @@
 
       u = 0.0 ; v = 0.0
-      u0 = 2*pi*radius/tau
     
-      SELECT CASE(test) 
-        CASE(0)  
-          u=u0*(cos(lat)*cos(alpha)+sin(lat)*sin(alpha)*cos(lon))
-          v=-u0*sin(lon)*sin(alpha)
-        CASE(1) 
-          lon = lon - 2*pitbytau 
-          u= kk*sin(lon)*sin(lon)*sin(2*lat)*cos(pitbytau)+ u0*cos(lat) 
-          v= kk*sin(2*lon)*cos(lat)*cos(pitbytau) 
-        CASE(2)
-          coff_hadly = -radius*w0_hadly*pi/(5.0*ztop)  
-          u= u0_hadly*cos(lat) 
-          v= coff_hadly*cos(lat)*sin(5.*lat)*cos(pi*zrl(l)/ztop)*cos(pitbytau)
-        CASE DEFAULT 
-          PRINT*,"not valid choice of wind" 
-       END SELECT 
-
+      SELECT CASE(case_wind) 
+      CASE(0)  ! Solid-body rotation
+         u=u0*(cos(lat)*cos(alpha)+sin(lat)*sin(alpha)*cos(lon))
+         v=-u0*sin(lon)*sin(alpha)
+      CASE(1)  ! 3D Deformational flow - 
+         lon = lon-2*pitbytau
+         u = kk*sin(lon)*sin(lon)*sin(2*lat)*cos(pitbytau)+ u0*cos(lat)
+         u = u + u1*cos(lon)*cos(lat)**2
+         v = kk*sin(2*lon)*cos(lat)*cos(pitbytau) 
+      CASE(2) ! Hadley-like flow
+         u = u0_hadley*cos(lat)
+         v = v0*cos(lat)*sin(5.*lat)    ! Eq. 37 p. 19
+      CASE DEFAULT 
+         PRINT*,"not valid choice of wind" 
+      END SELECT
+      
       Velocity(:)=(u*e_lon(:)+v*e_lat(:)+1e-50)
-
+      
     END SUBROUTINE compute_velocity
     
-    SUBROUTINE compute_zrl(zrl)
-    USE disvert_mod
-    IMPLICIT NONE
-      REAL(rstd),INTENT(OUT) :: zrl(llm)
-      INTEGER :: l
-      REAL(rstd) :: zr(llm+1)
-      REAl(rstd) :: pr
-      REAL(rstd), PARAMETER :: T0=300
-      REAL(rstd), PARAMETER :: R_d=287.0   
-      REAL(rstd), PARAMETER :: ps0=100000.0
-
-      
-      DO    l    = 1, llm+1
-        pr = ap(l) + bp(l)*ps0
-        zr(l) = -R_d*T0/g*log(pr/ps0) 
-      ENDDO    
-
-      DO l = 1, llm 
-        zrl(l) = 0.5*(zr(l) + zr(l+1)) 
-      END DO
-    
-    END SUBROUTINE compute_zrl
-  
   END SUBROUTINE  wind_profile
 
Index: /codes/icosagcm/trunk/src/icosa_mod.f90
===================================================================
--- /codes/icosagcm/trunk/src/icosa_mod.f90	(revision 24)
+++ /codes/icosagcm/trunk/src/icosa_mod.f90	(revision 25)
@@ -15,3 +15,7 @@
   USE transfert_mod
   
+  ! Variables defined by run.def
+
+  REAL(rstd) :: ncar_dz, ncar_p0, ncar_T0 ! read from run.def by disvert
+
 END MODULE icosa
Index: /codes/icosagcm/trunk/src/timeloop_gcm.f90
===================================================================
--- /codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 24)
+++ /codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 25)
@@ -90,5 +90,5 @@
   CALL init_caldyn(dt)
   CALL init_guided(dt)
-!  CALL init_advect_tracer(dt)
+  CALL init_advect_tracer(dt)
   
   CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
@@ -97,7 +97,7 @@
     PRINT *,"It No :",It
 
-    CALL guided(it,f_ps,f_theta_rhodz,f_u,f_q)
+    CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q)
     CALL caldyn(it,f_phis,f_ps,f_theta_rhodz,f_u, f_dps, f_dtheta_rhodz, f_du)
-!    CALL advect_tracer(f_ps,f_u,f_q)
+    CALL advect_tracer(f_ps,f_u,f_q)
     
     SELECT CASE (TRIM(scheme))
