MODULE tropical_cyclone

!=======================================================================
!
!  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
          if (n.gt.40) then

              !PRINT *,'FPI did not converge after 20 interations in q & T!!!'
PRINT *,'FPI did not converge after 40 interations in q & T!!!'
          else if ( abs(zn-zhere)/abs(zn) > deltaz) then
              zhere = zn
              goto 20
          end if
          height = zn
       end if
       z=height
    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 tropical_cyclone 
