Merge pull request #4 from BrodiePearson/ocean/AddStokesTerms_ADC

Adds namelist parameters for Stokes drift/Langmuir forcing

src/core_ocean/adc_mixing/Registry_adc_mixing_fields_opts.xml

@@ -19,6 +19,22 @@
 	      description="number of decimal places to keep, inverse is position of decimal rounding"
 	      possible_values="large positive reals"
 	      />
+  <nml_option name="config_adc_langmuir" type="logical" default_value=".false." units="unitless"
+          description="Switch to turn on wave effects (Langmuir turbulence) in ADC scheme. If false the Stokes drift of waves is set to zero. "
+          possible_values=".true. or .false."
+          />
+  <nml_option name="config_adc_La_t" type="real" default_value="0.3" units="unitless"
+          description="Turbulent Langmuir number La_t = SQRT(friction velocity / surface Stokes drift magnitude). Default value is that of an 'equilibrium wind-sea'."
+          possible_values="positive real numbers, increasing values indicate weaker waves"
+          />
+  <nml_option name="config_adc_stokes_depth" type="real" default_value="4.8" units="m"
+          description="Exponential decay depth of Stokes drift profile. Default value is common in LES studies."
+          possible_values="positive real numbers"
+          />
+  <nml_option name="config_adc_stokes_angle" type="real" default_value="0.0" units="radians"
+          description="Angle between wind and waves (surface stress and surface Stokes drift vector respectively). Positive values denote waves to the left of the wind. NOTE: units are RADIANS."
+          possible_values="positive or negative real numbers"
+          />
   <nml_option name="config_adc_tau_o" type="real" default_value="1800" units="s^{-1}"
 	      description="characterstic eddy turnover timescale"
 	      possible_values="positive real numbers"

src/core_ocean/shared/mpas_ocn_adcReconstruct.F

@@ -178,7 +178,7 @@ subroutine compute_ADC_tends(nCells,nVertLevels, nTracers, dt,activeTracers, uve
       real(kind=RKIND) :: B, Cval, diff, wtav, dzmid, Ksps, Sz, Tz, w4k, w4kp1, w2k, w2kp1
       real(kind=RKIND) :: lareaFraction, wstar, Q, w3av, tempMoment, frictionVelocity
       real(kind=RKIND) :: sfcFrictionVelocitySquared, wtSumUp, wtSumDn, wsSumUp, wsSumDn
-      real(kind=RKIND) :: uStokes_surf, vStokes_surf, stokes_decay_depth, La_t
+      real(kind=RKIND) :: uStokes_surf, vStokes_surf
 
       real(kind=RKIND),dimension(nVertLevels,nCells) :: Swumd
       real(kind=RKIND),dimension(nVertLevels,nCells) :: tauw3, tau_tracer, tauVel, tauvVel
@@ -219,20 +219,25 @@ subroutine compute_ADC_tends(nCells,nVertLevels, nTracers, dt,activeTracers, uve
               frictionVelocity = sqrt(sfcFrictionVelocitySquared + config_adc_bc_wstar * wstar * wstar)
               frictionVelocity = max( config_adc_frictionVelocityMin, frictionVelocity )
 
+              ! Initialize Stokes drift surface values as zero, overwrite these if config_adc_langmuir = .True.
+              uStokes_surf = 0.0_RKIND
+              vStokes_surf = 0.0_RKIND
+              if (config_adc_langmuir) then
+                  uStokes_surf = COS(stokes_angle)*frictionVelocity / (La_t**2.0_RKIND)
+                  vStokes_surf = SIN(stokes_angle)*frictionVelocity / (La_t**2.0_RKIND)
+              endif
+
 !             Calculate exponential Stokes drift profiles from surface values & decay depth
 !             [stokes_decay_depth]. This is temporary; eventually MPAS will pass profiles
-!             For now this is calculated using an assumed Langmuir number of La_t and waves in
-!             x/u-direction only
-!             La_t = (friction velocity / surface stokes drift)^(0.5)
-              La_t = 0.3_RKIND
-              stokes_decay_depth = 4.8_RKIND
-              uStokes_surf = 0.0_RKIND !frictionVelocity / (La_t**2.0_RKIND)
-              vStokes_surf = 0.0_RKIND
+!             For now this is calculated using the turbulent Langmuir number La_t
+!             La_t = (friction velocity / surface stokes drift magnitude)^(0.5) and the
+!             stokes_angle (angle between surface stress and Stokes drift vectors)
+
               uStokes(1,iCell) = uStokes_surf
               vStokes(1,iCell) = vStokes_surf
               do k=1,nVertLevels
-                 uStokes(k,iCell) = uStokes_surf*EXP(ze(k,iCell)/stokes_decay_depth)
-                 vStokes(k,iCell) = vStokes_surf*EXP(ze(k,iCell)/stokes_decay_depth)
+                 uStokes(k,iCell) = uStokes_surf*EXP(ze(k,iCell)/stokes_depth)
+                 vStokes(k,iCell) = vStokes_surf*EXP(ze(k,iCell)/stokes_depth)
               enddo
 
               do k=1,3

src/core_ocean/shared/mpas_ocn_turbulence.F

@@ -48,7 +48,8 @@ module ocn_turbulence
                        alpha_st_tracer1, alpha_st_tracer2, alpha_st_0, alpha_st_1, &
                        alpha_st_2, alpha_0, alpha_1, alpha_2, B1, Kt, grav, c_mom, &
                        c_therm, c_mom_w3, c_epsilon, c_slow, c_slow_tracer, slow_w_factor, &
-                       kappa_FL, kappa_w3, kappa_VAR, Cww_D, Cww_E, adcRound
+                       kappa_FL, kappa_w3, kappa_VAR, Cww_D, Cww_E, adcRound, &
+                       La_t, stokes_depth, stokes_angle
 
    integer,public :: iterCount, nCells, nVertLevels
 
@@ -182,6 +183,9 @@ subroutine ocn_turbulenceCreate(domain)
       kappa_FL = config_adc_kappaFL
       kappa_VAR = config_adc_kappaVAR
       kappa_w3 = config_adc_kappaW3
+      La_t = config_adc_La_t
+      stokes_depth = config_adc_stokes_depth
+      stokes_angle = config_adc_stokes_angle
       adcRound = 10.0_RKIND**config_adc_decimals_to_keep
       iterCount = 1
 
@@ -1243,6 +1247,9 @@ subroutine ocn_turbulenceDestroy(ierr)
       kappa_FL = 0.0_RKIND
       kappa_VAR = 0.0_RKIND
       kappa_w3 = 0.0_RKIND
+      La_t = 0.0_RKIND
+      stokes_depth = 0.0_RKIND
+      stokes_angle = 0.0_RKIND
 
       !$acc exit data delete(KspsD,    &
       !$acc                   nCells,   &