sin_DCCM2012 Function

pure elemental function sin_DCCM2012(spectrum,wspd,wdir,input_height,ustar,&
  vonkarman) result(tendency)

  !! TODO implement currents averaged over the effective depth layer for
  !! modulation of phase speed.

  !! Returns a spectrum instance with the wave growth ($S_{in}$) tendency
  !! formulated by Donelan et al. (2012) and based on the sheltering hypothesis
  !! by Jeffreys (1924, 1925).
  !!
  !! The result instance has the units of 1/s. This source function must be
  !! re-evaluated if any of the input parameters change.
  !!
  !! References:
  !!
  !! Donelan, M. A., M. Curcic, S. S. Chen, and A. K. Magnusson, 2012: Modeling
  !! waves and wind stress, *J. Geophys. Res. Oceans*, **117**, C00J23,
  !! doi:10.1029/2011JC007787.
  !!
  !! Jeffreys, H., 1924: On the formation of waves by wind, *Proc. R. Soc. A*,
  !! **107**, 189–206.
  !!
  !! Jeffreys, H., 1925: On the formation of waves by wind, II, *Proc. R. Soc.
  !! A*, **110**, 341–347.

  use mod_spectrum,only:spectrum_type
  use mod_aerodynamic_drag,only:windAtReferenceHeight
  use mod_const,only:twopi

  type(spectrum_type),intent(in) :: spectrum
    !! Input spectrum instance
  real(kind=rk),intent(in) :: wspd
    !! Input wind speed [m/s]
  real(kind=rk),intent(in) :: wdir
    !! Input wind direction [rad], mathematical convention
  real(kind=rk),intent(in) :: input_height
    !! Height of input wind speed [m/s]
  real(kind=rk),intent(in) :: ustar
    !! Air-side friction velocity [m/s]
  real(kind=rk),intent(in) :: vonkarman
    !! Von Karman constant

  type(spectrum_type) :: tendency

  real(kind=rk),dimension(:,:),allocatable :: s_in

  real(kind=rk),dimension(:),allocatable :: f
  real(kind=rk),dimension(:),allocatable :: th
  real(kind=rk),dimension(:),allocatable :: k
  real(kind=rk),dimension(:),allocatable :: cp
  real(kind=rk),dimension(:),allocatable :: omega
  real(kind=rk),dimension(:),allocatable :: half_wavelength
  real(kind=rk),dimension(:),allocatable :: wspd_input

  real(kind=rk) :: grav
  real(kind=rk) :: rho_air
  real(kind=rk) :: rho_water

  real(kind=rk),parameter :: a1_windsea = 0.11_rk
  real(kind=rk),parameter :: a1_swell = 0.01_rk
  real(kind=rk),parameter :: a1_opposed = 0.10_rk

  real(kind=rk),parameter :: field_scale_negative = a1_opposed / a1_windsea
  real(kind=rk),parameter :: field_scale_swell = a1_swell / a1_opposed

  real(kind=rk),dimension(:,:),allocatable :: sheltering_coefficient

  integer :: nfreq,nfreqs
  integer :: ndir,ndirs

  tendency = spectrum

  grav = spectrum % getGravity()
  rho_air = spectrum % getAirDensity()
  rho_water = spectrum % getWaterDensity()
  f = spectrum % getFrequency()
  th = spectrum % getDirections()
  k = spectrum % getWavenumber()
  cp = spectrum % getPhaseSpeed()
  half_wavelength = 0.5_rk*spectrum % getWavelength()
  omega = twopi*f

  nfreqs = size(f)
  ndirs = size(th)

  ! Evaluate wind speed at height of half-wavelength of each wave component
  wspd_input = windAtReferenceHeight(wspd,input_height,half_wavelength,ustar,&
    vonkarman)

  allocate(s_in(nfreqs,ndirs))
  allocate(sheltering_coefficient(nfreqs,ndirs))

  ! Set the initial sheltering coefficient to a1_windsea everywhere
  sheltering_coefficient = a1_windsea
  do concurrent(nfreq = 1:nfreqs,ndir = 1:ndirs)
    ! If input is negative, adjust the sheltering coefficient to a1_opposed
    if(wspd_input(nfreq)*cos(wdir-th(ndir))-cp(nfreq) < 0)then
      sheltering_coefficient(nfreq,ndir) = sheltering_coefficient(nfreq,ndir)&
        * field_scale_negative
      ! If input is negative but has positive misalignment, adjust the
      ! sheltering coefficient to a1_swell
      if(cos(wdir-th(ndir)) > 0)then
        sheltering_coefficient(nfreq,ndir) = sheltering_coefficient(nfreq,ndir)&
        * field_scale_swell
      endif
    endif
  enddo

  do concurrent(ndir=1:ndirs)
    s_in(:,ndir) = sheltering_coefficient(:,ndir)*rho_air/rho_water            &
      *(wspd_input*cos(wdir-th(ndir))-cp)*abs(wspd_input*cos(wdir-th(ndir))-cp)&
      *omega*k/grav
  enddo

  tendency = s_in
  deallocate(s_in,sheltering_coefficient)

endfunction sin_DCCM2012