Import README and NCL script to Git repository

README

@@ -0,0 +1,15 @@
+NCL Scripts to draw the Gill patterns
+=====================================
+
+Reference
+--------
+
+* Gill A. E., 1980: Some simple solutions for heat-induced tropical circulation. *Quart. J. Roy. Meteor. Soc.*, **106 (449)**, 447--462.
+
+Files
+-----
+* gill.ncl: parameters and functions representing solutions
+* plot.ncl: plot (u, v) in vectors w in shades and pressure in contours
+* plot_panel.ncl: w and p superimposed by (u, v) in vector in seperate panels
+* plot_hadley.ncl: plot Hadley circulation
+* plot_walker.ncl: plot Walker circulation

gill.ncl

@@ -0,0 +1,154 @@
+pi = acos(-1.0)
+pir = 1.0/pi
+L = 2.0
+k = pi / (2.0 * L)
+eps = 0.1
+nx = 251
+ny =  81
+nz =  30
+xmin = -10.0
+ymin =  -4.0
+ymax =   4.0
+xmax =  15.0
+zmin =   0.0
+zmax =    pi
+I = 4.0 * L * pir
+
+function lon()
+begin
+  return fspan(xmin, xmax, nx)
+end
+
+function lat()
+begin
+  return fspan(ymin, ymax, ny)
+end
+
+function lev()
+begin
+  return fspan(zmin, zmax, nz)
+end
+
+function q0(x)
+local ekr, kx
+begin
+  ekr = 1.0 / (eps * eps + k * k)
+  kx = k * x
+  if (x .le. L) then
+    if (x .le. -L) then
+      return 0.0
+    else
+      return -ekr * (eps * cos(kx) + k * (sin(kx) + exp(-eps * (x + L))))
+    end if
+  else
+    return -ekr * k * (1.0 + exp(-2.0 * eps * L)) * exp(eps * (L - x))
+  end if
+end
+
+function qn(n,x)
+local ekr, kx
+begin
+  nn = 2 * n + 1
+  ekr = 1.0 / (nn * nn * eps * eps + k * k)
+  kx = k * x
+  if (x .le. L) then
+    if (x .le. -L) then
+      return -ekr * k * (1.0 + exp(-2 * nn * eps * L)) * exp(nn * eps * (x + L))
+    else
+      return ekr * (-nn * eps * cos(kx) + k * (sin(kx) - exp(nn * eps * (x - L))))
+    end if
+  else
+    return 0.0
+  end if
+end
+
+function F(x)
+begin
+  if (fabs(x) .lt. L) then
+    return cos(k * x)
+  else
+    return 0.0
+  end if
+end
+
+procedure symmetric(p, u, v, w, x, y)
+local i, j, q0x, q2x
+begin
+  do j = 0, ny - 1
+    do i = 0, nx - 1
+      q0x = q0(x(i))
+      q2x = qn(1, x(i))
+      Fx = F(x(i))
+      ey2 = exp(-0.25 * y(j)^2)
+      p(j, i) = 0.5 * (q0x + q2x * (1.0 + y(j) * y(j))) * ey2
+      u(j, i) = 0.5 * (q0x + q2x * (y(j) * y(j) - 3.0)) * ey2
+      v(j, i) = (Fx + 4.0 * eps * q2x) * y(j) * ey2
+      w(j, i) = (0.5 * eps * q0x + Fx + 0.5 * eps * q2x * (1.0 + y(j) * y(j))) * ey2
+    end do
+  end do
+end
+
+procedure antisymmetric(p, u, v, w, x, y)
+local i, j, q3x, y3, y2, ey2
+begin
+  do j = 0, ny - 1
+    y2 = y(j) * y(j)
+    y3 = y2 * y(j)
+    do i = 0, nx - 1
+      q3x = qn(2,x(i))
+      Fx = F(x(i))
+      ey2 = exp(-0.25 * y(j)^2)
+      p(j, i) = 0.5 * q3x * y3 * ey2
+      u(j, i) = 0.5 * q3x *(y3 - 6 * y(j)) * ey2
+      v(j, i) = (6.0 * eps * q3x * (y2 - 1.0) + Fx * y2) * ey2
+      w(j, i) = (0.5 * eps * q3x * y3 + Fx * y(j)) * ey2
+    end do
+  end do
+end
+
+;procedure hadley_symmetric(p, u, v, w, y, z)
+; kai: clockwise + (meteorology, oceanography as opposed to fluid dynamics)
+procedure hadley_symmetric(p, u, kai, y, z)
+local j, k, yy, eyy
+begin
+  do j = 0, ny-1
+    yy = y(j) * y(j)
+    eyy = I * exp(-0.25 * yy)
+    p(j) = -(4.0 + yy) / (6.0 * eps) * eyy
+    do k = 0, nz-1
+      u(k, j) = -yy / (6.0 * eps) * eyy * cos(z(k))
+;      v(k, j) = -y(j) / 3.0 * eyy * cos(z(k))
+;      w(k, j) = -(2.0 - yy) / 6.0 * eyy * sin(z(k))
+      kai(k, j) = y(j) / 3.0 * eyy * sin(z(k))
+    end do
+  end do
+end
+
+procedure hadley_antisymmetric(p, u, kai, y, z)
+local j, k, yy, eyy
+begin
+  do j = 0, ny-1
+    yy = y(j) * y(j)
+    yyy = yy * y(j)
+    eyy = I * exp(-0.25 * yy)
+    p(j) = -yyy / (10.0 * eps) * eyy
+    do k = 0, nz-1
+      u(k, j) = (6.0 * y(j) - yyy) / (10.0 * eps) * eyy * cos(z(k))
+      kai(k, j) = (yy - 6.0) / 5.0 * eyy * sin(z(k))
+    end do
+  end do
+end
+
+procedure walker(p, kai, x, z)
+local i, k, q0x, q2x, sqrtpi
+begin
+  sqrtpi = sqrt(pi)
+  do i = 0, nx-1
+    q0x = q0(x(i))
+    q2x = qn(1, x(i))
+    p(i) = (q0x + 3.0 * q2x) * sqrtpi
+    do k = 0, nz-1
+      kai(k, i) = (q2x - q0x) * sqrtpi * sin(z(k))
+    end do
+  end do
+end

