Merge branch 'fortran-kepert' into ensemble

wind/fwind.f90

@@ -1,15 +1,17 @@
-subroutine fkerpert(R, lam, V, Z, f, rMax, vFm, thetaFm, Vm, Ux, Uy, n)
+subroutine fkerpert(R, lam, f, rMax, vFm, thetaFm, Vm, Ux, Uy, n)
    !$ use omp_lib
 
    integer, intent(in) :: n
    doubleprecision, intent(in) :: f, rMax, vFm, Vm
-   doubleprecision, dimension(n), intent(in) :: R, lam, V, Z
+   doubleprecision, dimension(n), intent(in) :: R, lam
    doubleprecision, dimension(n), intent(inout) :: Ux, Uy
 
    doubleprecision :: Umod, Vt, al, be, gam , K, Cd, u0s, v0s, chi, eta, psi, albe, b
    doubleprecision :: ups, ums, vps, vms, usf, vsf, phi, us, vs
    complex(8) :: A0, j, Am, Ap
    logical :: ind
+   doubleprecision :: aa, bb, cc, delta, edelta, V, Z
+   logical :: icore
 
    b = 1.0
    K = 50.0
@@ -19,6 +21,30 @@ subroutine fkerpert(R, lam, V, Z, f, rMax, vFm, thetaFm, Vm, Ux, Uy, n)
    !$OMP PARALLEL DO shared(Ux, Uy)
    do i = 1, n
 
+      aa = ((d2Vm / 2. - (dVm - vMax / rMax) / rMax) / rMax)
+      bb = (d2Vm - 6 * aa * rMax) / 2.
+      cc = dVm - 3 * aa * rMax ** 2 - 2 * bb * rMax
+
+      delta = (rMax / R(i)) ** beta
+      edelta = exp(-delta)
+      icore = (R(i) <= rMax)
+      if (icore) then
+         Z = R(i) * (R(i) * 4 * aa + 3 * bb) + 2 * cc
+         V = (R(i) * (R(i) * (R(i) * aa + bb) + cc))
+      else
+         V = sqrt((dP * beta / rho) * delta * edelta + (R(i) * f / 2.) ** 2) - R(i) * abs(f) / 2.
+         Z = abs(f) + &
+            (beta**2 * dP * (delta**2) * edelta / &
+             (2 * rho * R(i)) - beta**2 * dP * delta * edelta / &
+             (2 * rho * R(i)) + R(i) * f**2 / 4) / &
+            sqrt(beta * dP * delta * edelta / &
+                    rho + (R(i) * f / 2)**2) + &
+            (sqrt(beta * dP * delta * edelta / &
+                     rho + (R(i) * f / 2)**2)) / R(i)
+      end if
+      V = sign(V, f)
+      Z = sign(Z, f)
+
       if ((vFm > 0) .and. (Vm/vFm < 5.)) then
          Umod = vFm * abs(1.25*(1. - (vFm/Vm)))
       else
@@ -31,18 +57,18 @@ subroutine fkerpert(R, lam, V, Z, f, rMax, vFm, thetaFm, Vm, Ux, Uy, n)
          Vt = Umod
       end if
 
-      al = ((2. * V(i) / R(i)) + f) / (2. * K)
-      be = (f + Z(i)) / (2. * K)
-      gam = (V(i) / (2. * K * R(i)))
+      al = ((2. * V / R(i)) + f) / (2. * K)
+      be = (f + Z) / (2. * K)
+      gam = (V / (2. * K * R(i)))
 
       albe = sqrt(al / be)
 
       ind = abs(gam) > sqrt(al * be)
-      chi = abs((Cd / K) * V(i) / sqrt(sqrt(al * be)))
-      eta = abs((Cd / K) * V(i) / sqrt(sqrt(al * be) + abs(gam)))
-      psi = abs((Cd / K) * V(i) / sqrt(abs(sqrt(al * be) - abs(gam))))
+      chi = abs((Cd / K) * V / sqrt(sqrt(al * be)))
+      eta = abs((Cd / K) * V / sqrt(sqrt(al * be) + abs(gam)))
+      psi = abs((Cd / K) * V / sqrt(abs(sqrt(al * be) - abs(gam))))
 
