Langmuir

src/gotm/diagnostics.F90

@@ -92,9 +92,10 @@ subroutine do_diagnostics(nlev)
    use density,      only: rho0,cp
    use meanflow,     only: gravity,drag
    use meanflow,     only: h,u,v,s,t,NN,SS,buoy,rad
+   use stokes_drift, only: dusdz, dvsdz
    use turbulence,   only: turb_method
    use turbulence,   only: kappa
-   use turbulence,   only: num,nuh
+   use turbulence,   only: num,nuh, nucl
    use turbulence,   only: tke,P,B
    use turbulence,   only: Rig   
    use observations, only: tprof_input,b_obs_sbf
@@ -151,8 +152,10 @@ subroutine do_diagnostics(nlev)
    taux(nlev) = -tx
    tauy(nlev) = -ty
    do i=nlev-1,1,-1
-      taux(i)=-num(i)*(u(i+1)-u(i))/(0.5*(h(i+1)+h(i)))
-      tauy(i)=-num(i)*(v(i+1)-v(i))/(0.5*(h(i+1)+h(i)))
+      taux(i)=-num(i)*(u(i+1)-u(i))/(0.5*(h(i+1)+h(i)))                &
+              -nucl(i)*dusdz%data(i)
+      tauy(i)=-num(i)*(v(i+1)-v(i))/(0.5*(h(i+1)+h(i)))                &
+              -nucl(i)*dvsdz%data(i)
    end do
    taux(0)=-drag(1)*u(1)*sqrt(u(1)**2+v(1)**2)
    tauy(0)=-drag(1)*v(1)*sqrt(u(1)**2+v(1)**2)

src/gotm/gotm.F90

@@ -56,7 +56,7 @@ module gotm
 
    use turbulence,  only: turb_method
    use turbulence,  only: init_turbulence,post_init_turbulence,do_turbulence
-   use turbulence,  only: num,nuh,nus
+   use turbulence,  only: num,nuh,nus, nucl
    use turbulence,  only: const_num,const_nuh
    use turbulence,  only: gamu,gamv,gamh,gams
    use turbulence,  only: Rig
@@ -792,8 +792,8 @@ subroutine integrate_gotm()
       call coriolis(nlev,dt)
 
 !     update velocity
-      call uequation(nlev,dt,cnpar,tx,num,gamu,ext_press_mode)
-      call vequation(nlev,dt,cnpar,ty,num,gamv,ext_press_mode)
+      call uequation(nlev,dt,cnpar,tx,num, nucl, gamu,ext_press_mode)
+      call vequation(nlev,dt,cnpar,ty,num, nucl, gamv,ext_press_mode)
       call external_pressure(ext_press_mode,nlev)
       call internal_pressure(nlev)
       call friction(nlev,kappa,avmolu,tx,ty,plume_type)
@@ -889,10 +889,10 @@ subroutine integrate_gotm()
 !        update one-point models
 # ifdef SEAGRASS
          call do_turbulence(nlev,dt,depth,u_taus,u_taub,z0s,z0b,h,      &
-                            NN,SS,xP, SSCSTK=SSCSTK)
+                            NN,SS,xP, SSCSTK=SSCSTK, SSSTK=SSSTK)
 # else
          call do_turbulence(nlev,dt,depth,u_taus,u_taub,z0s,z0b,h,      &
-                            NN,SS, SSCSTK=SSCSTK)
+                            NN,SS, SSCSTK=SSCSTK, SSSTK=SSSTK)
 # endif
       end select
 

src/gotm/register_all_variables.F90

@@ -513,6 +513,7 @@ subroutine register_turbulence_variables(nlev)
    call fm%register('num', 'm2/s', 'turbulent diffusivity of momentum', standard_name='??', dimensions=(/id_dim_zi/), data1d=num(0:nlev),category='turbulence', part_of_state=.true.)
    call fm%register('nuh', 'm2/s', 'turbulent diffusivity of heat', standard_name='??', dimensions=(/id_dim_zi/), data1d=nuh(0:nlev),category='turbulence', part_of_state=.true.)
    call fm%register('nus', 'm2/s', 'turbulent diffusivity of salt', standard_name='??', dimensions=(/id_dim_zi/), data1d=nus(0:nlev),category='turbulence', part_of_state=.true.)
