Index: codes/icosagcm/trunk/src/dcmip2016_cyclone.f90
===================================================================
--- codes/icosagcm/trunk/src/dcmip2016_cyclone.f90	(revision 414)
+++ codes/icosagcm/trunk/src/dcmip2016_cyclone.f90	(revision 414)
@@ -0,0 +1,269 @@
+MODULE dcmip2016_cyclone_mod
+
+!=======================================================================
+!
+!  Date:  July 29, 2015
+!
+!  Function for setting up idealized tropical cyclone initial conditions
+!
+!  SUBROUTINE tropical_cyclone_sample(
+!    lon,lat,p,z,zcoords,u,v,t,thetav,phis,ps,rho,q)
+!
+!  Given a point specified by: 
+!      lon    longitude (radians) 
+!      lat    latitude (radians) 
+!      p/z    pressure (Pa) / height (m)
+!  zcoords    1 if z is specified, 0 if p is specified
+!
+!  the functions will return:
+!        p    pressure if z is specified (Pa)
+!        z    geopotential height if p is specified (m)
+!        u    zonal wind (m s^-1)
+!        v    meridional wind (m s^-1)
+!        t    temperature (K)
+!   thetav    virtual potential temperature (K)
+!     phis    surface geopotential (m^2 s^-2)
+!       ps    surface pressure (Pa)
+!      rho    density (kj m^-3)
+!        q    specific humidity (kg/kg)
+!
+!  Initial data are currently identical to:
+!
+!       Reed, K. A., and C. Jablonowski, 2011: An analytic
+!       vortex initialization technique for idealized tropical
+!       cyclone studies in AGCMs. Mon. Wea. Rev., 139, 689-710. 
+!
+!  Author: Kevin A. Reed
+!          Stony Brook University
+!          Email: kevin.a.reed@stonybrook.edu
+!
+!=======================================================================
+
+  IMPLICIT NONE
+
+!=======================================================================
+!    Physical constants
+!=======================================================================
+
+  REAL(8), PARAMETER ::               &
+       a     = 6371220.0d0,           & ! Reference Earth's Radius (m)
+       Rd    = 287.0d0,               & ! Ideal gas const dry air (J kg^-1 K^1)
+       g     = 9.80616d0,             & ! Gravity (m s^2)
+       cp    = 1004.5d0,              & ! Specific heat capacity (J kg^-1 K^1)
+       Lvap  = 2.5d6,                 & ! Latent heat of vaporization of water
+       Rvap  = 461.5d0,               & ! Ideal gas constnat for water vapor
+       Mvap  = 0.608d0,               & ! Ratio of molar mass of dry air/water
+       pi    = 3.14159265358979d0,    & ! pi
+       p0    = 100000.0d0,            & ! surface pressure (Pa)
+       kappa = 2.d0/7.d0,             & ! Ratio of Rd to cp
+       omega = 7.29212d-5,            & ! Reference rotation rate of the Earth (s^-1)
+       deg2rad  = pi/180.d0             ! Conversion factor of degrees to radians
+
+!=======================================================================
+!    Test case parameters
+!=======================================================================
+  REAL(8), PARAMETER ::         &
+       rp         = 282000.d0,  & ! Radius for calculation of PS
+       dp         = 1115.d0,    & ! Delta P for calculation of PS
+       zp         = 7000.d0,    & ! Height for calculation of P
+       q0         = 0.021d0,    & ! q at surface from Jordan
+       gamma      = 0.007d0,    & ! lapse rate
+       Ts0        = 302.15d0,   & ! Surface temperature (SST)
+       p00        = 101500.d0,  & ! global mean surface pressure
+       cen_lat    = 10.d0,      & ! Center latitude of initial vortex
+       cen_lon    = 180.d0,     & ! Center longitufe of initial vortex
+       zq1        = 3000.d0,    & ! Height 1 for q calculation
+       zq2        = 8000.d0,    & ! Height 2 for q calculation
+       exppr      = 1.5d0,      & ! Exponent for r dependence of p
+       exppz      = 2.d0,       & ! Exponent for z dependence of p
+       ztrop      = 15000.d0,   & ! Tropopause Height
+       qtrop      = 1.d-11,     & ! Tropopause specific humidity
+       rfpi       = 1000000.d0, & ! Radius within which to use fixed-point iter.
+       constTv    = 0.608d0,    & ! Constant for Virtual Temp Conversion
+       deltaz     = 2.d-13,     & ! Small number to ensure convergence in FPI
+       epsilon    = 1.d-25,     & ! Small number to aviod dividing by zero in wind calc
+       exponent = Rd*gamma/g,   & ! exponent
+       T0    = Ts0*(1.d0+constTv*q0),             & ! Surface temp
+       Ttrop = T0 - gamma*ztrop,                  & ! Tropopause temp
+       ptrop = p00*(Ttrop/T0)**(1.d0/exponent)      ! Tropopause pressure
+
+CONTAINS 
+
+!=======================================================================
+!    Evaluate the tropical cyclone initial conditions
+!=======================================================================
+  SUBROUTINE tropical_cyclone_test(lon,lat,p,z,zcoords,u,v,t,thetav,phis,ps,rho,q) &
+    BIND(c, name = "tropical_cyclone_test")
+
+    IMPLICIT NONE
+
+    !------------------------------------------------
+    !   Input / output parameters
+    !------------------------------------------------
+
+    REAL(8), INTENT(IN) ::     &
+              lon,             &     ! Longitude (radians)
+              lat                    ! Latitude (radians)
+
+    REAL(8), INTENT(INOUT) ::  &
+              p,               &     ! Pressure (Pa)
+              z                      ! Height (m)
+
+    INTEGER, INTENT(IN) :: zcoords     ! 1 if z coordinates are specified
+                                     ! 0 if p coordinates are specified
+
+    REAL(8), INTENT(OUT) ::    &
+              u,               &     ! Zonal wind (m s^-1)
+              v,               &     ! Meridional wind (m s^-1)
+              t,               &     ! Temperature (K)
+              thetav,          &     ! Virtual potential temperature (K)
+              phis,            &     ! Surface Geopotential (m^2 s^-2)
+              ps,              &     ! Surface Pressure (Pa)
+              rho,             &     ! Density (kg m^-3)
+              q                      ! Specific Humidity (kg/kg)
+
+    !------------------------------------------------
+    !   Local variables
+    !------------------------------------------------
+    real(8)  :: d1, d2, d, vfac, ufac, height, zhere, gr, f, zn
+
+    integer  n
+
+    !------------------------------------------------
+    !   Define Great circle distance (gr) and
+    !   Coriolis parameter (f)
+    !------------------------------------------------
+    f  = 2.d0*omega*sin(cen_lat*deg2rad)           ! Coriolis parameter
+    gr = a*acos(sin(cen_lat*deg2rad)*sin(lat) + &  ! Great circle radius
+         (cos(cen_lat*deg2rad)*cos(lat)*cos(lon-cen_lon*deg2rad)))
+
+    !------------------------------------------------
+    !   Initialize PS (surface pressure)
+    !------------------------------------------------
+    ps = p00-dp*exp(-(gr/rp)**exppr) 
+
+    !------------------------------------------------
+    !   Initialize altitude (z) if pressure provided
+    !   or pressure if altitude (z) is provided
+    !------------------------------------------------
+    if (zcoords .eq. 1) then
+
+       height = z
+ 
+       if (height > ztrop) then
+          p = ptrop*exp(-(g*(height-ztrop))/(Rd*Ttrop))
+       else
+          p = (p00-dp*exp(-(gr/rp)**exppr)*exp(-(height/zp)**exppz)) &
+              * ((T0-gamma*height)/T0)**(1/exponent)
+       end if
+ 
+    else
+
+       height = (T0/gamma)*(1.d0-(p/ps)**exponent)
+
+       ! If inside a certain distance of the center of the storm
+       ! perform a Fixed-point iteration to calculate the height
+       ! more accurately
+
+       if (gr < rfpi ) then
+          zhere = height 
+          n = 1
+          20 continue
+          n = n+1
+          zn = zhere - fpiF(p,gr,zhere)/fpidFdz(gr,zhere)
+          if (n.gt.20) then
+              PRINT *,'FPI did not converge after 20 interations in q & T!!!'
+          else if ( abs(zn-zhere)/abs(zn) > deltaz) then
+              zhere = zn
+              goto 20
+          end if
+          height = zn
+       end if
+    end if
+
+    !------------------------------------------------
+    !   Initialize U and V (wind components)
+    !------------------------------------------------
+    d1 = sin(cen_lat*deg2rad)*cos(lat) - &
+         cos(cen_lat*deg2rad)*sin(lat)*cos(lon-cen_lon*deg2rad)
+    d2 = cos(cen_lat*deg2rad)*sin(lon-cen_lon*deg2rad)
+    d  = max(epsilon, sqrt(d1**2.d0 + d2**2.d0))
+    ufac = d1/d
+    vfac = d2/d
+    
+    if (height > ztrop) then
+        u = 0.d0
+        v = 0.d0
+    else
+        v = vfac*(-f*gr/2.d0+sqrt((f*gr/2.d0)**(2.d0) &
+            - exppr*(gr/rp)**exppr*Rd*(T0-gamma*height) &
+            /(exppz*height*Rd*(T0-gamma*height)/(g*zp**exppz) &
+            +(1.d0-p00/dp*exp((gr/rp)**exppr)*exp((height/zp)**exppz)))))
+        u = ufac*(-f*gr/2.d0+sqrt((f*gr/2.d0)**(2.d0) &
+            - exppr*(gr/rp)**exppr*Rd*(T0-gamma*height) &
+            /(exppz*height*Rd*(T0-gamma*height)/(g*zp**exppz) &
+            +(1.d0-p00/dp*exp((gr/rp)**exppr)*exp((height/zp)**exppz)))))
+    end if
+
+    !------------------------------------------------
+    !   Initialize water vapor mixing ratio (q)
+    !------------------------------------------------
+    if (height > ztrop) then
+        q = qtrop
+    else
+        q = q0*exp(-height/zq1)*exp(-(height/zq2)**exppz)
+    end if
+
+    !------------------------------------------------
+    !   Initialize temperature (T)
+    !------------------------------------------------
+    if (height > ztrop) then
+        t = Ttrop
+    else
+        t = (T0-gamma*height)/(1.d0+constTv*q)/(1.d0+exppz*Rd*(T0-gamma*height)*height &
+            /(g*zp**exppz*(1.d0-p00/dp*exp((gr/rp)**exppr)*exp((height/zp)**exppz))))
+    end if
+
+    !-----------------------------------------------------
+    !   Initialize virtual potential temperature (thetav)
+    !-----------------------------------------------------
+    thetav = t * (1.d0+constTv*q) * (p0/p)**(Rd/cp)
+
+    !-----------------------------------------------------
+    !   Initialize surface geopotential (PHIS)
+    !-----------------------------------------------------
+    phis = 0.d0  ! constant
+
+    !-----------------------------------------------------
+    !   Initialize density (rho)
+    !-----------------------------------------------------
+    rho = p/(Rd*t*(1.d0+constTv*q))
+
+  END SUBROUTINE tropical_cyclone_test
+
+!-----------------------------------------------------------------------
+!    First function for fixed point iterations
+!-----------------------------------------------------------------------
+  REAL(8) FUNCTION fpiF(phere, gr, zhere)
+    IMPLICIT NONE
+    REAL(8), INTENT(IN) :: phere, gr, zhere
+
+      fpiF = phere-(p00-dp*exp(-(gr/rp)**exppr)*exp(-(zhere/zp)**exppz)) &
+             *((T0-gamma*zhere)/T0)**(g/(Rd*gamma))
+
+  END FUNCTION fpiF
+
+!-----------------------------------------------------------------------
+!    Second function for fixed point iterations
+!-----------------------------------------------------------------------
+  REAL(8) FUNCTION fpidFdz(gr, zhere) 
+    IMPLICIT NONE
+    REAL(8), INTENT(IN) :: gr, zhere
+
+      fpidFdz =-exppz*zhere*dp*exp(-(gr/rp)**exppr)*exp(-(zhere/zp)**exppz)/(zp*zp)*((T0-gamma*zhere)/T0)**(g/(Rd*gamma)) &
+               +g/(Rd*T0)*(p00-dp*exp(-(gr/rp)**exppr)*exp(-(zhere/zp)**exppz))*((T0-gamma*zhere)/T0)**(g/(Rd*gamma)-1.d0)
+
+  END FUNCTION fpidFdz
+
+END MODULE dcmip2016_cyclone_mod
+
Index: codes/icosagcm/trunk/src/etat0_dcmip2016_cyclone.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0_dcmip2016_cyclone.f90	(revision 413)
+++ codes/icosagcm/trunk/src/etat0_dcmip2016_cyclone.f90	(revision 414)
@@ -3,21 +3,5 @@
   IMPLICIT NONE
   PRIVATE
