[SCM] Gerris Flow Solver branch, upstream, updated. e8f73a07832050124d2b8bf6c6f35b33180e65a8

[Stephane Popinet]

The following commit has been merged in the upstream branch:
commit 8d7b04ce68553e8e5e01ef12bef4dd2f3943a168
Author: Stephane Popinet <popinet at users.sf.net>
Date:   Wed Aug 5 19:56:50 2009 +1000

    Added stability correction for wind wave model
    
    darcs-hash:20090805095650-d4795-f7da99ed29884576c7ac08ece613df5a047737b0.gz

diff --git a/modules/wavewatch.mod b/modules/wavewatch.mod
index 133f82b..75b770c 100644
--- a/modules/wavewatch.mod
+++ b/modules/wavewatch.mod
@@ -39,7 +39,7 @@ static double theta (guint ith, guint ntheta)
 
 typedef struct {
   GfsWave * wave;
-  GfsVariable * ustar, * fpi, * u10, * v10;
+  GfsVariable * ustar, * fpi, * u10, * v10, * as;
   REAL * A;     /* Actions (NK*NTHETA) */
   REAL * CG;    /* Group velocities (NK) */
   REAL * WN;    /* Wavenumbers (NK) */
@@ -70,17 +70,39 @@ static void action_to_energy (FttCell * cell, SourceParams * p)
     }
 }
 
+static void stability_correction (FttCell * cell, SourceParams * p, REAL * U10ABS, REAL * U10DIR)
+{
+  double u10 = p->u10 ? GFS_VALUE (cell, p->u10) : 0.;
+  double v10 = p->v10 ? GFS_VALUE (cell, p->v10) : 0.;
+  *U10ABS = sqrt (u10*u10 + v10*v10);
+  *U10DIR = atan2 (v10, u10);
+  if (p->as) {
+    /* see p.24 of wavewatch manual version 3.12 */
+    double shstab = 1.4, ofstab = -0.01, ccng = -0.1, ccps = 0.1, ffng = -150., ffps = 150.;
+    double zwind = 10.;  /* height at which the wind is defined */
+    double grav = 9.81; /* acceleration of gravity */
+    double stab0 = zwind*grav/273.;
+    double max = MAX (5., *U10ABS);
+    double stab = stab0*GFS_VALUE (cell, p->as)/(max*max);
+    stab = MAX (-1., MIN (1., stab));
+    double tharg1 = MAX (0., ffng*(stab - ofstab));
+    double tharg2 = MAX (0., ffps*(stab - ofstab));
+    double cor1 = ccng*tanh (tharg1);
+    double cor2 = ccps*tanh (tharg2);
+    double cor = sqrt ((1. + cor1 + cor2)/shstab);
+    *U10ABS /= cor;
+  }
+}
+
 static void source (FttCell * cell, SourceParams * p)
 {
   energy_to_action (cell, p);
 
   INTEGER IX, IY, IMOD = 1;
   REAL DEPTH = 1000.; /* fixme: depth is fixed at 1000 m for now */
+  REAL U10ABS, U10DIR;
 
-  double u10 = p->u10 ? GFS_VALUE (cell, p->u10) : 0.;
-  double v10 = p->v10 ? GFS_VALUE (cell, p->v10) : 0.;
-  REAL U10ABS = sqrt (u10*u10 + v10*v10);
-  REAL U10DIR = atan2 (v10, u10);
+  stability_correction (cell, p, &U10ABS, &U10DIR);
 
   REAL USTAR = GFS_VALUE (cell, p->ustar);
   REAL FPI = GFS_VALUE (cell, p->fpi);
@@ -218,6 +240,7 @@ static void wavewatch_source (GfsWave * wave)
   p.fpi = gfs_variable_from_name (domain->variables, "Fpi");
   p.u10 = gfs_variable_from_name (domain->variables, "U10");
   p.v10 = gfs_variable_from_name (domain->variables, "V10");
+  p.as = gfs_variable_from_name (domain->variables, "AS");
 
   guint i;
   for (i = 0; i < wave->nk; i++) {

-- 
Gerris Flow Solver