+   call fm%register('nucl', 'm2/s', 'turbulent diffusivity of momentum down Stokes gradient', standard_name='??', dimensions=(/id_dim_zi/), data1d=nucl(0:nlev),category='turbulence', part_of_state=.true.)
    call fm%register('gamu', 'm2/s2', 'non-local flux of u-momentum', standard_name='??', dimensions=(/id_dim_zi/), data1d=gamu(0:nlev),category='turbulence')
    call fm%register('gamv', 'm2/s2', 'non-local flux of v-momentum', standard_name='??', dimensions=(/id_dim_zi/), data1d=gamv(0:nlev),category='turbulence')
    call fm%register('gamh', 'K m/s', 'non-local heat flux', standard_name='??', dimensions=(/id_dim_zi/), data1d=gamh(0:nlev),category='turbulence')

src/meanflow/meanflow.F90

@@ -60,7 +60,7 @@ module meanflow
    REALTYPE, public, dimension(:), allocatable  :: SS,SSU,SSV
 
 !  Stokes-Eulerian cross-shear and Stokes shear squared
-   REALTYPE, public, dimension(:), allocatable  :: SSCSTK
+   REALTYPE, public, dimension(:), allocatable  :: SSCSTK, SSSTK
 
 !  buoyancy, short-wave radiation,
 !  extra production of tke by see-grass etc
@@ -364,6 +364,10 @@ subroutine post_init_meanflow(nlev,latitude)
    if (rc /= 0) STOP 'init_meanflow: Error allocating (SSCSTK)'
    SSCSTK = _ZERO_
 
+   allocate(SSSTK(0:nlev),stat=rc)
+   if (rc /= 0) STOP 'init_meanflow: Error allocating (SSSTK)'
+   SSSTK = _ZERO_
+
    allocate(xP(0:nlev),stat=rc)
    if (rc /= 0) STOP 'init_meanflow: Error allocating (xP)'
    xP = _ZERO_
@@ -464,6 +468,7 @@ subroutine clean_meanflow()
    if (allocated(SSU)) deallocate(SSU)
    if (allocated(SSV)) deallocate(SSV)
    if (allocated(SSCSTK)) deallocate(SSCSTK)
+   if (allocated(SSSTK)) deallocate(SSSTK)
    if (allocated(xP)) deallocate(xP)
    if (allocated(buoy)) deallocate(buoy)
    if (allocated(rad)) deallocate(rad)
@@ -523,6 +528,7 @@ subroutine print_state_meanflow()
    if (allocated(SSU)) LEVEL2 'SSU',SSU
    if (allocated(SSV)) LEVEL2 'SSV',SSV
    if (allocated(SSCSTK)) LEVEL2 'SSCSTK',SSCSTK
+   if (allocated(SSSTK)) LEVEL2 'SSSTK',SSSTK
    if (allocated(buoy)) LEVEL2 'buoy',buoy
    if (allocated(rad)) LEVEL2 'rad',rad
    if (allocated(xp))  LEVEL2 'xP',xP

src/meanflow/shear.F90

@@ -111,7 +111,7 @@ subroutine shear(nlev,cnpar)
 ! !USES:
    use meanflow,   only: h,u,v,uo,vo
    use meanflow,   only: SS,SSU,SSV
-   use meanflow,     only: SSCSTK
+   use meanflow,     only: SSCSTK, SSSTK
    use stokes_drift, only: dusdz, dvsdz
 
    IMPLICIT NONE
@@ -164,6 +164,9 @@ subroutine shear(nlev,cnpar)
       SSCSTK(i) =   dusdz%data(i) * (u(i+1)-u(i)) / (0.5*(h(i+1)+h(i))) &
                   + dvsdz%data(i) * (v(i+1)-v(i)) / (0.5*(h(i+1)+h(i)))
 
+      ! Stokes shear squared
+      SSSTK(i) = dusdz%data(i)**2 + dvsdz%data(i)**2
+
    end do
 
    SSU(0   ) = SSU(1    )
@@ -178,6 +181,9 @@ subroutine shear(nlev,cnpar)
    SSCSTK(0   ) = SSCSTK(1     )
    SSCSTK(nlev) = SSCSTK(nlev-1)
 