-      A0 = -(chi * (1.0 + j * (1.0 + chi)) * V(i)) / (2.0 * chi**2 + 3.0 * chi + 2.0)
+      A0 = -(chi * (1.0 + j * (1.0 + chi)) * V) / (2.0 * chi**2 + 3.0 * chi + 2.0)
 
       u0s = realpart(A0) * albe * sign(b, f)
       v0s = imagpart(A0)
@@ -69,7 +95,7 @@ subroutine fkerpert(R, lam, V, Z, f, rMax, vFm, thetaFm, Vm, Ux, Uy, n)
 
       ! Total surface wind in (moving coordinate system)
       us = u0s + ups + ums
-      vs = v0s + vps + vms + V(i)
+      vs = v0s + vps + vms + V
 
       usf = us + Vt * cos(lam(i) - thetaFm)
       vsf = vs - Vt * sin(lam(i) - thetaFm)
@@ -154,4 +180,44 @@ subroutine fhollandvort(Z, R, d2Vm, dVm, rMax, vMax, beta, dP, rho, f, n)
    end do
    !$OMP END PARALLEL DO
 
-end subroutine fhollandvort
+end subroutine fhollandvort
+
+
+subroutine fcomb(V, Z, R, d2Vm, dVm, rMax, vMax, beta, dP, rho, f, n)
+   !$ use omp_lib
+   doubleprecision, intent(in) :: d2Vm, dVm, rMax, beta, dP, rho, f, vMax
+   integer, intent(in) :: n
+   doubleprecision, intent(in), dimension(n) :: R
+   doubleprecision, intent(inout), dimension(n) :: V, Z
+   doubleprecision :: aa, bb, cc, delta, edelta
+   logical :: icore
+
+   aa = (d2Vm / 2 - (dVm - vMax / rMax) / rMax) / rMax
+   bb = (d2Vm - 6 * aa * rMax) / 2.
+   cc = dVm - 3 * aa * rMax ** 2 - 2 * bb * rMax
+
+   !$OMP PARALLEL DO shared(Z)
+   do i = 1, n
+      delta = (rMax / R(i)) ** beta
+      edelta = exp(-delta)
+      icore = (R(i) <= rMax)
+      if (icore) then
+         Z(i) = R(i) * (R(i) * 4 * aa + 3 * bb) + 2 * cc
+         V(i) = (R(i) * (R(i) * (R(i) * aa + bb) + cc))
+      else
+         V(i) = sqrt((dP * beta / rho) * delta * edelta + (R(i) * f / 2.) ** 2) - R(i) * abs(f) / 2.
+         Z(i) = abs(f) + &
+            (beta**2 * dP * (delta**2) * edelta / &
+             (2 * rho * R(i)) - beta**2 * dP * delta * edelta / &
+             (2 * rho * R(i)) + R(i) * f**2 / 4) / &
+            sqrt(beta * dP * delta * edelta / &
+                    rho + (R(i) * f / 2)**2) + &
+            (sqrt(beta * dP * delta * edelta / &
+                     rho + (R(i) * f / 2)**2)) / R(i)
+      end if
+
+      Z(i) = sign(Z(i), f)
+   end do
+   !$OMP END PARALLEL DO
+
+end subroutine fcomb

wind/windmodels.py

@@ -1076,8 +1076,6 @@ def field(self, R, lam, vFm, thetaFm, thetaMax=0.):
 
         """
 
-        V = self.velocity(R)
-        Z = self.vorticity(R)
         K = 50.  # Diffusivity
         Cd = 0.002  # Constant drag coefficient
         Vm = self.profile.vMax
@@ -1087,13 +1085,15 @@ def field(self, R, lam, vFm, thetaFm, thetaMax=0.):
             Ux, Vy = np.empty_like(R), np.empty_like(R)
             n = Ux.size
             fkerpert(
-                R.ravel(), lam.ravel(), V.ravel(), Z.ravel(), self.f, self.rMax, vFm, thetaFm,
+                R.ravel(), lam.ravel(), self.f, self.rMax, vFm, thetaFm,
                 Vm, Ux.ravel(), Vy.ravel(), n
             )
             return Ux, Vy
         except ImportError:
             pass
 
+        V = self.velocity(R)
+        Z = self.vorticity(R)
         if (vFm > 0) and (Vm/vFm < 5.):
             Umod = vFm * np.abs(1.25*(1. - (vFm/Vm)))
         else: