Index: codes/icosagcm/trunk/src/dcmip2016_baroclinic_wave.f90
===================================================================
--- codes/icosagcm/trunk/src/dcmip2016_baroclinic_wave.f90	(revision 410)
+++ codes/icosagcm/trunk/src/dcmip2016_baroclinic_wave.f90	(revision 410)
@@ -0,0 +1,443 @@
+MODULE dcmip2016_baroclinic_wave_mod
+
+!=======================================================================
+!
+!  Date:  July 29, 2015
+!
+!  Functions for setting up idealized initial conditions for the
+!  Ullrich, Melvin, Staniforth and Jablonowski baroclinic instability.
+!
+!  SUBROUTINE baroclinic_wave_sample(
+!    deep,moist,pertt,X,lon,lat,p,z,zcoords,u,v,w,t,phis,ps,rho,q)
+!
+!  Options:
+!     deep    deep atmosphere (1 = yes or 0 = no)
+!    moist    include moisture (1 = yes or 0 = no)
+!    pertt    type of perturbation (0 = exponential, 1 = stream function)
+!        X    Earth scaling factor
+!
+!  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 and zcoords = 1 (Pa)
+!        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    water vapor mixing ratio (kg/kg)
+!
+!
+!  Author: Paul Ullrich
+!          University of California, Davis
+!          Email: paullrich@ucdavis.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 ::               &
+       T0E        = 310.d0     ,      & ! temperature at equatorial surface (K)
+       T0P        = 240.d0     ,      & ! temperature at polar surface (K)
+       B          = 2.d0       ,      & ! jet half-width parameter
+       K          = 3.d0       ,      & ! jet width parameter
+       lapse      = 0.005d0             ! lapse rate parameter
+
+  REAL(8), PARAMETER ::               &
+       pertu0     = 0.5d0      ,      & ! SF Perturbation wind velocity (m/s)
+       pertr      = 1.d0/6.d0  ,      & ! SF Perturbation radius (Earth radii)
+       pertup     = 1.0d0      ,      & ! Exp. perturbation wind velocity (m/s)
+       pertexpr   = 0.1d0      ,      & ! Exp. perturbation radius (Earth radii)
+       pertlon    = pi/9.d0    ,      & ! Perturbation longitude
+       pertlat    = 2.d0*pi/9.d0,     & ! Perturbation latitude
+       pertz      = 15000.d0   ,      & ! Perturbation height cap
+       dxepsilon  = 1.d-5               ! Small value for numerical derivatives
+ 
+  REAL(8), PARAMETER ::               &
+       moistqlat  = 2.d0*pi/9.d0,     & ! Humidity latitudinal width
+       moistqp    = 34000.d0,         & ! Humidity vertical pressure width
+       moisttr    = 0.1d0,            & ! Vertical cut-off pressure for humidity
+       moistqs    = 1.d-12,           & ! Humidity above cut-off
+       moistq0    = 0.018d0,          & ! Maximum specific humidity
+       moistqr    = 0.9d0,            & ! Maximum saturation ratio
+       moisteps   = 0.622d0,          & ! Ratio of gas constants
+       moistT0    = 273.16d0,         & ! Reference temperature (K)
+       moistE0Ast = 610.78d0            ! Saturation vapor pressure at T0 (Pa) 
+
+CONTAINS
+
+!=======================================================================
+!    Generate the baroclinic instability initial conditions
+!=======================================================================
+  SUBROUTINE baroclinic_wave_test(deep,moist,pertt,X,lon,lat,p,z,zcoords,u,v,t,thetav,phis,ps,rho,q) &
+    BIND(c, name = "baroclinic_wave_test")
+ 
+    IMPLICIT NONE
+
+!-----------------------------------------------------------------------
+!     input/output params parameters at given location
+!-----------------------------------------------------------------------
+    INTEGER, INTENT(IN)  :: &
+                deep,       & ! Deep (1) or Shallow (0) test case
+                moist,      & ! Moist (1) or Dry (0) test case
+                pertt         ! Perturbation type
+
+    REAL(8), INTENT(IN)  :: &
+                lon,        & ! Longitude (radians)
+                lat,        & ! Latitude (radians)
+                X             ! Earth scaling parameter
+
+    REAL(8), INTENT(INOUT) :: &
+                p,            & ! Pressure (Pa)
+                z               ! Altitude (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             ! water vapor mixing ratio (kg/kg)
+
+    !------------------------------------------------
+    !   Local variables
+    !------------------------------------------------
+    REAL(8) :: aref, omegaref
+    REAL(8) :: T0, constH, constC, scaledZ, inttau2, rratio
+    REAL(8) :: inttermU, bigU, rcoslat, omegarcoslat
+    REAL(8) :: eta, qratio, qnum, qden
+
+    !------------------------------------------------
+    !   Pressure and temperature
+    !------------------------------------------------
+    if (zcoords .eq. 1) then
+      CALL evaluate_pressure_temperature(deep, X, lon, lat, z, p, t)
+    else
+      CALL evaluate_z_temperature(deep, X, lon, lat, p, z, t)
+    end if
+
+    !------------------------------------------------
+    !   Compute test case constants
+    !------------------------------------------------
+    aref = a / X
+    omegaref = omega * X
+
+    T0 = 0.5d0 * (T0E + T0P)
+
+    constH = Rd * T0 / g
+
+    constC = 0.5d0 * (K + 2.d0) * (T0E - T0P) / (T0E * T0P)
+
+    scaledZ = z / (B * constH)
+
+    inttau2 = constC * z * exp(- scaledZ**2)
+
+    ! radius ratio
+    if (deep .eq. 0) then
+      rratio = 1.d0
+    else
+      rratio = (z + aref) / aref;
+    end if
+
+    !-----------------------------------------------------
+    !   Initialize surface pressure
+    !-----------------------------------------------------
+    ps = p0
+
+    !-----------------------------------------------------
+    !   Initialize velocity field
+    !-----------------------------------------------------
+    inttermU = (rratio * cos(lat))**(K - 1.d0) - (rratio * cos(lat))**(K + 1.d0)
+    bigU = g / aref * K * inttau2 * inttermU * t
+    if (deep .eq. 0) then
+      rcoslat = aref * cos(lat)
+    else
+      rcoslat = (z + aref) * cos(lat)
+    end if
+
+    omegarcoslat = omegaref * rcoslat
+    
+    u = - omegarcoslat + sqrt(omegarcoslat**2 + rcoslat * bigU)
+    v = 0.d0
+
+    !-----------------------------------------------------
+    !   Add perturbation to the velocity field
+    !-----------------------------------------------------
+
+    ! Exponential type
+    if (pertt .eq. 0) then
+      u = u + evaluate_exponential(lon, lat, z)
+
+    ! Stream function type
+    elseif (pertt .eq. 1) then
+      u = u - 1.d0 / (2.d0 * dxepsilon) *                       &
+          ( evaluate_streamfunction(lon, lat + dxepsilon, z)    &
+          - evaluate_streamfunction(lon, lat - dxepsilon, z))
+
+      v = v + 1.d0 / (2.d0 * dxepsilon * cos(lat)) *            &
+          ( evaluate_streamfunction(lon + dxepsilon, lat, z)    &
+          - evaluate_streamfunction(lon - dxepsilon, lat, z))
+    end if
+
+    !-----------------------------------------------------
+    !   Initialize surface geopotential
+    !-----------------------------------------------------
+    phis = 0.d0
+
+    !-----------------------------------------------------
+    !   Initialize density
+    !-----------------------------------------------------
+    rho = p / (Rd * t)
+
+    !-----------------------------------------------------
+    !   Initialize specific humidity
+    !-----------------------------------------------------
+    if (moist .eq. 1) then
+      eta = p/p0
+
+      if (eta .gt. moisttr) then
+        q = moistq0 * exp(- (lat/moistqlat)**4)          &
+                    * exp(- ((eta-1.d0)*p0/moistqp)**2)
+      else
+        q = moistqs
+      end if
+
+      ! Convert virtual temperature to temperature
+      t = t / (1.d0 + Mvap * q)
+
+    else
+      q = 0.d0
+    end if
+
+    !-----------------------------------------------------
+    !   Initialize virtual potential temperature
+    !-----------------------------------------------------
+    thetav = t * (1.d0 + 0.61d0 * q) * (p0 / p)**(Rd / cp)
+
+  END SUBROUTINE baroclinic_wave_test
+
+!-----------------------------------------------------------------------
+!    Calculate pointwise pressure and temperature
+!-----------------------------------------------------------------------
+  SUBROUTINE evaluate_pressure_temperature(deep, X, lon, lat, z, p, t)
+
+    INTEGER, INTENT(IN)  :: deep ! Deep (1) or Shallow (0) test case
+
+    REAL(8), INTENT(IN)  :: &
+                X,          & ! Earth scaling ratio
+                lon,        & ! Longitude (radians)
+                lat,        & ! Latitude (radians)
+                z             ! Altitude (m)
+
+    REAL(8), INTENT(OUT) :: &
+                p,          & ! Pressure (Pa)
+                t             ! Temperature (K)
+
+    REAL(8) :: aref, omegaref
+    REAL(8) :: T0, constA, constB, constC, constH, scaledZ
+    REAL(8) :: tau1, tau2, inttau1, inttau2
+    REAL(8) :: rratio, inttermT
+
+    !--------------------------------------------
+    ! Constants
+    !--------------------------------------------
+    aref = a / X
+    omegaref = omega * X
+
+    T0 = 0.5d0 * (T0E + T0P)
+    constA = 1.d0 / lapse
+    constB = (T0 - T0P) / (T0 * T0P)
+    constC = 0.5d0 * (K + 2.d0) * (T0E - T0P) / (T0E * T0P)
+    constH = Rd * T0 / g
+
+    scaledZ = z / (B * constH)
+
+    !--------------------------------------------
+    !    tau values
+    !--------------------------------------------
+    tau1 = constA * lapse / T0 * exp(lapse * z / T0) &
+         + constB * (1.d0 - 2.d0 * scaledZ**2) * exp(- scaledZ**2)
+    tau2 = constC * (1.d0 - 2.d0 * scaledZ**2) * exp(- scaledZ**2)
+
+    inttau1 = constA * (exp(lapse * z / T0) - 1.d0) &
+            + constB * z * exp(- scaledZ**2)
+    inttau2 = constC * z * exp(- scaledZ**2)
+
+    !--------------------------------------------
+    !    radius ratio
+    !--------------------------------------------
+    if (deep .eq. 0) then
+      rratio = 1.d0
+    else
+      rratio = (z + aref) / aref;
+    end if
+
+    !--------------------------------------------
+    !    interior term on temperature expression
+    !--------------------------------------------
+    inttermT = (rratio * cos(lat))**K &
+             - K / (K + 2.d0) * (rratio * cos(lat))**(K + 2.d0)
+
+    !--------------------------------------------
+    !    temperature
+    !--------------------------------------------
+    t = 1.d0 / (rratio**2 * (tau1 - tau2 * inttermT))
+
+    !--------------------------------------------
+    !    hydrostatic pressure
+    !--------------------------------------------
+    p = p0 * exp(- g / Rd * (inttau1 - inttau2 * inttermT))
+
+  END SUBROUTINE evaluate_pressure_temperature
+
+!-----------------------------------------------------------------------
+!    Calculate pointwise z and temperature given pressure
+!-----------------------------------------------------------------------
+  SUBROUTINE evaluate_z_temperature(deep, X, lon, lat, p, z, t)
+    
+    INTEGER, INTENT(IN)  :: deep ! Deep (1) or Shallow (0) test case
+
+    REAL(8), INTENT(IN)  :: &
+                X,          & ! Earth scaling ratio
+                lon,        & ! Longitude (radians)
+                lat,        & ! Latitude (radians)
+                p             ! Pressure (Pa)
+
+    REAL(8), INTENT(OUT) :: &
+                z,          & ! Altitude (m)
+                t             ! Temperature (K)
+
+    INTEGER :: ix
+
+    REAL(8) :: z0, z1, z2
+    REAL(8) :: p0, p1, p2
+
+    z0 = 0.d0
+    z1 = 10000.d0
+
+    CALL evaluate_pressure_temperature(deep, X, lon, lat, z0, p0, t)
+    CALL evaluate_pressure_temperature(deep, X, lon, lat, z1, p1, t)
+
+    DO ix = 1, 100
+      z2 = z1 - (p1 - p) * (z1 - z0) / (p1 - p0)
+
+      CALL evaluate_pressure_temperature(deep, X, lon, lat, z2, p2, t)
+
+      IF (ABS((p2 - p)/p) .lt. 1.0d-13) THEN
+        EXIT
+      END IF
+
+      z0 = z1
+      p0 = p1
+
+      z1 = z2
+      p1 = p2
+    END DO
+
+    z = z2
+
+    CALL evaluate_pressure_temperature(deep, X, lon, lat, z, p0, t)
+
+  END SUBROUTINE evaluate_z_temperature
+
+!-----------------------------------------------------------------------
+!    Exponential perturbation function
+!-----------------------------------------------------------------------
+  REAL(8) FUNCTION evaluate_exponential(lon, lat, z)
+
+    REAL(8), INTENT(IN)  :: &
+                lon,        & ! Longitude (radians)
+                lat,        & ! Latitude (radians)
+                z             ! Altitude (meters)
+
+    REAL(8) :: greatcircler, perttaper
+
+    ! Great circle distance
+    greatcircler = 1.d0 / pertexpr &
+      * acos(sin(pertlat) * sin(lat) + cos(pertlat) * cos(lat) * cos(lon - pertlon))
+
+    ! Vertical tapering of stream function
+    if (z < pertz) then
+      perttaper = 1.d0 - 3.d0 * z**2 / pertz**2 + 2.d0 * z**3 / pertz**3
+    else
+      perttaper = 0.d0
+    end if
+
+    ! Zonal velocity perturbation
+    if (greatcircler < 1.d0) then
+      evaluate_exponential = pertup * perttaper * exp(- greatcircler**2)
+    else
+      evaluate_exponential = 0.d0
+    end if
+
+  END FUNCTION evaluate_exponential
+
+!-----------------------------------------------------------------------
+!    Stream function perturbation function
+!-----------------------------------------------------------------------
+  REAL(8) FUNCTION evaluate_streamfunction(lon, lat, z)
+
+    REAL(8), INTENT(IN)  :: &
+                lon,        & ! Longitude (radians)
+                lat,        & ! Latitude (radians)
+                z             ! Altitude (meters)
+
+    REAL(8) :: greatcircler, perttaper, cospert
+
+    ! Great circle distance
+    greatcircler = 1.d0 / pertr &
+      * acos(sin(pertlat) * sin(lat) + cos(pertlat) * cos(lat) * cos(lon - pertlon))
+
+    ! Vertical tapering of stream function
+    if (z < pertz) then
+      perttaper = 1.d0 - 3.d0 * z**2 / pertz**2 + 2.d0 * z**3 / pertz**3
+    else
+      perttaper = 0.d0
+    end if
+
+    ! Horizontal tapering of stream function
+    if (greatcircler .lt. 1.d0) then
+      cospert = cos(0.5d0 * pi * greatcircler)
+    else
+      cospert = 0.d0
+    end if
+
+    evaluate_streamfunction = &
+        (- pertu0 * pertr * perttaper * cospert**4)
+
+  END FUNCTION evaluate_streamfunction
+
+END MODULE dcmip2016_baroclinic_wave_mod
+