+   SSSTK (0   ) = SSSTK (1     )
+   SSSTK (nlev) = SSSTK (nlev-1)
+
    return
    end subroutine shear
 !EOC

src/meanflow/uequation.F90

@@ -5,7 +5,7 @@
 ! !ROUTINE: The U-momentum equation\label{sec:uequation}
 !
 ! !INTERFACE:
-   subroutine uequation(nlev,dt,cnpar,tx,num,gamu,ext_method)
+   subroutine uequation(nlev,dt,cnpar,tx,num, nucl, gamu,ext_method)
 !
 ! !DESCRIPTION:
 !  This subroutine computes the transport of momentum in
@@ -66,6 +66,7 @@ subroutine uequation(nlev,dt,cnpar,tx,num,gamu,ext_method)
    use meanflow,     only: gravity,avmolu
    use meanflow,     only: h,u,uo,v,w,avh
    use meanflow,     only: drag,SS,runtimeu
+   use stokes_drift, only: dusdz
    use observations, only: w_adv_input,w_adv_discr
    use observations, only: uprof_input,vel_relax_tau,vel_relax_ramp
    use observations, only: int_press_type
@@ -94,6 +95,9 @@ subroutine uequation(nlev,dt,cnpar,tx,num,gamu,ext_method)
 !  diffusivity of momentum (m^2/s)
    REALTYPE, intent(in)                :: num(0:nlev)
 
+!  eddy coefficient of momentum flux down the Stokes gradient (m^2/s)
+   REALTYPE, intent(in)                :: nucl(0:nlev)
+
 !  non-local flux of momentum (m^2/s^2)
    REALTYPE, intent(in)                :: gamu(0:nlev)
 
@@ -175,6 +179,10 @@ subroutine uequation(nlev,dt,cnpar,tx,num,gamu,ext_method)
 !     add non-local fluxes
 !      Qsour(i) = Qsour(i) - ( gamu(i) - gamu(i-1) )/h(i)
 
+!     add down Stokes gradient fluxes
+      Qsour(i) = Qsour(i) + ( nucl(i  )*dusdz%data(i  )                 &
+                             -nucl(i-1)*dusdz%data(i-1) )/h(i)
+
    end do
 
 !  implement bottom friction as source term

src/meanflow/vequation.F90

@@ -5,7 +5,7 @@
 ! !ROUTINE: The V-momentum equation\label{sec:vequation}
 !
 ! !INTERFACE:
-   subroutine vequation(nlev,dt,cnpar,ty,num,gamv,ext_method)
+   subroutine vequation(nlev,dt,cnpar,ty,num, nucl, gamv,ext_method)
 !
 ! !DESCRIPTION:
 !  This subroutine computes the transport of momentum in
@@ -45,6 +45,7 @@ subroutine vequation(nlev,dt,cnpar,ty,num,gamv,ext_method)
    use meanflow,     only: gravity,avmolu
    use meanflow,     only: h,v,vo,u,w,avh
    use meanflow,     only: drag,SS,runtimev
+   use stokes_drift, only: dvsdz
    use observations, only: w_adv_input,w_adv_discr
    use observations, only: vprof_input,vel_relax_tau,vel_relax_ramp
    use observations, only: int_press_type
@@ -73,6 +74,9 @@ subroutine vequation(nlev,dt,cnpar,ty,num,gamv,ext_method)
 !  diffusivity of momentum (m^2/s)
    REALTYPE, intent(in)                :: num(0:nlev)
 
+!  eddy coefficient of momentum flux down the Stokes gradient (m^2/s)
+   REALTYPE, intent(in)                :: nucl(0:nlev)
+
 !  non-local flux of momentum (m^2/s^2)
    REALTYPE, intent(in)                :: gamv(0:nlev)
 
@@ -154,6 +158,10 @@ subroutine vequation(nlev,dt,cnpar,ty,num,gamv,ext_method)
 !     add non-local fluxes
 !      Qsour(i) = Qsour(i) - ( gamv(i) - gamv(i-1) )/h(i)
 
+!     add down Stokes gradient fluxes
+      Qsour(i) = Qsour(i) + ( nucl(i  )*dvsdz%data(i  )                 &
+                             -nucl(i-1)*dvsdz%data(i-1) )/h(i)
+
    end do
 
 !  implement bottom friction as source term