plot.ncl

@@ -0,0 +1,88 @@
+lsymmetric = True
+;lsymmetric = False
+lboth = False
+;lboth = True
+load "gill.ncl"
+
+begin
+  x = lon()
+  y = lat()
+  p = new((/ny, nx/), "float")
+  p!0 = "y"
+  p&y = y
+  p!1 = "x"
+  p&x = x
+  u = p
+  v = p
+  w = p
+  if (lsymmetric.or.lboth) then
+    symmetric(p, u, v, w, x, y)
+  else
+    antisymmetric(p, u, v, w, x, y)
+  end if
+  if (lboth) then
+    pa = p
+    ua = p
+    va = p
+    wa = p
+    antisymmetric(pa, ua, va, wa, x, y)
+    p = p + pa
+    u = u + ua
+    v = v + va
+    w = w + wa
+  end if
+
+  wks = gsn_open_wks("x11", "gill")
+;  gsn_define_colormap(wks, "BlAqGrWh2YeOrReVi22")
+  gsn_define_colormap(wks, "gs16")
+
+  resp = True
+  resp@vpWidthF = 0.8
+  resp@vpHeightF = 0.4
+  resp@vpXF = 0.1
+  resp@vpYF = 0.7
+  resp@gsnDraw = False
+  resp@gsnFrame = False
+  resv = resp
+  resp@cnLevelSelectionMode = "ExplicitLevels"
+  resp@cnInfoLabelOn = False
+  resw = resp
+
+;  resw@cnLevels = (/      -0.3, -0.1, 0.1, 0.3, 0.6, 0.9/)
+;  resw@cnFillColors = (/ 6,    4,   -1,  8,  9,  10,  11/)
+  resw@cnLevels = (/      -0.3, -0.1, 0.1, 0.3, 0.6, 0.9/)
+  resw@cnFillColors = (/ 4,    4,   -1,  9,  9,  9,  9/)
+  resw@cnFillOn = True
+;  resw@cnLinesOn = False
+  resw@lbLabelFontHeightF = 0.02
+  plotw = gsn_csm_contour(wks, w, resw)
+
+  resp@cnLevels = (/-1.5, -1.2, -0.9, -0.6, -0.3, 0.3, 0.6, 0.9, 1.2, 1.5/)
+  resp@cnFillOn = False
+  resp@cnLineThicknessF = 3
+  resp@cnLineLabelBackgroundColor = -1
+  plotp = gsn_csm_contour(wks, p, resp)
+
+  resv@vcRefAnnoOrthogonalPosF = -1.45
+  resv@vcRefMagnitudeF = 1.0
+  resv@vcRefLengthF = 0.02
+; Causes Segmentation Fault
+;  resv@vcGlyphStyle = "CurlyVector"
+  resv@vcGlyphStyle = "FillArrow"
+  resv@vcMinDistanceF = 0.02
+  resv@vcRefAnnoString2On = False
+;  resv@vcLineArrowThicknessF = 2.0
+
+;  res@cnLevelSelectionMode = "AutomaticLevels"
+;  res@gsnContourNegLineDashPattern = 2
+;  plot = gsn_csm_contour(wks, u, res)
+;  plot = gsn_csm_contour(wks, v, res)
+  plot = gsn_csm_vector(wks, u, v, resv)
+  overlay(plotw, plotp)
+  overlay(plotw, plot)
+
+  resp = True
+
+  draw(plotw)
+  frame(wks)
+end 