-  REAL(rstd),PARAMETER :: zt=15000
-  REAL(rstd),PARAMETER :: q0=0.021
-  REAL(rstd),PARAMETER :: qt=1e-11
-  REAL(rstd),PARAMETER :: T0=302.15
-  REAL(rstd),PARAMETER :: Ts=302.15
-  REAL(rstd),PARAMETER :: zq1=3000
-  REAL(rstd),PARAMETER :: zq2=8000
-  REAL(rstd),PARAMETER :: Gamma=0.007
-  REAL(rstd),PARAMETER :: pb=101500
-  REAL(rstd),PARAMETER :: Deltap=1115
-  REAL(rstd),PARAMETER :: rp=282000
-  REAL(rstd),PARAMETER :: zp=7000
-  REAL(rstd),PARAMETER :: epsilon=1e-25
-  REAL(rstd),PARAMETER :: Rd=287
-  
-  REAL(rstd), PARAMETER :: lonc=Pi, latc=Pi/18, &
-       Tv0=T0*(1+0.608*q0), Tvt=Tv0-Gamma*zt
+
 
   PUBLIC getin_etat0, compute_etat0
@@ -28,7 +12,7 @@
     USE mpipara, ONLY : is_mpi_root
 
-    IF(nqtot<1) THEN
+    IF(nqtot<5) THEN
        IF (is_mpi_root)  THEN
-          PRINT *, "nqtot must be at least 1 for test case DCMIP5"
+          PRINT *, "nqtot must be at least 5 for test case DCMIP2016 tropical cyclone"
        END IF
        STOP
@@ -38,6 +22,9 @@
 
   SUBROUTINE compute_etat0(ngrid,lon,lat, phis, ps, Temp, ulon, ulat, q)
-    USE disvert_mod
-    USE omp_para
+  USE disvert_mod
+  USE omp_para
+  USE dcmip2016_cyclone_mod, ONLY : tropical_cyclone_test
+  USE terminator, ONLY: initial_value_Terminator
+  IMPLICIT NONE
     INTEGER, INTENT(IN)    :: ngrid
     REAL(rstd),INTENT(IN)  :: lon(ngrid)
@@ -50,95 +37,33 @@
     REAL(rstd),INTENT(OUT) :: q(ngrid,llm,nqtot)
 
+    INTEGER :: zcoords
+    REAL :: p,z
+    REAL :: rho, thetav
     INTEGER :: l,ij
-    REAL(rstd) :: r, d, d1,d2, zz, Tave, Tp, vt, aa, bb
 
-    phis(:)=0.
-
+   
     DO ij=1,ngrid
-       r=radius*arc(lonc,latc, lon(ij),lat(ij))
-       ps(ij)=pb-DeltaP*exp(-(r/rp)**1.5)
-    END DO
-
-    DO l=ll_begin,ll_end
-       aa=.5*(ap(l)+ap(l+1))
-       bb=.5*(bp(l)+bp(l+1))
-       DO ij=1,ngrid
-          r=radius*arc(lonc,latc, lon(ij),lat(ij))
-          CALL compute_z(aa+bb*ps(ij),ps(ij),r,zz)
-          IF (zz<=zt) THEN
-             q(ij,l,1)=q0*exp(-zz/zq1)*exp(-(zz/zq2)**2)
-             Tave=Tv0-Gamma*zz
-             Tp=(Tv0-Gamma*zz) * ( ( 1 + 2*Rd*(Tv0-Gamma*zz)*zz /                          &
-                  (g*zp**2*(1-pb/Deltap*exp((r/rp)**1.5+(zz/zp)**2))))**(-1) - 1) 
-             vt=-omega*sin(latc)*r+sqrt( (omega*sin(latc)*r)**2 -  1.5*(r/rp)**1.5 * (Tv0-Gamma*zz)*Rd    &
-                  / (1 + 2*Rd*(Tv0-Gamma*zz)*zz/(g*zp**2)        &
-                  - pb/Deltap * exp((r/rp)**1.5 + (zz/zp)**2 )))
-          ELSE
-             q(ij,l,1)=qt
-             Tave=Tvt
-             Tp=0
-             vt=0
-          ENDIF
-          Temp(ij,l)=Tave+Tp 
-
-          d1=sin(latc)*cos(lat(ij))-cos(latc)*sin(lat(ij))*cos(lon(ij)-lonc)
-          d2=cos(latc)*sin(lon(ij)-lonc)
-          d=MAX(1e-25,sqrt(d1**2+d2**2))
-          ulon(ij,l)=vt*d1/d
-          ulat(ij,l)=vt*d2/d
+       z=0.
+       zcoords=1                 
+       CALL tropical_cyclone_test(lon(ij),lat(ij),p,z,zcoords,ulon(ij,1),ulat(ij,1), &
+                                 temp(ij,1),thetav,phis(ij),ps(ij),rho,q(ij,1,1))
+       
+       zcoords=0
+       DO l=ll_begin,ll_end
+         p=0.5*(ap(l)+ap(l+1) + (bp(l)+bp(l+1)) * ps(ij))
+         CALL tropical_cyclone_test(lon(ij),lat(ij),p,z,zcoords,ulon(ij,l),ulat(ij,l), &
+                                   temp(ij,l),thetav,phis(ij),ps(ij),rho,q(ij,l,1))
+                                   
+          IF (physics_thermo==thermo_fake_moist) temp(ij,l)=Temp(ij,l)*(1+0.608*q(ij,l,1))
+          q(ij,l,2)=0.
+          q(ij,l,3)=0.
+          CALL initial_value_Terminator(lat(ij),lon(ij),q(ij,l,4),q(ij,l,5))
        END DO
     ENDDO
+
+
   END SUBROUTINE compute_etat0
   
-  SUBROUTINE compute_z(pmodel,ps,r,z)
-  USE icosa
-  IMPLICIT NONE
-     REAL(rstd),PARAMETER   :: epsilon0=2e-13
-     REAL(rstd),INTENT(IN)  :: pmodel
-     REAL(rstd),INTENT(IN)  :: ps
-     REAL(rstd),INTENT(IN)  :: r
-     REAL(rstd),INTENT(OUT) :: z
-     
-     REAL(rstd) :: pt, pave, pp, dpdz
-     REAL(rstd) :: znp1
-     REAL(rstd) :: epsilon
-     REAL(rstd) :: p
-     INTEGER  :: n     
 
-     pt=pb*(Tvt/Tv0)**(g/(Rd*Gamma))
-
-     IF (pmodel>pt) THEN
-       z=Tv0/Gamma*(1-(pmodel/ps)**(Rd*Gamma/g))
-     ELSE
-       z=zt+Rd*Tvt/g*log(Pt/pmodel)
-     ENDIF
-     
-     IF (z>zt .OR. r>1000000) return
-     
-     epsilon=1
-     n=0
-     
-     DO WHILE(epsilon>epsilon0 .AND. n<20)
-  
-       IF (z<zt) THEN
-         pave=pb*((Tv0-Gamma*z)/Tv0)**(g/(Rd*Gamma))
-         pp= -Deltap * exp(-(r/rp)**1.5-(z/zp)**2) *  ( (Tv0-Gamma*z)/Tv0 )**(g/(Rd*Gamma))
-       ELSE
-         pave=pt*exp(g*(zt-z)/(Rd*Tvt))
-         pp=0
-       ENDIF
-       p=pave+pp
-     
-       dpdz =  2*Deltap*z/zp**2 * exp(-(r/rp)**1.5) * exp(-(z/zp)**2) * ((Tv0-Gamma*z)/Tv0)**(g/(Rd*Gamma))  &
-             - g/(Rd*Tv0) * (preff - Deltap*exp(-(r/rp)**1.5) * exp(-(z/zp)**2)) * ( (Tv0-Gamma*z)/Tv0 )**(g/(Rd*Gamma)-1)
-           
-       znp1 = z - ( pmodel - p ) / (-dpdz)
-       
-       epsilon=ABS( (znp1-z)/znp1 )
-       z=znp1
-       n=n+1
-     ENDDO  
-   
-   END SUBROUTINE compute_z
 
 END MODULE etat0_dcmip2016_cyclone_mod