src/stokes_drift/stokes_drift_theory.F90

@@ -45,8 +45,11 @@ subroutine stokes_drift_theory(nlev,z,zi,u10,v10,gravity)
    vs0%value = _ZERO_
    ds%value = _ZERO_
 
-!  wind direction
    wind_speed = sqrt(u10**2+v10**2)
+
+   if (wind_speed .gt. _ZERO_) then
+
+!  wind direction
    xcomp = u10 / wind_speed
    ycomp = v10 / wind_speed
 
@@ -66,6 +69,8 @@ subroutine stokes_drift_theory(nlev,z,zi,u10,v10,gravity)
    vs0%value = us0_to_u10 * v10
    ds%value = stokes_srf(0)*abs(zi(0))/max(SMALL, sqrt(us0%value**2.+vs0%value**2.))
 
+   end if
+
    return
 
    end subroutine stokes_drift_theory

src/turbulence/CMakeLists.txt

@@ -5,6 +5,7 @@ add_library(turbulence
    cmue_b.F90
    cmue_c.F90
    cmue_d.F90
+   cmue_d_h15.F90
    cmue_ma.F90
    cmue_sg.F90
    compute_cpsi3.F90

src/turbulence/alpha_mnb.F90

@@ -5,7 +5,7 @@
 ! !ROUTINE: Update dimensionless alpha's\label{sec:alpha}
 !
 ! !INTERFACE:
-   subroutine alpha_mnb(nlev,NN,SS)
+   subroutine alpha_mnb(nlev,NN,SS, SSCSTK, SSSTK)
 !
 ! !DESCRIPTION:
 ! This subroutine updates the dimensionless numbers $\alpha_M$, $\alpha_N$,
@@ -22,11 +22,18 @@ subroutine alpha_mnb(nlev,NN,SS)
 ! !USES:
   use turbulence,  only:     tke,eps,kb
   use turbulence,  only:     as,an,at
+  use turbulence,  only:     av, aw
   IMPLICIT NONE
 !
 ! !INPUT PARAMETERS:
   integer,  intent(in)      :: nlev
   REALTYPE, intent(in)      :: NN(0:nlev),SS(0:nlev)
+
+!  Stokes-Eulerian cross-shear (1/s^2)
+   REALTYPE, intent(in), optional      :: SSCSTK(0:nlev)
+
+!  Stokes shear squared (1/s^2)
+   REALTYPE, intent(in), optional      :: SSSTK (0:nlev)
 !
 ! !REVISION HISTORY:
 !  Original author(s): Lars Umlauf
@@ -35,21 +42,31 @@ subroutine alpha_mnb(nlev,NN,SS)
 !-----------------------------------------------------------------------
 ! !LOCAL VARIABLES:
   integer              :: i
-  REALTYPE             :: tau2
+  REALTYPE             :: tau2(0:nlev)
 
 !-----------------------------------------------------------------------
 !BOC
 
-  do i=0,nlev
-     tau2   = tke(i)*tke(i) / ( eps(i)*eps(i) )
-     as(i)  = tau2 * SS(i)
-     an(i)  = tau2 * NN(i)
-     at(i)  = tke(i)/eps(i) * kb(i)/eps(i)
+   do i=0,nlev
+      tau2(i) = tke(i)*tke(i) / ( eps(i)*eps(i) )
+      as(i)   = tau2(i) * SS(i)
+      an(i)   = tau2(i) * NN(i)
+      at(i)   = tke(i)/eps(i) * kb(i)/eps(i)
+
+!     clip negative values
+      as(i) = max(as(i),1.e-10*_ONE_)
+      at(i) = max(at(i),1.e-10*_ONE_)
+   end do
+
+   if (present(SSCSTK) .and. present(SSSTK)) then
+      do i=0,nlev
+         av(i)  = tau2(i) * SSCSTK(i)
+         aw(i)  = tau2(i) * SSSTK(i)
 
-!    clip negative values
-     as(i) = max(as(i),1.e-10*_ONE_)
-     at(i) = max(at(i),1.e-10*_ONE_)
-  end do
+!        clip negative values
+         aw(i) = max(aw(i),1.e-10*_ONE_)
+      end do
+   end if
 
   return
 end subroutine alpha_mnb