plot_hadley.ncl

@@ -0,0 +1,64 @@
+load "gill.ncl"
+;lsymmetric = True
+lsymmetric = False
+lboth = False
+;lboth = True
+begin
+  y = lat()
+  z = lev()
+  p = new((/ny/), "float")
+  p!0 = "y"
+  p&y = y
+  u = new((/nz, ny/), "float")
+  u!0 = "z"
+  u&z = z
+  u!1 = "y"
+  u&y = y
+  kai = u
+
+  wks = gsn_open_wks("x11", "hadley")
+
+  resu = True
+  resu@vpWidthF = 0.8
+  resu@vpHeightF = 0.4
+  resu@vpXF = 0.1
+  resu@vpYF = 0.7
+  resu@gsnDraw = False
+  resu@gsnFrame = False
+  resp = resu
+
+;  hadley_symmetric(p, u, v, w, y, z)
+  if (lsymmetric.or.lboth) then
+    hadley_symmetric(p, u, kai, y, z)
+  else
+    hadley_antisymmetric(p, u, kai, y, z)
+  end if
+  if (lboth) then
+    pa = p
+    ua = u
+    kaia = u
+    hadley_antisymmetric(pa, ua, kaia, y, z)
+    p = p + pa
+    u = u + ua
+    kai = kai + kaia
+  end if
+
+  resu@tmYLMode = "Explicit"
+  resu@tmYLValues = (/z(0), 0.5 * pi, z(nz-1)/)
+  resu@tmYLLabels = (/"0.0", "D/2", "D"/)
+  resu@cnLineThicknessF = 3.0
+  resu@gsnContourNegLineDashPattern = 2
+  resu@gsnContourZeroLineThicknessF = 0.0
+  resk = resu
+
+  plot = new(3, "graphic")
+  plot(0) = gsn_csm_contour(wks, u, resu)
+  plot(1) = gsn_csm_contour(wks, kai, resk)
+  plot(2) = gsn_csm_xy(wks, y, p, resp)
+
+  pres = True
+  pres@gsnPanelFigureStrings = (/"(a)", "(b)", "(c)"/)
+  pres@gsnPanelFigureStringsPerimOn = False
+
+  gsn_panel(wks, plot, (/3, 1/), pres)
+end

plot_panel.ncl

@@ -0,0 +1,66 @@
+lsymmetric = True
+;lsymmetric = False
+load "gill.ncl"
+
+begin
+  x = lon()
+  y = lat()
+  p = new((/ny, nx/), "float")
+  p!0 = "y"
+  p&y = y
+  p!1 = "x"
+  p&x = x
+  u = p
+  v = p
+  w = p
+  if (lsymmetric) then
+    symmetric(p, u, v, w, x, y)
+  else
+    antisymmetric(p, u, v, w, x, y)
+  end if
+
+  wks = gsn_open_wks("x11", "gill_wp")
+
+  resp = True
+  resp@vpWidthF = 0.8
+  resp@vpHeightF = 0.4
+  resp@vpXF = 0.1
+  resp@vpYF = 0.7
+  resp@gsnDraw = False
+  resp@gsnFrame = False
+  resv = resp
+  resp@cnLineThicknessF = 3.0
+  resp@cnLevelSelectionMode = "ExplicitLevels"
+  resp@cnInfoLabelOn = False
+  resw = resp
+
+;  resw@gsnContourZeroLineThicknessF = 2
+;  resw@cnLevels = (/-0.6, 0.3, -0.1, 0.0, 0.1, 0.3, 0.6/)
+  resw@gsnContourNegLineDashPattern = 2
+  resw@cnLevels = (/-0.6, 0.3, -0.1, 0.1, 0.3, 0.6/)
+  plotw = gsn_csm_contour(wks, w, resw)
+
+  resp@gsnContourNegLineDashPattern = 0
+  resp@cnLevels = (/-1.5, -1.2, -0.9, -0.6, -0.3, 0.3, 0.6, 0.9, 1.2, 1.5/)
+  resp@cnFillOn = False
+  plotp = gsn_csm_contour(wks, p, resp)
+
+  resv@vcRefAnnoOrthogonalPosF = -1.15
+  resv@vcRefMagnitudeF = 1.0
+  resv@vcRefLengthF = 0.02
+  resv@vcMinDistanceF = 0.02
+  resv@vcRefAnnoString2On = False
+;  resv@vcMinMagnitudeF = 0.3
+
+  plot = new(2, "graphic")
+
+  plot(0) = gsn_csm_vector(wks, u, v, resv)
+  overlay(plot(0), plotw)
+
+  plot(1) = gsn_csm_vector(wks, u, v, resv)
+  overlay(plot(1), plotp)
+
+  pres = True
+
+  gsn_panel(wks, plot, (/2, 1/), pres)
+end