Index: codes/icosagcm/devel/src/dcmip/guided_ncar_mod.f90
===================================================================
--- codes/icosagcm/devel/src/dcmip/guided_ncar_mod.f90	(revision 533)
+++ codes/icosagcm/devel/src/dcmip/guided_ncar_mod.f90	(revision 533)
@@ -0,0 +1,164 @@
+MODULE guided_ncar_mod
+  USE icosa
+  PRIVATE
+  
+  INTEGER,SAVE :: case_wind
+!$OMP THREADPRIVATE(case_wind)
+  
+  REAL(rstd), PARAMETER :: alpha=0.0 ! tilt of solid-body rotation
+  REAL(rstd), PARAMETER :: tau_hadley=daysec, 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
+    IMPLICIT NONE
+    CHARACTER(LEN=255) :: wind
+    wind='deform'
+    CALL getin('dcmip1_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(tt, f_ps, f_theta_rhodz, f_u, f_q)
+  USE icosa
+    IMPLICIT NONE
+    REAL(rstd), INTENT(IN):: tt
+    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 :: ue(:,:)
+    INTEGER :: ind
+    
+    DO ind = 1 , ndomain
+      IF (.NOT. assigned_domain(ind)) CYCLE
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind) 
+      ue = f_u(ind)  
+      CALL wind_profile(tt,ue)
+    END DO    
+
+  END SUBROUTINE guided
+  
+
+  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) :: rr1,rr2,bb,cc,aa,hmx
+    REAL(rstd) :: v1(3),v2(3),ny(3)
+    INTEGER :: i,j,n,l
+    REAL(rstd) :: pitbytau,kk, pr, zr, u0, u1, v0
+
+    pitbytau = pi*tt/tau
+    kk = 10*radius/tau
+    u0 = 2*pi*radius/tau ! for solid-body rotation
+!---------------------------------------------------------
+    DO l = 1,llm 
+       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(pi*tt/tau_hadley) ! 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
+    END DO
+    
+  CONTAINS
+	   
+    SUBROUTINE compute_velocity(x,l,velocity)
+      IMPLICIT NONE
+      REAL(rstd),INTENT(IN)  :: x(3)
+      INTEGER,INTENT(IN)::l
+      REAL(rstd),INTENT(OUT) :: velocity(3)
+      REAL(rstd) :: e_lat(3), e_lon(3)
+      REAL(rstd) :: lon,lat
+      REAL(rstd) :: u,v
+	      
+      CALL xyz2lonlat(x/radius,lon,lat)
+      e_lat(1) = -cos(lon)*sin(lat)
+      e_lat(2) = -sin(lon)*sin(lat)
+      e_lat(3) = cos(lat)
+        
+      e_lon(1) = -sin(lon)
+      e_lon(2) = cos(lon)
+      e_lon(3) = 0
+
+      u = 0.0 ; v = 0.0
+    
+      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
+    
+  END SUBROUTINE  wind_profile
+
+
+END MODULE guided_ncar_mod
Index: codes/icosagcm/devel/src/dcmip/terminator.f90
===================================================================
--- codes/icosagcm/devel/src/dcmip/terminator.f90	(revision 533)
+++ codes/icosagcm/devel/src/dcmip/terminator.f90	(revision 533)
@@ -0,0 +1,124 @@
+!================================================================================================
+! This is the "toy" chemistry module.
+!================================================================================================
+
+MODULE Terminator
+
+implicit none
+private
+save
+
+public :: initial_value_Terminator       ! initialize cl and cl2
+public :: tendency_Terminator             ! interface to tendency computation
+
+!integer, parameter :: 8 = selected_real_kind (12)
+
+real(8), parameter :: cly_constant = 4.e-6_8
+real(8), parameter :: pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164_8
+real(8), parameter :: half_pi = pi*0.5_8
+real(8), parameter :: degrees_to_radians = pi/180.0_8
+real(8), parameter :: k1_lat_center =   20.d0*degrees_to_radians
+real(8), parameter :: k1_lon_center =  300.d0*degrees_to_radians
+
+contains
+
+!===============================================================================
+!  Solar photolysis rate and recombination rate
+!===============================================================================
+
+subroutine k_vals( lat, lon, k1, k2 )
+
+!-----------------------------------------------------------------------
+! Arguments:
+!-----------------------------------------------------------------------
+real(8), intent(in)    :: lat, lon  ! latitude and longitude, radians
+real(8), intent(out)   :: k1, k2    ! reaction rates
+
+k1 = 1.0_8*max(0.d0,sin(lat)*sin(k1_lat_center) + cos(lat)*cos(k1_lat_center)*cos(lon-k1_lon_center))
+k2 = 1._8
+
+return
+
+end subroutine k_vals
+
+!===============================================================================
+!  Tendencies of cl and cl2
+!===============================================================================
+
+subroutine tendency_Terminator( lat, lon, cl, cl2, dt, cl_f, cl2_f )
+
+!-----------------------------------------------------------------------
+! Arguments:
+!-----------------------------------------------------------------------
+
+real(8), intent(in)    :: lat, lon  ! latitude and longitude, degrees
+real(8), intent(in)    :: cl, cl2   ! molar mixing ratio of cl and cl2
+real(8), intent(in)    :: dt        ! size of physics time step
+
+real(8), intent(out)   :: cl_f, cl2_f  ! time rate of change of cl and cl2
+
+!-----------------------------------------------------------------------
+! Local variables
+!-----------------------------------------------------------------------
+
+real(8) :: r, det, expdt, el ! useful algebaic quantities used in the computation
+real(8) :: k1, k2            ! reaction rates
+real(8) :: cly               ! quantity that should be conseved
+
+call k_vals( lat, lon, k1, k2 )
+
+r = k1 / (4._8*k2)
+cly = cl + 2._8* cl2
+
+det = sqrt( r*r + 2._8*r*cly )
+expdt = exp( -4._8*k2*det*dt )
+
+if ( abs(det * k2 * dt) .gt. 1e-16 ) then
+el = (1._8 - expdt) /det /dt
+else
+el = 4._8*k2
+endif
+
+cl_f  = -el * (cl - det + r)*(cl + det + r) / (1._8 + expdt + dt*el*(cl + r))
+cl2_f = -cl_f / 2._8
+
+return
+
+end subroutine tendency_Terminator
+
+!===============================================================================
+!  Compute initial values
+!===============================================================================
+
+subroutine initial_value_Terminator( lat, lon, cl, cl2 )
+
+!-----------------------------------------------------------------------
+! Arguments:
+!-----------------------------------------------------------------------
+
+real(8), intent(in)  :: lat, lon  ! latitude and longitude, degrees
+real(8), intent(out) :: cl, cl2   ! molar mixing ratio of cl and cl2
+
+!-----------------------------------------------------------------------
+! Local variables
+!-----------------------------------------------------------------------
+
+real(8) :: r, det  ! useful algebraic forms
+real(8) :: k1, k2  ! reaction rates
+
+call k_vals( lat, lon, k1, k2 )
+
+r = k1 / (4._8*k2)
+det = sqrt(r*r + 2._8*cly_constant*r)
+
+cl  = (det-r)
+cl2 = cly_constant/2._8 - (det-r)/2._8
+
+return
+
+end subroutine initial_value_Terminator
+
+end module Terminator
+
+
+
