!   This file contains subroutines for the Babanin physics according to Rogers et al (JTECH 2012)
!   based on work of Babanin, Young, Tsagareli, Ardhuin and others
MODULE SDSBABANIN

CONTAINS

  SUBROUTINE CALC_SDS(NFREQ,EDENS,F,KDS,ANAR_IN,TESTFL,KWAVE,CG)
    !
    USE SWCOMM1, ONLY: CHTIME
    USE SWCOMM3, ONLY: A1SDS,A2SDS,P1SDS,P2SDS,UPWARDS,GRAV,PI

    ! REAL   A1SDS  : coefficient on T1
    ! REAL   A2SDS  : coefficient on T2
    ! REAL   P1SDS  : power on T1
    ! REAL   P2SDS  : power on T2
    ! P1SDS,P2SDS is the power on (F-Fth)/Fth or (F-Fth)/F which is part of
    ! how I formulate "X" given by Young and Babanin eq 27 :
    ! P1SDS: power on threshold exceedence in T1 (precise definition depends
    ! on whether using "U" vs "D" method)
    ! P2SDS: power on threshold exceedence in T2 (precise definition depends on
    ! whether using "U" vs "D" method)
    ! LOGICAL UPWARDS : true if concave up

    IMPLICIT NONE
    !
!
!     SWAN (Simulating WAves Nearshore); a third generation wave model
!     Copyright (C) 1993-2017  Delft University of Technology
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation; either version 2 of
!     the License, or (at your option) any later version.
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     A copy of the GNU General Public License is available at
!     http://www.gnu.org/copyleft/gpl.html#SEC3
!     or by writing to the Free Software Foundation, Inc.,
!     59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
!
    !
    ! INPUT: EDENS(f), f, OUTPUT: calculate Kds(f) [Kds(f)=Sds(f)/EDENS(f)]
    ! EDENS is given wave spectrum m2/Hz
    ! f is frequency in Hz example: f=0.01:0.01:0.5
    ! @EDENS/@t=Sin+Snl+Sds
    !
    ! Kds=T1+T2
    !
    ! T1=a1*f*A*DEDENS
    !   (inherent breaking dissipation)
    ! T2=a2*cumsum(A*DEDENS)   an integration of A*DEDENS from fp to f
    !   (induced breaking dissipation)
    !
    ! A is directional narrowness, calculated in calling subroutine, but set
    !   to unity in this version
    !
    ! DEDENS=EDENS-EDENST
    ! EDENST=(2*g.^2)*((2*pi).^-4)*(f.^-5).*(A.^-1)*Bnt
    ! a1 and a2 are function of U/cp given by Babanin
    ! Bnt is an empirical coefficient related to the spectral density in the
    !   saturation range (Babanin et al 2007)
    !
    ! Purpose: to calculate the spectral dissipation Sds based on given wave
    !   spectrum
    ! References:
    ! 1) Babanin and Young ("WAVES" conf. 2005)
    ! 2) Young and Babanin (JPO 2006)
    ! 3) Babanin et al. (H/F conf. 2007), "Implementation of new ..."
    ! 4) WISE presentations by Rogers in 2009, 2010
    ! 5) Babanin et al. (JPO 2010,) paper re: Sds term
    ! 6) Tsagareli et al. (JPO 2010), paper re: Sin term
    ! 7) Rogers et al. (JTECH 2012)
    ! The code is "unified" insofar as depending on the setting of the
    !   logical "UPWARDS", one gets something similar to the method of :
    !   =true: Fabrice et al. (more similar to the early, nondirectional
    !      version of Fabrice's Sds published in a conference paper, as
    !      opposed to the Ardhuin et al. JPO 2010 directional Sds). With L=2
    !      and large exceedence, the dissipation will get very strong
    !     NDEDENS(is)=DEDENS(is)/EDENST(is)
    !   =false: Tsagareli, Babanin, Young. (but it differs from theirs
    !      insofar as our normalization for the cumulative term is more correct
    !      mathemetically, so we do not have to created a set of a1, b1 for each
    !      wave age...with the normalization used here, one calibration works
    !      for all  wave age values. Another improvement is in the exponent L :
    !      in the Tsagereli formula, the only dimensionally correct form of the
    !      equation is for L=1. This means that the normalized value is always
    !      less than 1, which means that using L=2 instead of L=1 will
    !      actually make the dissipation weaker (though since a1 a2 would need
    !      to be recalibrated, this relation is a bit more complicated)
    !
    ! Notes on calibration, last updated Sep 19 2013.
    ! Because these physics are "observation-based", and in Rogers et al (2012)
    !   we apply a number of rules and physical constraints based on
    !   observations, this means that the new physics are *less* tunable than
    !   prior source term packages, like KHH1984. This is a bit of a paradox,
    !   since we probably have more free parameters than KHH. Most important
    !   example, we can easily control SWH (total energy) but have little
    !   control over Tp, Tm02, Tm01, Tm-1,0
    ! Possibilities for calibrations are noted in the SWAN user's manual.
    !
    ! Further documentation with more up-to-date information than what you get
    !   from Rogers et al. (2012) :
    ! If you are inside the .navy.mil domain, see these wiki pages:
    ! https://www7320.nrlssc.navy.mil/Alvin/index.php/SWAN
    !   (especially the pages about * Babanin physics * Improvements * and
    !   * Suggested Settings * )
    ! If you are not inside the .navy.mil domain, you can ask us for
    !   this info: erick.rogers@nrlssc.navy.mil
    !
    ! Subroutine arguments:
    LOGICAL          , INTENT(IN)  ::  TESTFL
    ! TESTFL: true if output requested for this point (use sparingly, can result
    ! in large files)
    REAL             , INTENT(IN)  ::  EDENS(:)   ! E(f) m2/Hz
    REAL             , INTENT(IN)  ::  F(:)       ! frequency in Hz
    REAL             , INTENT(IN)  ::  ANAR_IN(:)
    ! ANAR_IN: directional narrowness as defined in Babanin publications
    ! (input value)
    INTEGER          , INTENT(IN)  ::  NFREQ      ! # freqs
    REAL             , INTENT(IN)  ::  KWAVE(:,:) ! KWAVE(MSC,MICMAX)
    REAL             , INTENT(IN)  ::  CG(:,:)    ! CG(MSC,MICMAX)
    REAL             , INTENT(OUT) ::  KDS(:)     ! Kds(f)=Sds(f)/E(f)

    ! Local variables:
    REAL             , ALLOCATABLE ::  SDS(:)     ! Sds(f), the source term
    REAL             , ALLOCATABLE ::  NDEDENS(:)
    ! NDEDENS(f)=DEDENS(f)/EDENST(f)
    REAL             , ALLOCATABLE ::  DEDENS(:)  ! DEDENS(f)=EDENS(f)-EDENST(f)
    REAL             , ALLOCATABLE ::  EDENST(:)  ! E(f) threshold for breaking
    REAL             , ALLOCATABLE ::  T1(:)      ! inherent dissipation/E(f)
    REAL             , ALLOCATABLE ::  T2(:)      ! induced dissipation/E(f)
    REAL             , ALLOCATABLE ::  ST1(:)     ! inherent dissipation
    REAL             , ALLOCATABLE ::  ST2(:)     ! induced dissipation
    REAL             , ALLOCATABLE ::  ANAR(:)
    ! ANAR = directional narrowness as defined in Babanin publications
    REAL             , ALLOCATABLE ::  XFF(:),ADF(:) ! temporary arrays
    REAL              :: ASUM   ! temporary variable for integration
    REAL              :: BNT
    ! BNT is an empirical coefficient related to the spectral density in the
    ! saturation range (Babanin et al 2007)
    INTEGER           :: II,IS ! counters
    ! *_INT: integrated values (for test output only)
    REAL              :: ST1_INT,ST2_INT,SDS_INT
    ! IMAX,I3FP, fp,  FD(:) re: T1 and T2 at 3fp (for test output only)
    INTEGER           :: IMAX,I3FP
    REAL              :: FP
    REAL, ALLOCATABLE :: FD(:)
    REAL ::  ELIM  ! needed for UPWARDS=.FALSE.
    REAL ::  CTMP1 ! temporary variable

    ! ------- START SUBROUTINE -------------------------------------------------

    ALLOCATE(SDS(NFREQ))
    ALLOCATE(XFF(NFREQ))
    ALLOCATE(NDEDENS(NFREQ))
    ALLOCATE(DEDENS(NFREQ))
    ALLOCATE(EDENST(NFREQ))
    ALLOCATE(T1(NFREQ))
    ALLOCATE(T2(NFREQ))
    ALLOCATE(ANAR(NFREQ))
    ALLOCATE(ST1(NFREQ))
    ALLOCATE(ST2(NFREQ))
    ALLOCATE(ADF(NFREQ))
    ALLOCATE(FD(NFREQ))

    BNT=(0.035**2)    !  Bnt value given by Babanin et al (2007)

    ! --------------------------------------------------------------------------
    ! get a1 and a2 as a function of U/cp
    ! --------------------------------------------------------------------------

    ! needed: A(f), see Young and Babanin eq 19.
    ! Originally, it was read in and applied, but per recommendation by Alex
    ! that A(f) is unnecessary/obsolete, I set it to 1.0 here:
    ANAR=1.0 ! option 1
!   ANAR=ANAR_IN ! option 2

    ! (point output write location 1)
    !IF(TESTFL)WRITE(*,*)'calc_Sds : a1,a2 = ',A1SDS,A2SDS

!   new calculation, depends on k and Cg instead of f :
    CTMP1=2.0*PI*BNT
    DO IS=1,NFREQ
       EDENST(IS)=CTMP1/(CG(IS,1)*KWAVE(IS,1)**3)
    END DO
!   potential optimization: if ANAR is always unity, then EDENST does not vary
!     with time step and only needs to be calculated once, rather than every
!     time this routine is called.

    DEDENS=EDENS-EDENST ! matrix operation

    ! if DEDENS < 0 set to zero
    DEDENS=MAX(0.0,DEDENS)

! notes: Mar 22 2011, Stefan has noticed sensitivity to ELIM.  (this variable
!   only applies to concave down case). I have experimented with ELIM=0.0 and
!   noticed no sensitivity for the U10=12 m/s point model case.

! ELIM is needed for "concave down"
!   ELIM=maxval(Edens)*1.0e-5 ! option 1
    ELIM=0.0 ! option 2

!   See notes above re: UPWARDS

    IF(UPWARDS)THEN
       NDEDENS(1:NFREQ)=DEDENS(1:NFREQ)/EDENST(1:NFREQ)
    ELSE
       DO IS=1,NFREQ
          IF(EDENS(IS).GT.ELIM)THEN
             NDEDENS(IS)=DEDENS(IS)/EDENS(IS)
          ELSE
             NDEDENS(IS)=0.0
          END IF
       ENDDO
    ENDIF

    NDEDENS=MAX(0.0,NDEDENS)
    ! --------------------------------------------------------------------------
    !                      calculate T1
    !  -------------------------------------------------------------------------

    DO  IS=1,NFREQ
       T1(IS)=A1SDS*F(IS)*ANAR(IS)*NDEDENS(IS)**P1SDS
    END DO

    ! --------------------------------------------------------------------------
    !                      calculate  T2  an integration from fp to f
    ! --------------------------------------------------------------------------

    ! Note that we do not worry about starting at fp (which can be difficult to
    !  define), since stuff below fp is typically not breaking, so it is not
    !  part of the calculation anyway.

    IF(A2SDS.GT.0.)THEN
      XFF=0.0
      DO  IS=1,NFREQ
         ASUM=0.0
         DO II=1,IS
            XFF(II)=F(II)
            ADF(II)=ANAR(II)*NDEDENS(II)**P2SDS
         ENDDO
         ! THE "INTEGRATE" ROUTINE IS BASED ON FINITE DIFFERENCING, BUT WE COULD
         ! REPLACE WITH AN OPERATION THAT USES FRINTF
         CALL INTEGRATE(ASUM,XFF,ADF,IS)
         T2(IS)=A2SDS*ASUM
      ENDDO
    ELSE
      T2=0.
    ENDIF

    T1=MAX(0.0,T1)
    IF(A2SDS.GT.0.)THEN
      T2=MAX(0.0,T2)
      KDS=(T1+T2)
    ELSE
      KDS=T1
    ENDIF

    ! (test output for NRL purposes)
!NRL    IF(TESTFL)THEN
!NRL       WRITE(411,*)CHTIME,' % CHTIME'
!NRL       WRITE(412,*)'% f(is),EDENS(is),EDENST(is),T1(is),T2(is)'
!NRL       DO  IS=1,NFREQ
!NRL          WRITE(412,205)F(IS),EDENS(IS),EDENST(IS),T1(IS),T2(IS)
!NRL       END DO
!NRL205    FORMAT(5(1X,E11.5))
!NRL
!NRL       ! calculate integrated T1 and T2 (for output purposes only)
!NRL
!NRL       ! Since we have already non-dimensionalized by using NDEDENS instead of
!NRL       !   DEDENS, we do it backwards: Sds=Kds*Edens.
!NRL       !   Note that in most cases, the calling routine will not use Sds.
!NRL       DO IS=1,NFREQ
!NRL          SDS(IS)=KDS(IS)*EDENS(IS)
!NRL          ST1(IS)=T1(IS)*EDENS(IS)
!NRL          ST2(IS)=T2(IS)*EDENS(IS)
!NRL       END DO
!NRL
!NRL       CALL INTEGRATE(ST1_INT,F,ST1,NFREQ)
!NRL       CALL INTEGRATE(ST2_INT,F,ST2,NFREQ)
!NRL       CALL INTEGRATE(SDS_INT,F,SDS,NFREQ)
!NRL
!NRL       ! calculate T1 and T2 at 3fp
!NRL       IMAX=MAXLOC(EDENS,1)
!NRL       FP=F(IMAX)
!NRL       FD=ABS(F-FP*3.0)
!NRL       I3FP=MINLOC(FD,1)
!NRL
!NRL       ! (POINT OUTPUT WRITE LOCATION 3)
!NRL       WRITE(*,208)ST1_INT,ST2_INT,SDS_INT,ST1(I3FP),ST2(I3FP)
!NRL208    FORMAT('integral of T1,T2,Sds = ',5(1X,E14.8))
!NRL    ENDIF

! Deallocate and return
    DEALLOCATE(SDS)
    DEALLOCATE(XFF)
    DEALLOCATE(NDEDENS)
    DEALLOCATE(DEDENS)
    DEALLOCATE(EDENST)
    DEALLOCATE(T1)
    DEALLOCATE(T2)
    DEALLOCATE(ANAR)
    DEALLOCATE(ST1)
    DEALLOCATE(ST2)
    DEALLOCATE(ADF)
    DEALLOCATE(FD)

  END SUBROUTINE CALC_SDS

  !****************************************************************************
  SUBROUTINE SWIND_DBYB ( SPCSIG,THETAW,KWAVE,MEMSINA,MEMSINB, &
                          AC2,UFRIC,WIND10,SPCDIR,ANYWND,CG    &
                         ,ZELEN )
  !****************************************************************************

    USE SWCOMM1, ONLY: CHTIME
    USE SWCOMM3
    USE SWCOMM4  ! includes TESTFL
    USE OCPCOMM4
!ESMF    USE M_GENARR, ONLY: SAVE_SINBAC, SINBAC

    IMPLICIT NONE
!
!
!     SWAN (Simulating WAves Nearshore); a third generation wave model
!     Copyright (C) 1993-2017  Delft University of Technology
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation; either version 2 of
!     the License, or (at your option) any later version.
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     A copy of the GNU General Public License is available at
!     http://www.gnu.org/copyleft/gpl.html#SEC3
!     or by writing to the Free Software Foundation, Inc.,
!     59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
!
    !
    !  SUBROUTINE SWIND_DBYB follows the structure of SWIND3 of SWAN.
    !     As such, it contains features that are not general to all wave models.
    !     In particular, the quadrant sweeping creates special challenges that
    !     are not general.
    !
    !  0. Authors
    !
    !     W.E. Rogers, NRL 7320
    !
    !  1. Updates
    !
    !  2. PURPOSE
    !
    !     Computation of the source term for the wind input for a
    !     third generation wind growth model:
    !
    !     1)  Exponential input term, (Donelan et al. JPO 2006)
    !
    !  3. METHOD
    !
    !     Sin (s,d) =  B*E(s,d)
    !
    !  4. Argument variables
    !
    ! i   SPCDIR: (*,1); spectral directions (radians)
    !             (*,2); cosine of spectral directions
    !             (*,3); sine of spectral directions
    !             (*,4); cosine^2 of spectral directions
    !             (*,5); cosine*sine of spectral directions
    !             (*,6); sine^2 of spectral directions
    ! i   SPCSIG: Relative frequencies in computational domain in sigma-space
    !
    !        IS          Counter of relative frequency band
    !        ID          Counter of directional distribution
    !        MSC         Maximum counter of relative frequency
    !        MDC         Maximum counter of directional distribution
    !
    !        REALS:
    !        ---------
    !        THETA       Spectral direction
    !        THETAW      Mean direction of the relative wind vector
    !        UFRIC       Wind friction velocity
    !
    !        one and more dimensional arrays:
    !        ---------------------------------
    !        KWAVE     2D    Wavenumber
    !        PWIND     1D    Wind coefficients
    !        ANYWND    1D    Wind input for bin considered
    !
    !     5. SUBROUTINES CALLING
    !
    !        SOURCE
    !
    !     6. SUBROUTINES USED
    !
    !        calc_Lfactor
    !
    !     7. ERROR MESSAGES
    !
    !        ---
    !
    !     8. REMARKS
    !
    !   Here is what I had prior to Oct 4 2011:
    !   U10PROXY=28.0*UFRIC
    !
    !   U10=28U* implies a fixed Cd. In order for the model to scale with U*,
    !   we must use a fixed Cd to convert the Donelan formula from U10 to U*.
    !   If  we use a proper, realistic Cd to do this, then we are basically
    !   keeping the Donelan formula in terms of U10. The 28U* comes from Komen
    !   et al. (1984). The physical constraint, which is in terms of U*, make
    !   the resulting scaling behavior less obvious.
    !
    !   Updated Aug 9 2012: When the new logical variable "TRUE_U10" is true,
    !   SWAN uses U10 instead of 28Ustar in SWIND_DBYB.  Ustar is used only
    !   for the physical constraint in this case.
    !
    !   Updated Feb 7 2014: "28" in the formula U10PROXY=28*UFRIC is made a
    !   user-defined variable. Note that "28" corresponds to C10=1.28e-3,
    !   which one might see for U10~5 m/s (pertinent to Snyder measurements
    !   and therefore used by Komen et al. (1984)) whereas a typical Lake
    !   George / DBYB value is C10=1.41e-3 to 1.43e-3. So U10=26.5*Ustar,
    !   and the max value of Cd would give U10~23.5*Ustar. Or, if we go by
    !   wind speed, a typical U10 during AUSWEX was 11 m/s, so with Hwang Cd,
    !   that comes out to a factor of 24 to 25. Via numerical experiments, I
    !   have found that using a factor (which I call "WINDSCALING") smaller
    !   than 28, and recalibrating the model, gives a stronger tendency to
    !   overpredict high frequency energy (already a problem mentioned in
    !   Rogers et al. JTECH 2012). Using WINDSCALING>28 and recalibrating
    !   the model improves the problem with overprediction of high frequency
    !   energy, but unfortunately, is making our manipulation of the DBYB
    !   formula more severe (noting that if we want to use DBYB without
    !   manipulation, we should use the "TRUE_U10" setting). Using larger factor
    !   also tends to reduce waveheights during tropical cyclones, e.g.
    !   maximum SWH during my Ivan simulation at 42040 is reduced by ~60 cm,
    !   which is a slight improvement. Suprisingly, using the larger factor
    !   does not make the "tail reduction" physical constraint more active.
    !   This is because energy in the tail is reduced, so Sin in the tail is
    !   reduced, so "first guess" of tau_wave is not significantly increased.
    !
    !   To summarize, the two options are as follows:
    !       a) U10PROXY=WINDSCALING*USTAR (recommended)
    !          (default WINDSCALING=28)
    !       b) U10=U10 (activated by using "TRUE_U10" setting)
    !
    !   Additional remarks about the tendency to overpredict high frequency
    !   energy (particularly with WINDSCALING=28) :
    !   This is evidenced by tendency to underpredict Tm02 even when
    !   SWH is overpredicted (in simulations of windsea-only situations). It
    !   can also be seen in direct comparison of E(f), of course. Another clear
    !   indicator is the tendency of E/ET to exceed 5, for example at 0.3 Hz
    !   when U10~12 m/s (see my slides from NRL External Review July 2009).
    !   E/ET, according to Babanin and Soloviev (1998) should not exceed 5, and
    !   this tendency of the "real ocean" is confirmed from buoy analyses by
    !   DW in these slides. Using larger values of p1 p2 also improves the
    !   overprediction, but not much. For example changing [p1,p2] from [4,4]
    !   to [8,8] improves the problem less than changing WINDSCALING from 28
    !   32 (with recalibration). Essentially, the Sds is like a spring that pulls
    !   E(f) toward ET(f). The values a1 a2 p1 p2 control the tension in the
    !   spring. Going from 28 to 32 requires a recalibration, an increase in
    !   a1 a2, i.e. an increase in the tension of the spring. Changing [p1,p2]
    !   from [4,4] to [8,8] is another type of increase in the tension of the
    !   spring. Another piece of evidence of overprediction of high frequency
    !   energy: it can be shown that the slope of the high-frequency tail in
    !   the saturation region is between -4 and -5, where we expect it to be -5.
    !
    !   Notes on method of calculating AINV:
    !
    !    1) Using single point model, U10=12 m/s
    !      CIRCLE 36 0.0271 10.0 62
    !      GEN3 BABANIN 6.24E-7 8.74E-6 4.0 4.0 1.2 0.0000 UP VECTAU AGROW
    !
    !      Using Babanin AINV as before
    !      max(HSKn1),max(HSKn2),max(HS) =  2.29 ;  2.77 ;  2.55
    !
    !    2) Using        AINV=4.0112*((RMSDIR(IS)-0.2) ** (0.81)) ! ER FORMULA
    !      (used in swancom2 in 2009, but precise origin unknown)
    !      (tends to give a high estimate of AINV)
    !      max(HSKn1),max(HSKn2),max(HS) =  2.29 ;  2.77 ;  2.52
    !
    !    3) Using        AINV=(RMSDIR(IS)/0.38) ** (1.0/0.94) ;   ! DW FORMULA
    !      (used in swancom2 in 2009, but precise origin unknown)
    !      max(HSKn1),max(HSKn2),max(HS) =  2.29 ;  2.77 ;  2.60
    !
    !    4) Using       AINV=(RMSDIR(IS)*3.4)-0.25 ! eyeball from Kuik_vs_AINV.m
    !      max(HSKn1),max(HSKn2),max(HS) =  2.29 ;  2.77 ;  2.55
    !
    !    5) Using       AINV=(RMSDIR(IS)*3.5333)-0.34862 ! robustfit from
    !      Kuik_vs_AINV.m
    !      max(HSKn1),max(HSKn2),max(HS) =  2.29 ;  2.77 ;  2.55
    !
    !    6) Using        AINV=RMSDIR(IS)*2.51 ! from experiments with
    !      test_DSPR_normal_sigma.m
    !      (tends to give a low estimate of AINV)
    !      max(HSKn1),max(HSKn2),max(HS) =  2.29 ;  2.77 ;  2.63
    !
    !    10. SOURCE
    !
    !***********************************************************************
    !
    ! Subroutine arguments:
    REAL    SPCDIR(MDC,6)
    REAL    SPCSIG(MSC)
    REAL    THETAW, KWAVE(MSC,ICMAX)
    REAL    MEMSINA(MDC,MSC,MCGRD), MEMSINB(MDC,MSC,MCGRD)
    REAL    UFRIC ,AC2(MDC,MSC,MCGRD) , WIND10
    REAL    CG(MSC,ICMAX)
    REAL    ZELEN(MCGRD)
    LOGICAL ANYWND(MDC)

    ! Local variables:
    INTEGER IDDUM ,ID    ,IS
    REAL    THETA ,SIGMA    ,SWINEB, CTW   ,STW   ,COSDIF
    REAL    TEMP2 ,UoverC
    REAL    S_IN(MDC,MSC)
    REAL    DMAX,AINV
    ! Ktheta is like D(theta), except max value at each freq is unity
    REAL    KTHETA(MSC,MDC)
    REAL    ANAR(MSC),SIGDENS(MSC),SQRTBN(MSC),CINV(MSC)
    REAL    GAMMAD,GDONEL,TEMP4,WPSI,TEMP5,TEMP6,BN
    REAL    SIN1D(MSC),SWND,FREQ,STRESS
    INTEGER IENT
    REAL    RMSDIR(MSC)

    REAL EDENS2D ! for test calcs
    REAL CINV2,CTH,STH
    REAL TAUX_linear,TAUY_linear
!NRL    REAL TAUX_tmp,TAUY_tmp,ENCHECK
    REAL ZE   !roughness length

    SAVE IENT
    DATA IENT/0/
    IF (LTRACE) CALL STRACE (IENT,'SWIND_DBYB')

    CTW   = COS(THETAW)
    STW   = SIN(THETAW)

!   Calculate 1d spectrum E(sigma)
    DO  IS = 1, MSC
       SIGDENS(IS) = 0.
       DO  ID = 1, MDC
          SIGDENS(IS) = SIGDENS(IS) + SPCSIG(IS) * AC2(ID,IS,KCGRD(1))
          KTHETA(IS,ID)=AC2(ID,IS,KCGRD(1))
       END DO
       SIGDENS(IS)=SIGDENS(IS)*DDIR  ! units m^2/(radHz)
    END DO

!   Calculate K(sigma,theta)
    DO  IS = 1, MSC
       DMAX=-1.0
       DO  ID = 1, MDC
          IF(KTHETA(IS,ID).GT.DMAX)DMAX=KTHETA(IS,ID)
       END DO
       IF(DMAX.EQ.0.0)THEN
          ! a fix for freq bins (usually first two or so) that are empty
          DMAX=1.0
          DO  ID = 1, MDC
             KTHETA(IS,ID)=1.0
          END DO
       ENDIF

!      Optional: Calculate circular RMS directional spread (for comparison with
!      AINV) :
!      rmsdir(is)=kuik(ktheta(is,:),ddir,spcdir(:,2),spcdir(:,3),.true./.false)

       DO  ID = 1, MDC
          KTHETA(IS,ID)=KTHETA(IS,ID)/DMAX
       END DO
    END DO

!   Calculate SQRT(Bn(IS))
    DO IS=1,MSC
       AINV=0.0
       DO  ID = 1, MDC
          AINV=AINV+KTHETA(IS,ID)*DDIR
       END DO

!   Optional: use RMSDIR in place of AINV: to see how to do this, read notes in
!      rev 622 and prior

       ANAR(IS)=1.0/AINV
       BN=ANAR(IS)*SIGDENS(IS)*(KWAVE(IS,1)**3)*CG(IS,1)
       SQRTBN(IS)=SQRT(BN)

    END DO

! Calculate S_in for *ALL* directions (not just current sweep). Though this
!   is somewhat wasteful (S_in will end up being calculated 4 times for
!   each bin), it is nececessary in order to have the correct stress going
!   into CALC_LFACTOR
! Update May 2017: Above issue was addressed by Marcel Z. via "memsin".

! Note, if this setting for TEMP2 is modified, it should also be modified in
! subroutine "CALC_TAU_TOTAL"
    IF(TRUE_U10)THEN
       TEMP2 = WIND10  !  "TEMP2=WIND10" is basically saying "use DBYB as is"
    ELSE
       TEMP2 = WNDSCL * UFRIC
    ENDIF
    S_IN=0.0

    DO IS = 1, MSC
       SIGMA = SPCSIG(IS)
       CINV(IS)  = KWAVE(IS,1) / SIGMA
       UOVERC = TEMP2 * CINV(IS)
       DO ID=1,MDC
          IF ( ANYWND(ID) ) THEN
             COSDIF = SPCDIR(ID,2)*CTW + SPCDIR(ID,3)*STW

!-------------------------------------------------------------------------------
! Yalin Fan, 08/06/13: Below equations calculate W(f,theta) in eq6 of Rogers et al (2012)
! This equation limits the grow rate to be positve, while Tolman and Chalikov (1996)
! gives negative growth rate when wind and wave have large angle
! If the 0 constraint is removed and the sign of TEMP4 is maintained, then W(f,theta)
! is very similar to Beta in Tolman and Chalikov (1996)
!-------------------------------------------------------------------------------

             TEMP4=( UOVERC * COSDIF - 1.0 )
             TEMP4=MAX(0.0,TEMP4)
             WPSI=TEMP4**2
             TEMP5=10.0*SQRTBN(IS)*WPSI-11.0
             TEMP6=1.0+TANH(TEMP5)
             GDONEL=2.8-TEMP6
             GAMMAD=GDONEL*SQRTBN(IS)*WPSI
             SWINEB = GAMMAD * SIGMA * PWIND(9) ! NEW SWINEB

! In Donelan notation, SWINEB is BETA, so I use BETA=GAMMA*sigma*rhoa/rhow
!   (Donelan eq 3)

! Calculate actual wind input term Sin=Beta*Edens

             S_IN(ID,IS)=SWINEB*AC2(ID,IS,KCGRD(1))*SPCSIG(IS)

! Note that IMATRA calculation will be done after reduction operation

          END IF
       ENDDO
    ENDDO

! calculate stress for the linear wind input part
    TAUX_linear=0.0
    TAUY_linear=0.0
    DO  IS = 1, MSC
        CINV2=KWAVE(IS,1)/SPCSIG(IS)
        DO ID = 1, MDC
           CTH = SPCDIR(ID,2) ! cos(theta)
           STH = SPCDIR(ID,3) ! sin(theta)
           TAUX_linear=TAUX_linear+CTH*MEMSINA(ID,IS,KCGRD(1))*CINV2* SPCSIG(IS)**2*FRINTF*DDIR
           TAUY_linear=TAUY_linear+STH*MEMSINA(ID,IS,KCGRD(1))*CINV2* SPCSIG(IS)**2*FRINTF*DDIR
        END DO
    END DO
    TAUX_linear=TAUX_linear*PWIND(17)*GRAV
    TAUY_linear=TAUY_linear*PWIND(17)*GRAV
! Note: test output that was here has been deleted. See rev 622 or earlier.

    CALL CALC_LFACTOR(TAUX_linear,TAUY_linear,RDFSIN,S_IN,UFRIC,PWIND,DDIR,SPCSIG,FRINTF,CINV,GRAV, &
                      WIND10,TESTFL,SPCDIR,VECTOR_TAU,TRUE_U10,CTW,STW,ZE)

    ZELEN(KCGRD(1)) = ZE

!-------------------------------------------------------------------------------
!   Begin negative wind input
!-------------------------------------------------------------------------------

!-------------------------------------------------------------------------------
! Yalin, Aug30, 2013.    Add negative input when wind and wave have large angle
! RDCOEF is the user adjustable reduction coefficient for the negative energy.
! RDCOEF = 0.4 is the suggested value. It is specified in INPUT as "NEGATINP 0.4"
! If no keyword "NEGATINP" specified, RDCOEF is default to 0.0
! ------------------------------------------------------------------------------
! WER, June 8, 2017. "0.4" is from Donelan (1999, Book Chapter, "Wind-induced growth
! and attenuation of laboratory waves") (See also discussion in Donelan et al. (2006).)
! Zieger et al. (2015) use reduction_coef = 0.04 in conjunction with the WW3 v4.18
! swell dissipation (see their table 1). (At time of writing, the WW3 v5.08 swell
! dissipation mentioned in the same table is not implemented in SWAN)
! ------------------------------------------------------------------------------
! WER, June 9, 2017. For efficiency, only go into this loop if ZIEGER. (Changed)
! ------------------------------------------------------------------------------
! WER, Jun 12 2017: I confirmed that use of ANYWND here is a bug, since it disables
! negative wind input for case of component in the +/- 90 deg arc opposing the wind.
! I have removed the ANYWND operation.

    IF(ZIEGER)THEN

      DO IS = 1, MSC
         SIGMA = SPCSIG(IS)
         CINV(IS)  = KWAVE(IS,1) / SIGMA
         UoverC = TEMP2 * CINV(IS)
         DO ID=1,MDC
            COSDIF = SPCDIR(ID,2)*CTW + SPCDIR(ID,3)*STW
            TEMP4 = ( UoverC * COSDIF - 1.0 )                         ! modify Rogers_etal_2012 eq 6
            TEMP4 = MIN(0.0,TEMP4)                                    ! modify Rogers_etal_2012 eq 6
            WPSI = TEMP4**2                                           ! modify Rogers_etal_2012 eq 6
            TEMP5 = 10.0*SQRTBN(IS)*WPSI-11.0
            TEMP6 = 1.0+TANH(TEMP5)
            GDONEL = 2.8-TEMP6
            GAMMAD = GDONEL*SQRTBN(IS)*WPSI
            SWINEB = GAMMAD * SIGMA * PWIND(9)
            S_IN(ID,IS)= S_IN(ID,IS) - SWINEB*AC2(ID,IS,KCGRD(1))*SPCSIG(IS)*RDCOEF
         ENDDO
      ENDDO

    END IF

!-------------------------------------------------------------------------------
!   End negative wind input
!-------------------------------------------------------------------------------

! We want to add B*N to RHS. This is SWINEB*AC2, see SWIND3.
! Above, we had S_IN(ID,IS)=SWINEB*AC2(ID,IS,KCGRD(1))*SPCSIG(IS)
! Thus, we just need to divide out SPCSIG(IS).
! Store this wind input in MEMSINB array for every grid point,
! which can be filled in IMATRA in next four sweeps
! (see subroutine FILSIN).
! Also note that this part, B*N, is explicit and stored in PLWNDS,
! which represents total wind input (=A+BE)
! (see subroutine FILSIN)

    DO IS = 1, MSC
       DO ID = 1, MDC
          S_IN(ID,IS)= S_IN(ID,IS) * RDFSIN(IS)
          MEMSINB(ID,IS,KCGRD(1)) = S_IN(ID,IS) / SPCSIG(IS)
!ESMF          IF (SAVE_SINBAC) SINBAC(ID,IS,KCGRD(1)) = &
!ESMF             SINBAC(ID,IS,KCGRD(1)) + S_IN(ID,IS)/SPCSIG(IS)
       ENDDO
    ENDDO

! calculate stress (test point only)
!NRL      IF(TESTFL)THEN
!NRL         TAUX_tmp=0.0
!NRL         TAUY_tmp=0.0
!NRL         ENCHECK=0.0
!NRL         DO  IS = 1, MSC
!NRL            CINV2=KWAVE(IS,1)/SPCSIG(IS)
!NRL            DO ID = 1, MDC
!NRL               CTH = SPCDIR(ID,2) ! cos(theta)
!NRL               STH = SPCDIR(ID,3) ! sin(theta)
!NRL               TAUX_tmp=TAUX_tmp+CTH*CINV2*S_IN(ID,IS)*DDIR*FRINTF*SPCSIG(IS)
!NRL               TAUY_tmp=TAUY_tmp+STH*CINV2*S_IN(ID,IS)*DDIR*FRINTF*SPCSIG(IS)
!NRL               ENCHECK=ENCHECK+AC2(ID,IS,KCGRD(1))*DDIR*FRINTF*SPCSIG(IS)**2
!NRL            END DO
!NRL         END DO
!NRL         TAUX_tmp=TAUX_tmp*PWIND(17)*GRAV
!NRL         TAUY_tmp=TAUY_tmp*PWIND(17)*GRAV
!NRL         WRITE(*,*)'SWIND: TAUX,TAUY,HM0 = ',TAUX_tmp,TAUY_tmp,(4*SQRT(ENCHECK))
!NRL      ENDIF

    !     *** test output ***
    IF (ITEST.GE. 80.AND.TESTFL) THEN
       WRITE(PRTEST,6000) KCGRD(1), THETAW*180./PI
6000   FORMAT(' SWIND_DBYB: POINT  THETAW        :',I5,E12.4)
       WRITE(PRTEST,6100) TEMP2, UFRIC
6100   FORMAT(' SWIND_DBYB: TEMP2 UFRC     :',3E12.4, /, &
         '  IS ID1 ID2       Wind source term')
       DO IS = 1, MSC
          WRITE(PRTEST,6200) IS, 1, MDC,(MEMSINB(ID,IS,KCGRD(1)), &
             ID=1,MDC)
6200      FORMAT(3I4, 600e12.4)
       ENDDO
       WRITE(PRTEST,*)
    END IF

    RETURN
  END SUBROUTINE SWIND_DBYB

  SUBROUTINE CALC_LFACTOR(TAUX_linear,TAUY_linear,LFACTOR_L,S_IN,UFRIC,PWIND,DDIR_RAD,SIGMA_S,FRINTF, &
                          CINV_S,GRAV,WIND10,TESTFL,SPCDIR,VECTOR_TAU,TRUE_U10,CTHETA_WIND, &
                          STHETA_WIND, ZE)
!
    USE SWCOMM1, ONLY: CHTIME
!
    IMPLICIT NONE
!
!
!     SWAN (Simulating WAves Nearshore); a third generation wave model
!     Copyright (C) 1993-2017  Delft University of Technology
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation; either version 2 of
!     the License, or (at your option) any later version.
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     A copy of the GNU General Public License is available at
!     http://www.gnu.org/copyleft/gpl.html#SEC3
!     or by writing to the Free Software Foundation, Inc.,
!     59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
!
!
!   Objective : provide Lfactor (1-d array) by solving for optimal REDUC
!   (scalar). It uses a solver method that should work for all monotonic
!   relations. For non-monotonic relations, it may get confused. Here,
!   tau_normal definitely has a monotonic dependence on REDUC, so this
!   potential limitation on the solver is OK. (Monotonicity of S_in, or lack
!   thereof, is irrelevant.)

!   Points made here:
!   1) The normal stress plus the tangential/viscous stress cannot be greater
!      than the total stress.
!   2) We do not know what the contribution is beyond fmax, but we know that it
!      is not negative, so again, we say that the normal stress in our
!      prognostic range cannot be greater than the total stress...so
!      tau_total=tau_normal_prognostic+tau_normal_diagnostic+tau_tangential...
!      all are >=0, so maximum allowed value of tau_normal_prognostic is
!      tau_total.  Update: in this code, we partially address this issue by
!      extending the integration (fmax) to 10 Hz for the purposes of this
!      integration.
!   3) Similar to Kakha's argument, we use the exponential adjustment that is
!      stronger at the higher frequencies (more reduction) since this is where
!      our formula is less well-informed by observations (more wiggle room),
!      since typical observations are for the dominant waves, f=fp
!   4) Kakha uses f/fp to scale the reduction, which makes sense in the context
!      of item (3), but use of fp has disadvantages, so we use U/C instead.
!      This means that for young wind sea, S_in for the entire spectrum is
!      reduced...will this be a problem? It is probably worthwhile to try using
!      an fp value defined in a manner similar to Tolman and Chalikov, which
!      does not have the traditional problems with using fp.
!   5) S_in=S_in_linear+S_in_exponential....same argument as (1). Linear wind
!      input is assumed small and not included in the stress calculation. This
!      was in fact a major source of error in some cases, but I have limited
!      the linear wind input source function so that it cannot contribute more
!      than 1% of the total stress now. Thus, the error associated with
!      excluding it from the stress calculation here will be less than 1% now.
!      Having said that, it would be useful to include it, as a future code
!      improvement. (Update: Yalin has added linear part to stress calculation.
!      This was rev 669, Feb 26 2014.)

!    Notes re: "vector_tau":
!      In earlier versions of this code, to simplify calcs, I used S_in1D(f)
!      S_in1D_S=sum(S_in,1)*DDIR_rad*(2.0*PI) ! in units m2/Hz
!      this simplication implies that we will follow the Tsagareli et al.
!      (2010 JPO) method for calculating stress:
!      tau=integral (Sin(f)/C(f) df
!      This calculation is not so good, since it assumes all stresses are in
!      the same direction and so never counteract each other.
!      That is the method reported in Rogers et al. (JTECH 2012)
!    Updated August 2012:
!      We replace the non-directional
!      tau=integral of Sin(f)/C(f) df
!      with
!      tau_x=integral of Sin(f,theta)*k_x/(k*C(f)) df dtheta
!      tau_y=integral of Sin(f,theta)*k_y/(k*C(f)) df dtheta
!      essentially using the unit vector <k> = k_vector/k_scalar
!      and it can be written also
!      tau_vector=integral of (k_vector/omega) Sin(f,theta) df dtheta
!      Now, if vector_tau is true, we use a vector calculation.
!      If vector_tau is false, we do not need the directional Sin(f,theta).
!      However, to keep the code simple I pass Sin(f,theta) to calc_tau_total
!      routine regardless.

!   Note: "_S" denotes "short" (fewer freqs)
!         "_L" denotes "long" (freqs out to 10 Hz)
!         (1d version: was used for calculations, but now used only
!           for diagnostic output)
!         (2d version: added and now used for calculations)

    REAL             , INTENT(OUT) ::  LFACTOR_L(:)
    REAL             , INTENT(OUT) ::  ZE
    REAL             , INTENT(IN)  ::  S_IN(:,:)
    REAL             , INTENT(IN)  ::  SPCDIR(:,:)
    REAL             , INTENT(IN)  ::  UFRIC
    REAL             , INTENT(IN)  ::  DDIR_RAD
    REAL             , INTENT(IN)  ::  FRINTF
    REAL             , INTENT(IN)  ::  PWIND(:)
    REAL             , INTENT(IN)  ::  CINV_S(:)
    REAL             , INTENT(IN)  ::  SIGMA_S(:)
    REAL             , INTENT(IN)  ::  GRAV
    REAL             , INTENT(IN)  ::  TAUX_linear,TAUY_linear,WIND10,CTHETA_WIND,STHETA_WIND
    LOGICAL          , INTENT(IN)  ::  TESTFL
    LOGICAL          , INTENT(IN)  ::  VECTOR_TAU,TRUE_U10

    REAL             , ALLOCATABLE  ::  S_IN1D_L(:)
    REAL             , ALLOCATABLE  ::  S_IN_L(:,:)
    REAL             , ALLOCATABLE  ::  DF(:)
    REAL             , ALLOCATABLE  ::  CINV_L(:)
    REAL             , ALLOCATABLE  ::  SIGMA_L(:)

    REAL :: DSIGMA
    REAL :: REDUC
    REAL :: SIGN_NEW,SIGN_OLD
    REAL :: FRAT
    REAL :: ERR
    REAL :: RCHANGE
    REAL :: Z0
    REAL :: AFA

!   tau_total_max is the maximum allowed value for total wind stress
!   tau_total is the value for the total wind stress computed from tau_wave and
!     tau_visc (i.e. tau_normal and tau_tangential)
!   tau_sin_x,tau_sin_y are tau contribution from wind input: in this routine,
!     it is only used for test output

    REAL :: TAU_TOTAL, TAU_wav
    REAL :: TAU_TOTAL_MAX
    REAL :: TAU_SIN_X,TAU_SIN_Y
    REAL :: DDIR_DEG
    REAL :: RHOA,RHOW,PI
    REAL :: FREQ_TMP
    REAL :: U10PROXY,TAU_VISC,TAU_VISC_X,TAU_VISC_Y

    INTEGER :: NF_OLD,NF_NEW,NDIR
    INTEGER :: DIR
    INTEGER :: ITER
    INTEGER :: SLOW_DOWN
    INTEGER :: ISOL
    INTEGER :: IS,IDIR

!...S_in(theta,freq)
    NDIR=SIZE(S_IN,1)
    NF_OLD=SIZE(S_IN,2)
    PI=2.0*ACOS(0.0)

    AFA = 0.01

!...find nf_new
    FRAT=FRINTF+1.0 ! =freq(nf_old)/freq(nf_old-1)
    FREQ_TMP=SIGMA_S(NF_OLD)/(2.0*PI)
    NF_NEW=-99
    DO IS=(NF_OLD+1),(NF_OLD+100)
       FREQ_TMP=FREQ_TMP*FRAT
       IF(FREQ_TMP > 10.0)THEN ! 10.0 here is f=10 Hz
          NF_NEW=IS
          EXIT
       ENDIF
    ENDDO

!   shortcoming of this method: the maximum frequency changes slightly from
!   one simulation to the next, depending on the frequency bands the user
!   selects. To correct this, we should stop at 10 Hz and change df of last bin
!   df for the last bin would have to use finite differencing, instead of
!   FRINTF (as used below), e.g. DF=FREQ(nf_new)-FREQ(nf_new-1)

!... allocate arrays on nf_new
    ALLOCATE(S_IN_L(NDIR,NF_NEW))
    ALLOCATE(S_IN1D_L(NF_NEW)) ! for diagnostic output only
    ALLOCATE(DF(NF_NEW))
    ALLOCATE(CINV_L(NF_NEW))
    ALLOCATE(SIGMA_L(NF_NEW))

!... extend freqs to nf_new
    SIGMA_L(1:NF_OLD)=SIGMA_S(1:NF_OLD)
    CINV_L(1:NF_OLD)=CINV_S(1:NF_OLD)
    S_IN_L(:,1:NF_OLD)=S_IN(:,1:NF_OLD)

    DO IS=(NF_OLD+1),NF_NEW
       SIGMA_L(IS)=SIGMA_L(IS-1)*FRAT
       CINV_L(IS)=SIGMA_L(IS)*0.102 ! deep water assumption for hf tail is OK
!...   For Sin, we use an f-2 approximation
!...   ...Sin=Beta*E~sigma*(U/C)^2*E = f * f^2 * f^-5 = f^-2
       DO IDIR=1,NDIR
          S_IN_L(IDIR,IS)=S_IN_L(IDIR,NF_OLD)*(SIGMA_L(NF_OLD)/SIGMA_L(IS))**2
       ENDDO
    ENDDO

    S_IN_L=S_IN_L*(2.0*PI) ! was in units of m2/(rad-Hz)/rad, so put in
                           ! units m2/Hz/rad

    S_IN1D_L=SUM(S_IN_L,1)*DDIR_RAD ! note that variable is S_in(ID,IS)
!NRL    IF(TESTFL)THEN
!NRL       WRITE(416,*)CHTIME,' % CHTIME'
!NRL       WRITE(416,*)' % (freq , old S_in)'
!NRL       DO IS=1,NF_NEW
!NRL          WRITE(416,211)(SIGMA_L(IS)/(2.0*PI)),S_IN1D_L(IS)
!NRL       END DO
!NRL    ENDIF
    S_IN1D_L=0.0

    RHOA=PWIND(16)
    RHOW=PWIND(17)

    TAU_TOTAL_MAX=(UFRIC**2)*RHOA
    DDIR_DEG=DDIR_RAD*180.0/PI

!   Method: Tsagareli (thesis, 2009); Tsagareli et al. (JPO 2010, eq 3.7);
!     based on fitting to data of Banner and Peirson (JFM 1998).

!   New Aug 30 2013: Previously, tau_visc would go to zero for high wind speeds
!     (U10>22 m/s), which looks strange when plotted. This is addressed now by
!     reformulating such that tau_visc remains constant beyond the wind speed
!     at which tau_visc is maximum, i.e. U10 at which dtau/dU=0. This is 14.667
!     m/s, which also happens to be close to the maximum wind speed in the plot
!     of Banner and Peirson, which Tsagareli used. So we are also holding
!     tau_visc constant for wind speeds for which we have no guidance from
!     the observations.

    U10proxy=min(WIND10,14.66666666)
    tau_visc=(1.408e-3)*(U10proxy**2) - (6.4e-5)*(U10proxy**3)

!   Notes, Aug 15 2012: Stefan noticed that at low wind speeds, |Cv| is
!     slightly larger than |Cd|. This probably is not reasonable. To address
!     this, we say "tau_visc cannot be more than 90% of tau_total_max,
!     regardless". The 90% is semi-arbitrary, selected because it implies that
!     this "fix" has a very minor impact on computations, relative to, say,
!     60%.

    tau_visc=min((0.9*tau_total_max),tau_visc)

!   Note, Aug 30 2013: We may consider using a formulation based on "Law of the
!     Wall", (the subset that deals with  hydraulically smooth flow). This is
!     used by Makin et al. (BLM 1995), Donelan et al. (JGR 2012), Edson et al
!     (JPO 2013),  Hayes et al. (Icarus 2013), and George Mellor, among
!     others. The method results in tau_visc very similar to that of our
!     existing method in cases of low winds (U10<3 m/s) which is the regime
!     where tau_visc is expected to be important. However, it is very different
!     at high wind speeds, since it grows indefinitely with U10. The result is
!     that skin drag (tau_visc) is greater than 28% of total stress when U10
!     is greater than 25 m/s. See my matlab script Makin_iter.m. If we do this,
!     we may want to reduce tau_visc to account for sheltering. Donelan et al.
!     (JPO 2012) do this, though they use Cd, and I think using mss would be
!     better. But even if we reduce tau_visc to 15% of total stress in high
!     winds, this still implies a 15% reduction of wind stress, which would be
!     similar to a 8% reduction in Ustar for our hurricane modeling vs. our
!     current approach based on Banner and Peirson.
!     Example computations for U10=40 m/s:
!         tau_visc as implemented prior to Aug 30 2013 : zero
!         tau_visc as implemented after Aug 30 2013: 0.1 Pa
!         tau_visc from Makin_iter.m : 1.4 Pa
!         tau_total: 4.3 Pa

!   We will assume that viscous stress is in wind direction
    TAU_VISC_X=TAU_VISC*CTHETA_WIND
    TAU_VISC_Y=TAU_VISC*STHETA_WIND

!   DF is needed by calc_tau_total
    DO IS=1,NF_NEW
       DSIGMA=SIGMA_L(IS)*FRINTF  ! since DF=FREQ*FRINTF
       DF(IS)=DSIGMA/(2.0*PI)
    END DO

    ! calculate without reduction
    REDUC=0.0
    CALL CALC_TAU_TOTAL(TAU_TOTAL,LFACTOR_L,REDUC,S_IN_L,DF,CINV_L,WIND10,GRAV,&
      RHOW,SPCDIR,DDIR_RAD,VECTOR_TAU,TRUE_U10,UFRIC,TAU_VISC_X,TAU_VISC_Y, &
      TAU_SIN_X,TAU_SIN_Y,TAUX_linear,TAUY_linear,SIGMA_L,TESTFL)

206 FORMAT('tau_total_max = ',F9.6,' ; tau_total = ',F10.6,' ; tau_visc = ',&
           F9.6)
!NRL    IF(TESTFL)THEN
!NRL       WRITE(*,206)TAU_TOTAL_MAX,TAU_TOTAL,TAU_VISC
!NRL    ENDIF

    IF(TAU_TOTAL < TAU_TOTAL_MAX .OR. TAU_TOTAL < 1E-10)THEN
       DO IS=1,NF_OLD
          LFACTOR_L(IS)=1.0
       END DO
       DO IS=NF_OLD+1,NF_NEW
          LFACTOR_L(IS)=0.0
       END DO
! Deallocate and return
       DEALLOCATE(S_IN_L)
       DEALLOCATE(S_IN1D_L)
       DEALLOCATE(DF)
       DEALLOCATE(CINV_L)
       DEALLOCATE(SIGMA_L)
       RETURN
    ELSE
!NRL       IF(TESTFL)THEN
!NRL          WRITE(*,*)'REDUCING SIN TO MAKE TAU_TOTAL = TAU_TOTAL_MAX'
!NRL       ENDIF
    ENDIF

!   test output
!   e.g. 132+ for NGOM Katrina case, 502+ for GMEX HRD_buoy Ivan case...
!   note that if Edens is too large due to weak Sds, this will be large as
!   a result
!NRL    IF(TESTFL)THEN
!NRL       IF(TAU_TOTAL > 500.0)THEN
!NRL          WRITE(*,*)'WARNING: tau_total without reduction seems large: ', &
!NRL          TAU_TOTAL
!NRL       ENDIF
!NRL    ENDIF

    ! calculate *with* reduction, start value arbitrary
    REDUC=0.05
    CALL CALC_TAU_TOTAL(TAU_TOTAL,LFACTOR_L,REDUC,S_IN_L,DF,CINV_L,WIND10,GRAV,&
      RHOW,SPCDIR,DDIR_RAD,VECTOR_TAU,TRUE_U10,UFRIC,TAU_VISC_X,TAU_VISC_Y,&
      TAU_SIN_X,TAU_SIN_Y,TAUX_linear,TAUY_linear,SIGMA_L,TESTFL)
    ERR=TAU_TOTAL-TAU_TOTAL_MAX
    SIGN_NEW=SIGN(1.0,ERR)
!NRL    IF(TESTFL)THEN
!NRL       WRITE(*,207)0,REDUC,0.0,TAU_TOTAL,ERR
!NRL    ENDIF
    SLOW_DOWN=0
    RCHANGE=2.0
    ISOL=0

! Notes: Objective: find value of REDUC that yields desired tau_total.
!        Current method is to change REDUC sequentially, and slow down the
!        incrementing of REDUC if there is a sign change in the err value.
!        [incrementing of REDUC=RCHANGE]. It might be useful to experiment with
!        Newton-Raphson as an alternative

    DO ITER=1,50
       IF(TAU_TOTAL > TAU_TOTAL_MAX)THEN ! increase REDUC
          REDUC=REDUC*RCHANGE
       ELSE ! then reduce REDUC
          REDUC=REDUC/RCHANGE
       ENDIF
       CALL CALC_TAU_TOTAL(TAU_TOTAL,LFACTOR_L,REDUC,S_IN_L,DF,CINV_L,WIND10, &
            GRAV,RHOW,SPCDIR,DDIR_RAD,VECTOR_TAU,TRUE_U10,UFRIC,TAU_VISC_X, &
            TAU_VISC_Y,TAU_SIN_X,TAU_SIN_Y,TAUX_linear,TAUY_linear,SIGMA_L,TESTFL)
       ERR=TAU_TOTAL-TAU_TOTAL_MAX
       SIGN_OLD=SIGN_NEW
       SIGN_NEW=SIGN(1.0,ERR)
207    FORMAT('iter = ',I3,' ; REDUC = ',F9.6,' ; RCHANGE = ',F9.6, &
              ' ; tau_total = ', F9.6,' ; err = ',F7.4,' ; tau_sin_x = ',F9.6, &
              ' ; tau_sin_y = ',F9.6)
!NRL       IF(TESTFL)THEN
!NRL          WRITE(*,207)ITER,REDUC,RCHANGE,TAU_TOTAL,ERR,TAU_SIN_X,TAU_SIN_Y
!NRL       ENDIF

!     once "slow_down" is set for a given iteration, it stays in this state
!     for future iterations
       IF(SIGN_NEW /= SIGN_OLD)THEN
          SLOW_DOWN=1
       ENDIF
       IF(SLOW_DOWN==1)THEN
          RCHANGE=0.5*(1.0+RCHANGE)
       ENDIF

! For notes re: alternative method that was experimented with,
!   see rev 622 and prior.

       IF((ABS(ERR)/TAU_TOTAL_MAX) < 5.E-04)THEN
          ISOL=1
          EXIT
       ENDIF
    END DO ! DO ITER=...

!!  BEGIN: force stop in case of failure
    IF(ISOL==0)THEN
       IF((ABS(ERR)/TAU_TOTAL_MAX).GE.0.001)THEN
          WRITE(*,*)'warning: no solution found'
          WRITE(*,*)'  tau_total = ',TAU_TOTAL,' tau_total_max = ',TAU_TOTAL_MAX, &
               ' err = ',ERR,'(abs(err)/tau_total_max)  = ',(ABS(ERR)/TAU_TOTAL_MAX)
       ENDIF
!NRL       IF((ABS(ERR)/TAU_TOTAL_MAX).GE.1.0)THEN
!NRL          OPEN(411,FILE='DEBUG.DAT',FORM='UNFORMATTED')
!NRL          WRITE(411)NF_NEW
!NRL          WRITE(411)LFACTOR_L  ! LFACTOR_L(nf_new)
!NRL          WRITE(411)S_IN_L     ! S_in_L(ndir,nf_new)
!NRL          WRITE(411)DF         ! DF(nf_new)
!NRL          WRITE(411)CINV_L     ! CINV_L(nf_new)
!NRL          WRITE(411)TAU_TOTAL_MAX
!NRL          WRITE(411)TAU_TOTAL
!NRL          WRITE(411)REDUC
!NRL          WRITE(411)WIND10
!NRL          WRITE(411)RHOW
!NRL          CLOSE(411)
!NRL          WRITE(*,*)'stopping, see DEBUG.DAT'
!NRL          STOP
!NRL       END IF
       IF(TAU_TOTAL.NE.TAU_TOTAL)THEN ! check for NaN
          WRITE(*,*)'stopping due to tau_total=NaN'
          STOP
       ENDIF
    END IF
!!  END: force stop in case of failure

!!  BEGIN: test output for NRL purposes
!NRL210 FORMAT(F9.6,' ',E12.6,' % (freq , Lfactor)')
!NRL    IF(TESTFL)THEN
!NRL       DO IS=1,NF_NEW
!NRL          WRITE(413,210)(SIGMA_L(IS)/(2.0*PI)),LFACTOR_L(IS)
!NRL       END DO
!NRL    ENDIF

    S_IN1D_L=SUM(S_IN_L,1)*DDIR_RAD
!NRL211 FORMAT(F9.6,1X,E12.6)
!NRL    IF(TESTFL)THEN
!NRL       WRITE(414,*)CHTIME,' % CHTIME'
!NRL       WRITE(414,*)' % (freq , new S_in)'
!NRL       DO IS=1,NF_NEW
!NRL          WRITE(414,211)(SIGMA_L(IS)/(2.0*PI)),(S_IN1D_L(IS)*LFACTOR_L(IS))
!NRL       END DO
!NRL    ENDIF
!    S_IN1D_L=0.0
!!  END: test output

    TAU_wav=0.0
    DO IS=1,NF_NEW
       TAU_wav=TAU_wav+RHOW*GRAV*S_IN1D_L(IS)*LFACTOR_L(IS)   &
                    *CINV_L(IS)*DF(IS)
    END DO

    Z0 = AFA * UFRIC**2 / GRAV
    ZE = Z0/SQRT(1-TAU_wav/TAU_TOTAL)

! Deallocate and return
    DEALLOCATE(S_IN_L)
    DEALLOCATE(S_IN1D_L)
    DEALLOCATE(DF)
    DEALLOCATE(CINV_L)
    DEALLOCATE(SIGMA_L)

  END SUBROUTINE CALC_LFACTOR

  SUBROUTINE CALC_TAU_TOTAL(TAU_TOTAL,LFACTOR,REDUC,S_IN,DF,CINV,WIND10,GRAV, &
                            RHOW,SPCDIR,DDIR_RAD,VECTOR_TAU,TRUE_U10,UFRIC,TAU_VISC_X,TAU_VISC_Y, &
                            TAU_SIN_X,TAU_SIN_Y,TAUX_l,TAUY_l,SIGMA,TESTFL)

!   Objective : apply exponential reduction...more reduction at higher U/C
!        with factor smoothly intersecting value=1 at U/C=1. It is not applied
!         for U/C<1 since that would produce an increase...so factor 1
!         for U/C<=1

!     On methods used (Aug 14 2012):
!        Old method (reasonable enough if tau_wave, tau_total, and tau_visc are
!        all in the same direction) (does not make much sense if tau_wave is
!        opposite to wind direction, which can be the case if we have negative
!         wind input)
!          * tau_total=(UFRIC**2)*rhoa
!          * tau_wave_maximum=tau_total-tauv;
!          * pass Sin(f,theta) and reduction factor to subroutine
!            calc_tau_normal
!          * subroutine calc_tau_normal integrates Sin with reduction factor
!            and returns tau_normal
!          * repeat until tau_normal is less than tau_wave_maximum
!        New method:
!          * there is no tau_wave_maximum
!          * tau_total_maximum=(UFRIC**2)*rhoa (this is a scalar)
!          * tau_visc is now a vector, in the wind direction
!          * pass Sin(f,theta), reduction factor, x-component of tauv,
!            y-component of tauv to subroutine calc_tau_total
!          * subroutine calc_tau_total integrates Sin with reduction factor,
!            plus tau_visc components and returns tau_total (a scalar)
!          * repeat until tau_total (a scalar) is less than tau_total_maximum
!            (a scalar)

!     On computational efficiency, David Johnson, March 1 2013, writes: "... I
!        profiled your SWAN code and found out that fully 50% of the total run
!        time is spent in subroutine calc_tau_total, with an effective doubling
!        of runtimes with the new physics - much of this is probably the
!        repeated subroutine calls and memalloc etc.. By inlining the shear
!        stress calculation within calc_Lfactor, I managed to reduce run time
!        by about 40%. " <end quote> I have experimented with this and do not
!        find any improvement to run time (it actually increases run time by
!        1%).Since calc_tau_total is called in 3 separate places, inlining
!        makes the code itself less concise, less readable and more prone to
!        bugs. Therefore, I keep calc_tau_total as a separate subroutine for
!        now.

!     On run time:
!        I find that for a square domain, deep water, run in serial mode,
!        there is the following increase in computation time (expressed as a
!        factor between total run time with new physics and total run time
!        with old physics):
!            domain size 7x7: factor 1.74
!            domain size 101x101: factor 1.27 (with substantial i/o, this drops
!            to 1.21)

    USE SWCOMM3, ONLY: WNDSCL

    IMPLICIT NONE
!
!
!     SWAN (Simulating WAves Nearshore); a third generation wave model
!     Copyright (C) 1993-2017  Delft University of Technology
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation; either version 2 of
!     the License, or (at your option) any later version.
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     A copy of the GNU General Public License is available at
!     http://www.gnu.org/copyleft/gpl.html#SEC3
!     or by writing to the Free Software Foundation, Inc.,
!     59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
!
!
    REAL             , INTENT(OUT) ::  TAU_TOTAL
    REAL             , INTENT(OUT) ::  LFACTOR(:)
    REAL             , INTENT(IN)  ::  REDUC
    REAL             , INTENT(IN)  ::  S_IN(:,:)
    REAL             , INTENT(IN)  ::  SIGMA(:)
    REAL             , INTENT(IN)  ::  DF(:)
    REAL             , INTENT(IN)  ::  CINV(:)
    REAL             , INTENT(IN)  ::  WIND10,UFRIC
    REAL             , INTENT(IN)  ::  GRAV
    REAL             , INTENT(IN)  ::  RHOW
    REAL             , INTENT(IN)  ::  DDIR_RAD
    REAL             , INTENT(IN)  ::  TAU_VISC_X,TAU_VISC_Y,TAUX_l,TAUY_l
    REAL             , INTENT(IN)  ::  SPCDIR(:,:)
    REAL             , INTENT(OUT) ::  TAU_SIN_X,TAU_SIN_Y

!   SPCDIR explained:
!      SPCDIR: (*,1); spectral directions (radians)
!      ECOS  : =SPCDIR(*,2); cosine of spectral directions
!      ESIN  : =SPCDIR(*,3); sine of spectral directions

    LOGICAL          , INTENT(IN)  ::  VECTOR_TAU, TRUE_U10, TESTFL

    REAL             , ALLOCATABLE  :: AEXP(:)
    REAL             , ALLOCATABLE  :: UOVERC(:)
    REAL             , ALLOCATABLE  :: S_IN_RED(:,:)

    INTEGER :: NF,NDIR,IS,IDIR
    LOGICAL :: IS_THIS_A_NAN
    REAL    :: PI,CTHETA,STHETA,TAU_TOTAL_X,TAU_TOTAL_Y,TAU_SIN_SCALAR,TEMP2,TAU_TMP

    NF=SIZE(DF)
    NDIR=SIZE(S_IN,1)
    PI=2.0*ACOS(0.0)

    ALLOCATE(UOVERC(NF))
    ALLOCATE(AEXP(NF))
    ALLOCATE(S_IN_RED(NDIR,NF))

    IF(TRUE_U10)THEN
       TEMP2 = WIND10
    ELSE
       TEMP2 = WNDSCL * UFRIC
    ENDIF
    UOVERC=TEMP2 * CINV ! matrix op

    AEXP=EXP((1-UOVERC) * REDUC) ! MATRIX OP

    DO IS=1,NF
       LFACTOR(IS)=MIN(1.0,AEXP(IS))
       DO IDIR=1,NDIR
          S_IN_RED(IDIR,IS)=S_IN(IDIR,IS) * LFACTOR(IS)
       ENDDO
    ENDDO

!   see notes above about this calculation
    TAU_SIN_X=0.0
    TAU_SIN_Y=0.0
    TAU_SIN_SCALAR=0.0
    DO IS=1,NF
       DO IDIR=1,NDIR
          CTHETA=SPCDIR(IDIR,2)
          STHETA=SPCDIR(IDIR,3)
          TAU_SIN_SCALAR=TAU_SIN_SCALAR+RHOW*GRAV*S_IN_RED(IDIR,IS)*CINV(IS) &
                         *DF(IS)*DDIR_RAD
          TAU_SIN_X=TAU_SIN_X+RHOW*GRAV*S_IN_RED(IDIR,IS)*CINV(IS)*DF(IS) &
                         *CTHETA*DDIR_RAD
          TAU_SIN_Y=TAU_SIN_Y+RHOW*GRAV*S_IN_RED(IDIR,IS)*CINV(IS)*DF(IS) &
                         *STHETA*DDIR_RAD
       END DO
!NRL!      IF(TESTFL)THEN
!NRL!         output stress in cumulative-with-freq format
!NRL!         this generates a large file, so use only when needed
!NRL!         TAU_TMP=SQRT(TAU_SIN_X**2+TAU_SIN_Y**2)
!NRL!         WRITE(417,*)IS,(SIGMA(IS)/(2.0*PI)),TAU_TMP
!NRL!      END IF
    ENDDO

    IF(VECTOR_TAU)THEN
       TAU_TOTAL_X=TAU_SIN_X+TAU_VISC_X+TAUX_l
       TAU_TOTAL_Y=TAU_SIN_Y+TAU_VISC_Y+TAUY_l
       TAU_TOTAL=SQRT(TAU_TOTAL_X**2+TAU_TOTAL_Y**2)
    ELSE ! recover old method, based on scalar addition
       TAU_TOTAL=TAU_SIN_SCALAR+SQRT(TAU_VISC_X**2+TAU_VISC_Y**2)
    ENDIF

!   NaN/Inf check
    IS_THIS_A_NAN = .NOT. ( TAU_TOTAL .GE. -HUGE(TAU_TOTAL) .AND. TAU_TOTAL  &
                    .LE. HUGE(TAU_TOTAL) )
    IF(IS_THIS_A_NAN)THEN
       WRITE(*,*)'UFRIC = ',UFRIC
       WRITE(*,*)'stopping ; tau_total = ',TAU_TOTAL
       STOP
    ENDIF

! Deallocate and return
    DEALLOCATE(UOVERC)
    DEALLOCATE(AEXP)
    DEALLOCATE(S_IN_RED)

  END SUBROUTINE CALC_TAU_TOTAL

!****************************************************************
!
  SUBROUTINE SWIND0_NRL (SPCSIG,THETAW,ANYWND,         &
                         UFRIC,FPM,MEMSINA,SPCDIR,KWAVE)
!
!****************************************************************
!
      USE SWCOMM3
      USE SWCOMM4
      USE OCPCOMM4

      IMPLICIT NONE
!
!
!   --|-----------------------------------------------------------|--
!     | Delft University of Technology                            |
!     | Faculty of Civil Engineering                              |
!     | Environmental Fluid Mechanics Section                     |
!     | P.O. Box 5048, 2600 GA  Delft, The Netherlands            |
!     |                                                           |
!     | Programmers: The SWAN team                                |
!   --|-----------------------------------------------------------|--
!
!
!     SWAN (Simulating WAves Nearshore); a third generation wave model
!     Copyright (C) 1993-2017  Delft University of Technology
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation; either version 2 of
!     the License, or (at your option) any later version.
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     A copy of the GNU General Public License is available at
!     http://www.gnu.org/copyleft/gpl.html#SEC3
!     or by writing to the Free Software Foundation, Inc.,
!     59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
!
!
!     For background info on SWIND0_NRL, see SWIND0

!     Q: How is SWIND0_NRL different from SWIND0 ?
!     A: linear wind input term is reduced in some cases
!
!     See inline comments re:  methods of reducing the linear wind input term
!     Method 1: introduce a dependence on freq after f=1 Hz to draw it down
!     Method 2: limit the overall contribution to stress to 1% of the stress
!       that one would calculate from Ustar
!     Note that it is OK to use both together.
!     However, at high wind speeds (e.g. 38 m/s), Method 1 by itself is not
!       enough, since the contribution can grow to 37% of the total stress.
!     Note, Method 1 by itself actually causes a small *increase* in SWH for my
!       U=12 m/s case. I'm not sure why that is. Maybe via SNL.
!     Without either method, overall contribution to stress can be 1400% of
!       the stress calculated from Ufric! (factor 14)

      REAL    SPCDIR(MDC,6)
      REAL    SPCSIG(MSC)
      INTEGER ID      ,IS
      REAL    FPM     ,UFRIC   ,THETA   ,THETAW
      REAL    SWINEA(MDC,MSC)
      REAL    TEMP1   ,TEMP2
      REAL    CTW     ,STW     ,COSDIF
      REAL    TEMP3   ,FILTER
      REAL    MEMSINA(MDC,MSC,MCGRD)
      REAL    SIGMASIGMA1,REDUCTION
      REAL    KWAVE(MSC,MICMAX)
      LOGICAL ANYWND(MDC)
      REAL    TAUX_linear,TAUY_linear,CINV2,CTH,STH,TAU_FROM_UFRIC,TAU_MAGNITUDE
      INTEGER IENT
      REAL    ARGU

      SAVE IENT
      DATA IENT/0/
      IF (LTRACE) CALL STRACE (IENT,'SWIND0_NRL')
!
!     *** calculate linear wind input term ***
!
      CTW = COS(THETAW)
      STW = SIN(THETAW)
      FPM =  GRAV / ( 28.0 * UFRIC )
      TEMP1 = PWIND(31) / ( GRAV**2 * 2. * PI )
      SWINEA=HUGE(SWINEA) ! initialize for error catching

      DO IS = 1, MSC
        ARGU   = MIN (2., FPM / SPCSIG(IS))
        FILTER = EXP ( - ARGU**4 )

! note that this SIGMA is not in eq A-2 of Ris (1997). Thus, we are calculating
! "A/sigma" here, not "A"

        TEMP2  = TEMP1 / SPCSIG(IS)
        DO ID = 1, MDC
! Note that we include the entire range for ID, since we need to calc stress

           IF ( ANYWND(ID) .AND. SPCSIG(IS) .GE. (0.7 * FPM) ) THEN
              THETA  = SPCDIR(ID,1)
              COSDIF = SPCDIR(ID,2)*CTW + SPCDIR(ID,3)*STW
              TEMP3  = ( UFRIC *  MAX( 0. , COSDIF))**4
              SWINEA(ID,IS) = MAX( 0. , TEMP2 * TEMP3 * FILTER )

! begin method 1 for reduction
              REDUCTION=1.0
              SIGMASIGMA1=SPCSIG(IS)/(2.0*PI)
! SIGMASIGMA1 is SIGMA/SIGMA1, where SIGMA1 corresponds to f=1 Hz.
              IF(SIGMASIGMA1.GT.1.0)THEN
                 REDUCTION=SIGMASIGMA1**(-3)
              ENDIF
              SWINEA(ID,IS)=SWINEA(ID,IS)*REDUCTION
! end method 1 for reduction
           ELSE
              SWINEA(ID,IS) = 0.0
           END IF
        ENDDO
     ENDDO

! calculate total stress
! recall:   PWIND(16) 1D    Rho air (1.28)
      TAU_FROM_UFRIC=PWIND(16)*UFRIC**2

! calculate stress from linear wind input
      TAUX_linear=0.0
      TAUY_linear=0.0
      DO IS = 1, MSC
         CINV2=KWAVE(IS,1)/SPCSIG(IS)
         DO ID = 1, MDC
! Note that we include the entire range for ID, since we need to calc stress
            CTH = SPCDIR(ID,2) ! cos(theta)
            STH = SPCDIR(ID,3) ! sin(theta)
            TAUX_linear=TAUX_linear+CTH*CINV2*SWINEA(ID,IS)*SPCSIG(IS)*DDIR*FRINTF*SPCSIG(IS)
            TAUY_linear=TAUY_linear+STH*CINV2*SWINEA(ID,IS)*SPCSIG(IS)*DDIR*FRINTF*SPCSIG(IS)
         END DO
      END DO
! recall:    PWIND(17)   Rho water (1025)
      TAU_MAGNITUDE=SQRT(TAUX_linear**2+TAUY_linear**2)*PWIND(17)*GRAV
      IF (TAU_MAGNITUDE.GT.(0.01*TAU_FROM_UFRIC))THEN
         REDUCTION=(0.01*TAU_FROM_UFRIC)/TAU_MAGNITUDE ! method 2 for reduction
         SWINEA=SWINEA*REDUCTION
!NRL         IF(TESTFL)THEN
!NRL            WRITE(*,*)'swind0 : tau_magnitude,tau_from_ufric,reduction 2 = ' ,&
!NRL            TAU_MAGNITUDE,TAU_FROM_UFRIC,REDUCTION
!NRL         ENDIF
      ENDIF

      DO IS = 1, MSC
         IF ( SPCSIG(IS) .GE. (0.7 * FPM) ) THEN
            DO ID = 1, MDC
               MEMSINA(ID,IS,KCGRD(1)) = SWINEA(ID,IS)
!              *** test output ***
               IF (ITEST .GE. 80 .AND. TESTFL )  &
               WRITE (PRTEST, 333) ID, IS, FILTER, SWINEA(ID,IS)
 333           FORMAT (' ID IS FILTER  WIND SOURCE ',2I4, 1X, 2(1X,E11.4))
            ENDDO
         ENDIF
      ENDDO
!
! calculate stress (test point only)
!NRL      IF(TESTFL)THEN
!NRL         TAUX_linear=0.0
!NRL         TAUY_linear=0.0
!NRL         DO  IS = 1, MSC
!NRL            CINV2=KWAVE(IS,1)/SPCSIG(IS)
!NRL            DO ID = 1, MDC
!NRL! Note that we include the entire range for ID, since we need to calc stress
!NRL               CTH = SPCDIR(ID,2) ! cos(theta)
!NRL               STH = SPCDIR(ID,3) ! sin(theta)
!NRL! SPCSIG(IS) replaces "EN" as used in SWIND_Donelan: Thus, factor AC2 is
!NRL! removed, EN=sig*ac2
!NRL               TAUX_linear   =TAUX_linear +CTH*CINV2*SWINEA(ID,IS)*SPCSIG(IS) &
!NRL                              *DDIR*FRINTF*SPCSIG(IS)
!NRL               TAUY_linear   =TAUY_linear +STH*CINV2*SWINEA(ID,IS)*SPCSIG(IS) &
!NRL                              *DDIR*FRINTF*SPCSIG(IS)
!NRL            END DO
!NRL         END DO
!NRL! recall:    PWIND(17) 1D    Rho water (1025)
!NRL         TAUX_linear=TAUX_linear*PWIND(17)*GRAV
!NRL         TAUY_linear=TAUY_linear*PWIND(17)*GRAV
!NRL         WRITE(*,*)'SWIND0: TAUX,TAUY = ',TAUX_linear,TAUY_linear
!NRL      ENDIF
!
!     *** test output ***
!
      IF (ITEST.GE. 60.AND.TESTFL) THEN
        WRITE(PRINTF,400) KCGRD(1), THETAW*180./PI
 400    FORMAT(' SWIND0: POINT  THETAW       :',I5,E12.4)
        WRITE(PRINTF,500) TEMP1, FPM, UFRIC
 500    FORMAT(' SWIND0: TEMP1 FPM UFRC     :',3E12.4)
        WRITE(PRINTF,*)
        IF (ITEST.GE. 120.AND.TESTFL) THEN
          DO IS = 1, MSC
            DO ID = 1, MDC
              WRITE(PRINTF,100) IS,ID,ANYWND(ID)
 100          FORMAT(' IS ID ANYWND : ', 2I5,1X,L2)
            ENDDO
          ENDDO
        ENDIF
      END IF
!
      RETURN
!     end of subroutine SWIND0_NRL
      END SUBROUTINE SWIND0_NRL
!
subroutine filsin ( memsin, idcmin, idcmax, imatra, anywnd, plwnds, isstop, genc0 )
!
!   --|-----------------------------------------------------------|--
!     | Delft University of Technology                            |
!     | Faculty of Civil Engineering and Geosciences              |
!     | Environmental Fluid Mechanics Section                     |
!     | P.O. Box 5048, 2600 GA  Delft, The Netherlands            |
!     |                                                           |
!     | Programmer: Marcel Zijlema                                |
!   --|-----------------------------------------------------------|--
!
!
!     SWAN (Simulating WAves Nearshore); a third generation wave model
!     Copyright (C) 1993-2017  Delft University of Technology
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation; either version 2 of
!     the License, or (at your option) any later version.
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     A copy of the GNU General Public License is available at
!     http://www.gnu.org/copyleft/gpl.html#SEC3
!     or by writing to the Free Software Foundation, Inc.,
!     59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
!
!
!   Authors
!
!   40.88: Marcel Zijlema
!   40.88: Erick Rogers
!
!   Updates
!
!   40.88, May 2017: New subroutine
!   40.88, Jun 2017: Check for ZIEGER to allow case of negative wind input,
!                    wave direction opposing wind direction
!
!   Purpose
!
!   Fills IMATRA and GENC0 arrays with wind input term for every gridpoint per sweep direction
!
!   Modules used
!
    use swcomm3
    use swcomm4
    use ocpcomm4
!
    implicit none
!
!   Argument variables
!
    integer, intent(in)                         :: isstop ! maximum frequency that is propagated within a sweep
    !
    integer, dimension(MSC), intent(in)         :: idcmax ! maximum frequency-dependent counter in directional space
    integer, dimension(MSC), intent(in)         :: idcmin ! minimum frequency-dependent counter in directional space
    !
    real, dimension(MDC,MSC)      , intent(out) :: imatra ! coefficients of right hand side of action balance equation
    real, dimension(MDC,MSC,MGENR), intent(out) :: genc0  ! explicit part of generation in present vertex for output purposes
    real, dimension(MDC,MSC,MCGRD), intent(in)  :: memsin ! wind input stored in all active grid points
    real, dimension(MDC,MSC,NPTST), intent(out) :: plwnds ! explicit part of wind input for test output
    !
    logical, dimension(MDC)       , intent(in)  :: anywnd ! determine if wind input is active for bin
!
!   Local variables
!
    integer       :: id       ! loop counter over direction bins
    integer       :: iddum    ! counter in directional space for considered sweep
    integer, save :: ient = 0 ! number of entries in this subroutine
    integer       :: is       ! loop counter over frequency bins
!
!   Structure
!
!   Description of the pseudo code
!
!   Source text
!
    if (ltrace) call strace (ient,'filsin')
    !
    do is = 1, isstop
       do iddum = idcmin(is), idcmax(is)
          id = mod ( iddum - 1 + MDC , MDC ) + 1
          if ( anywnd(id).or.ZIEGER ) then
             if(TESTFL) plwnds(id,is,IPTST) = plwnds(id,is,IPTST) +  &
                                              memsin(id,is,KCGRD(1))
             imatra(id,is) = imatra(id,is)  + memsin(id,is,KCGRD(1))
             genc0(id,is,1)= genc0(id,is,1) + memsin(id,is,KCGRD(1))
          endif
       enddo
    enddo
    !
    if ( TESTFL .and. ITEST > 50 ) then
       write (PRINTF,101) idcmin(1), idcmax(1), MSC, isstop
       if ( ITEST > 100 ) then
          do is = 1, isstop
             do iddum = idcmin(is), idcmax(is)
                id = mod ( iddum - 1 + MDC , MDC ) + 1
                write (PRINTF,102) is, id, memsin(id,is,KCGRD(1))
             enddo
          enddo
       endif
    endif
    !
 101 format(' FILSIN: ID_MIN ID_MAX MSC ISTOP :',4i6)
 102 format(' FILSIN: IS ID MEMSIN()          :',2i6,e12.4)
    !
end subroutine filsin
!
  !****************************************************************
  SUBROUTINE SSWELL_ROGERS (SPCSIG  ,KWAVE   ,IDCMIN  ,IDCMAX,ISSTOP , DISSC1 ,ETOT &
                           ,IMATDA,URMSTOP, GRAV   , RHOAW  , MDC,TESTFL,IPTST,PLSWEL,CGO,CDSV,FESWELL)

    ! note that CGo is for diagnostic purposes only
    ! excluded : SPCDIR AC2 DEP2 IMATRA

    ! WARNING: This source term may have a large impact in stationary
    ! computations with SWAN, e.g. on scales of 2 deg x 2 deg

    !****************************************************************

    IMPLICIT NONE
!
!
!     SWAN (Simulating WAves Nearshore); a third generation wave model
!     Copyright (C) 1993-2017  Delft University of Technology
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation; either version 2 of
!     the License, or (at your option) any later version.
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     A copy of the GNU General Public License is available at
!     http://www.gnu.org/copyleft/gpl.html#SEC3
!     or by writing to the Free Software Foundation, Inc.,
!     59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
!
!
    LOGICAL, INTENT(IN) :: TESTFL
    INTEGER, INTENT(IN) :: ISSTOP,MDC,IPTST
    INTEGER, INTENT(IN) :: IDCMIN(:) ! IDCMIN(MSC)
    INTEGER, INTENT(IN) :: IDCMAX(:) ! IDCMAX(MSC)

    REAL   , INTENT(IN)    :: ETOT
    REAL   , INTENT(IN)    :: KWAVE(:,:) ! KWAVE(MSC,MICMAX)
    REAL   , INTENT(IN)    :: CGo(:,:)   ! CGo(MSC,MICMAX)
! CGo: group velocity without currents, but includes depth effects
    REAL   , INTENT(IN)    :: SPCSIG(:)  ! SPCSIG(MSC)
    REAL   , INTENT(IN)    :: URMSTOP
! URMSTOP:  this should be the "rms velocity amplitude" *not* the "rms velocity"

    REAL   , INTENT(IN)    :: RHOAW
    REAL   , INTENT(IN)    :: GRAV
    REAL   , INTENT(IN)    :: CDSV
    REAL   , INTENT(IN)    :: FESWELL

    REAL   , INTENT(INOUT) :: IMATDA(:,:)   ! IMATDA(MDC,MSC)
    REAL   , INTENT(INOUT) :: DISSC1(:,:,:) ! DISSC1(1:MDC,1:MSC,1:MDISP)
    REAL   , INTENT(OUT)   :: PLSWEL(:,:,:) ! PLSWEL(MDC,MSC,NPTST)

    ! parameters:
    REAL, PARAMETER   :: RE_CRIT=2.0E+5 ! Re_c or SWELLF4
    !          REAL, parameter   :: Cdsv=1.2 ! C_dsv or SWELLF5
    REAL, PARAMETER   :: NU_AIR=1.51E-5  ! 1.51E-5 m2/s at 20 deg C
    REAL, PARAMETER   :: NU_WATER=1.0E-6  ! 1.0E-6 m2/s at 20 deg C

    !   local variables:
    INTEGER           :: IS
    INTEGER           :: IDDUM,ID
    REAL              :: AORB,USIGTOP
    REAL              :: RE
    REAL              :: SWDIS(ISSTOP) ! this is S_SWELL / E(f,theta)
    REAL              :: KDS_PHILLIPS

! Q: What should fe be?
!   1)  0.004 < fe < 0.013 according to JPO 2010, near final version and
!     March 2009 version
!   2)  "for a smooth surface, fe is of the order 0.002 to 0.008" (Fabrice
!     proposes 0.0035), GRL 2009
!   3)  "Most of the values of fe fall in the range 0.005 to 0.01" (2008
!     Conf paper)
!   4)  Supplemental information for "Ocean swell evolution from distant
!     storms" : experiments with 3 constant values ( 0 , 0.0035 , 0.0070)
!   5)  In Ardhuin et al. (JPO 2010), it is calculated using a formula (eq. 10)
!   6)  Ardhuin et al. (JPO 2010) say: "Assuming a constant fe in (9),
!     Ardhuin et al. (2009b) found that swell observations are consistent
!     with 0.004 < fe < 0.013.

! Q: What should Re_crit be ?
!   1) my original code (up to and including v54):
!     ==>Re_crit=1e+5
!   2) change to this code in v55 :
!     ==>Re_crit=2e+5
!   3) JPO 2010, near-final version: "Semi-empirical dissipation source
!     functions for ocean waves: Part I, definition, calibration and
!     validation."
!     ==>Re_crit=2e+5
!   4) JPO 2010, March 2009 draft
!     ==>Re_crit=1e+5
!   5) GRL 2009: Observation of swell dissipation across oceans
!     ==>Re_crit=1e+5
!   6) conf paper 2008: Spectral wave dissipation based on observations: a
!     global validation
!     ==>Re_crit=1e+5
!   7) auxiliary material for "Strong decay of steep swells observed across
!     oceans"
!     ==>Re_crit=1e+5
!   8) Nature submission (?), "Ocean swell evolution from distant storms"
!     ==>Re_crit=2e+5
!   9) Supplemental information for "Ocean swell evolution from distant
!     storms"
!     ==>Re_crit=1e+5

!   my notes from phone call w/Fabrice:
! "....10^5 w/bug ....2x10^5 w/out bug.....sqrt orbital displacement aorb"
!
! Method:
!    @N/@t=S/sig=C(sig)*E(sig,theta)/sig=C(sig)*N(sig,theta)
!    where C(sig) has units of rad/sec
!    ...so if we are passing to IMATDA, we just need to pass C(sig)
!    Here, SWDIS(IS) is my C(sig)

! Notes about proportionality:
! Both S_{swell,Babanin} and S_{swell,Ardhuin} are proportional
!      to a^1 and omega^3.
! The only difference is in
!      1) the proportionality coefficient
!      2) the manner in which amplitude is converted to E(f)

! The dissipation by viscosity in the water is different. It is proportional to
!      a^0 and omega^4.

    AORB=2.0*SQRT(ETOT)
    USIGTOP=URMSTOP*SQRT(2.0)
    RE=4.0*USIGTOP*AORB/NU_AIR

    IF(RE > RE_CRIT)THEN
       DO IS=1, ISSTOP
          ! note that "-1" omitted since LHS
          SWDIS(IS) = RHOAW * 16.0 * FESWELL * SPCSIG(IS)**2  &
                      * USIGTOP / GRAV
       END DO
    ELSE
       DO IS=1, ISSTOP
          SWDIS(IS) = RHOAW * CDSV * 2.0 * KWAVE(IS,1)  &
                      * SQRT(2.0 * NU_AIR * SPCSIG(IS))
       END DO
    ENDIF

! Add dissipation by viscosity in the water.
! Method: Phillips pg 37 "the energy density of the wave field decreases
!   as exp(-4*nu*k^2*t)"
!   Thus, Kds_Phillips=4*nuw*(k^2), where Kds is (@E/@t)/E
! For other ideas about S_visc, see Pierson et al. (2013+),
!    "Rain-induced attenuation of water waves"

! Alternate reference: Dulov and Kosnnik (Izvestiya, Atmospheric and
! Oceanic Physics, 2009, vol 45, No 3, pp 380-391) :
!  gamma=-4*nu*k^2, where nu=1.3e-6 m2/s
! They also say, "This theory suggests that local three-wave
! interactions are a dominant mechanism for forming the spectrum in the
! entire capillary range except for its shortest wavelength portion,
! where viscous dissipation becomes important"
! The implication is that viscosity is only important for the shortest
! waves in the capillary range, for "forming the spectrum", and perhaps
! for dissipation also.

! Problem identified: if we use this for calculation of tau_wave_to_atm,
!   then we should not include Sds,visc...that part would go to
!   tau_wave_to_ocean

!   DO IS=1, ISSTOP
!      KDS_PHILLIPS=4.0*NU_WATER*(KWAVE(IS,1)**2)  ! UNITS OF RADIAN^2/S
!      SWDIS(IS) = SWDIS(IS) + KDS_PHILLIPS
!   END DO

    DO IS=1, ISSTOP
       !         Only fill the values for the current sweep
       DO IDDUM = IDCMIN(IS), IDCMAX(IS)
          ID = MOD(IDDUM - 1 + MDC, MDC) + 1
          IMATDA(ID,IS)   = IMATDA(ID,IS)   + SWDIS(IS)
          DISSC1(ID,IS,4) = DISSC1(ID,IS,4) + SWDIS(IS)
          IF (TESTFL) PLSWEL(ID,IS,IPTST) = -1.*(SWDIS(IS))
       END DO
    END DO

    RETURN
  END SUBROUTINE SSWELL_ROGERS

  SUBROUTINE SSWELL_ARDHUIN (SPCSIG, THETAW, KWAVE, IDCMIN, IDCMAX, ISSTOP, DISSC1, ETOT, IMATDA &
                            ,SPCDIR, UFRIC, URMSTOP, GRAV   , RHOAW  , TESTFL,IPTST,PLSWEL,MDC,CGO,CDSV)
!---------------------------------------------------------------------------------------------
!  SUBROUTINE SSWELL_ARDHUIN calculates swell dissipation using equations (8) and (9) in
!  Ardhuin et al (2010). These equations are the same as Rogers et al 2012 in Subroutine
!  SSWELL_ROGERS, except fe is calculated using wind-wave angle, friction velocity and
!  the significant surface orbital velocity.
!
!
!     SWAN (Simulating WAves Nearshore); a third generation wave model
!     Copyright (C) 1993-2017  Delft University of Technology
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation; either version 2 of
!     the License, or (at your option) any later version.
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     A copy of the GNU General Public License is available at
!     http://www.gnu.org/copyleft/gpl.html#SEC3
!     or by writing to the Free Software Foundation, Inc.,
!     59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
!
!
!  0. Authors
!     Yalin Fan
!
!  1. Updates
!
!  2. PURPOSE
!     Compute swell dissipation / attenuation for a third generation wave model
!
!  3. METHOD
!
!     Re > Recr
!     Sswell(f,theta) = 16 * (Rho_a/Rho_w) * fe * sigma^2 *uorb * E(f,theta) / g  (9)
!           fe = s1 * {feGM + [|s3| + s2 * cos(theta - theta_wnd)] * u* /uorb }
!
!     Ardhuin calculates using method by Grant and Madsen, which is quite involved,
!     so a shortcut taken here, based on what Ardhuin says: "fe,GM [is] of the order
!     of 0.003" (pg 1921, right column)
!
!     Re <= Recr
!     Sswell(f,theta) = 2 * (Rho_a/Rho_w) * Cdsv * k * sqrt(2*nv*sigma) * F(f,theta)  (8)
!
!  4. Argument variables
!
!     SPCSIG: Relative frequencies in computational domain in sigma-space
!     SPCDIR: (*,1); spectral directions (radians)                        30.82
!             (*,2); cosine of spectral directions                        30.82
!             (*,3); sine of spectral directions                          30.82
!             (*,4); cosine^2 of spectral directions                      30.82
!             (*,5); cosine*sine of spectral directions                   30.82
!             (*,6); sine^2 of spectral directions                        30.82
!        IS          Counter of relative frequency band
!        ID          Counter of directional distribution
!        MSC         Maximum counter of relative frequency
!        MDC         Maximum counter of directional distribution
!
!        one and more dimensional arrays:
!        ---------------------------------
!        KWAVE     2D    Wavenumber
!        IMATDA    2D    Coefficients of main diagonal of matrix
!        IDCMIN    1D    Frequency dependent counter
!        IDCMAX    1D    Frequency dependent counter
!
!     5. SUBROUTINES CALLING
!        SOURCE
!-------------------------------------------------------------------------------------

    IMPLICIT NONE
!
    LOGICAL, INTENT(IN) :: TESTFL
    INTEGER, INTENT(IN) :: ISSTOP,MDC,IPTST
    INTEGER, INTENT(IN) :: IDCMIN(:) ! IDCMIN(MSC)
    INTEGER, INTENT(IN) :: IDCMAX(:) ! IDCMAX(MSC)

    REAL   , INTENT(IN)    :: ETOT
    REAL   , INTENT(IN)    :: KWAVE(:,:)  ! KWAVE(MSC,MICMAX)
    REAL   , INTENT(IN)    :: CGo(:,:)    ! CGo(MSC,MICMAX)
! CGo: group velocity without currents, but includes depth effects
    REAL   , INTENT(IN)    :: SPCSIG(:)   ! SPCSIG(MSC)
    REAL   , INTENT(IN)    :: SPCDIR(:,:) ! SPCDIR(MDC,6)
    REAL   , INTENT(IN)    :: URMSTOP
! URMSTOP:  this should be the "rms velocity amplitude" *not* the "rms velocity"

    REAL   , INTENT(IN)    :: RHOAW
    REAL   , INTENT(IN)    :: GRAV
    REAL   , INTENT(IN)    :: CDSV
    REAL   , INTENT(IN)    :: THETAW
    REAL   , INTENT(IN)    :: UFRIC

    REAL   , INTENT(INOUT) :: IMATDA(:,:)   ! IMATDA(MDC,MSC)
    REAL   , INTENT(INOUT) :: DISSC1(:,:,:) ! DISSC1(1:MDC,1:MSC,1:MDISP)
    REAL   , INTENT(OUT)   :: PLSWEL(:,:,:) ! PLSWEL(MDC,MSC,NPTST)

    ! parameters:
    REAL, PARAMETER   :: RE_CRIT=2.0E+5 ! Re_c or SWELLF4
    !          REAL, parameter   :: Cdsv=1.2 ! C_dsv or SWELLF5
    REAL, PARAMETER   :: NU_AIR=1.51E-5  ! 1.51E-5 m2/s at 20 deg C
    REAL, PARAMETER   :: NU_WATER=1.0E-6  ! 1.0E-6 m2/s at 20 deg C

    !   local variables:
    INTEGER           :: IS
    INTEGER           :: IDDUM,ID
    REAL              :: AORB,USIGTOP
    REAL              :: RE
    REAL              :: SWDIS(MDC,ISSTOP) ! this is S_SWELL / E(f,theta)
    REAL              :: KDS_PHILLIPS
    REAL              :: FESWELL, CTW, STW, COSDIF
    REAL, PARAMETER   :: s1 = 0.8     ! Ardhuin, table A1
    REAL, PARAMETER   :: s2 = -0.018  ! fixed by Ardhuin (pg 1921, right column)
    REAL, PARAMETER   :: s3 = 0.015   ! fixed by Ardhuin (pg 1921, right column)
    REAL, PARAMETER   :: feGM = 0.003 ! fe Grant and Madsen, see notes above

    CTW = COS(THETAW)
    STW = SIN(THETAW)

    AORB=2.0*SQRT(ETOT)
    USIGTOP=URMSTOP*SQRT(2.0)
    RE=4.0*USIGTOP*AORB/NU_AIR

    IF(RE > RE_CRIT)THEN
       DO IS=1, ISSTOP
          DO IDDUM = IDCMIN(IS), IDCMAX(IS)
             ID = MOD(IDDUM - 1 + MDC, MDC) + 1
             COSDIF = SPCDIR(ID,2)*CTW + SPCDIR(ID,3)*STW
             FESWELL = s1*(feGM + (abs(s3) + s2*COSDIF)*UFRIC/USIGTOP)
             ! note that "-1" omitted since LHS
             SWDIS(ID,IS) = RHOAW * 16.0 * FESWELL * SPCSIG(IS)**2  &
                           * USIGTOP / GRAV
          END DO
       END DO
    ELSE
       DO IS=1, ISSTOP
          DO IDDUM = IDCMIN(IS), IDCMAX(IS)
             ID = MOD(IDDUM - 1 + MDC, MDC) + 1
             SWDIS(ID,IS) = RHOAW * CDSV * 2.0 * KWAVE(IS,1)  &
                            * SQRT(2.0 * NU_AIR * SPCSIG(IS))
          END DO
       END DO
    ENDIF

    DO IS=1, ISSTOP
       !         Only fill the values for the current sweep
       DO IDDUM = IDCMIN(IS), IDCMAX(IS)
          ID = MOD(IDDUM - 1 + MDC, MDC) + 1
          IMATDA(ID,IS)   = IMATDA(ID,IS)   + SWDIS(ID,IS)
          DISSC1(ID,IS,4) = DISSC1(ID,IS,4) + SWDIS(ID,IS)
          IF (TESTFL) PLSWEL(ID,IS,IPTST) = -1.*(SWDIS(ID,IS))
       END DO
    END DO

    RETURN
  END SUBROUTINE SSWELL_ARDHUIN

  SUBROUTINE SSWELL_ZIEGER (SPCSIG, KWAVE, AC2, CG, ISSTOP ,IDCMIN  ,IDCMAX, MDC, DISSC1, IMATDA, TESTFL, IPTST, PLSWEL)

!---------------------------------------------------------------------------------------------
!  SUBROUTINE SSWELL_ZIEGER calculates swell dissipation using equations (21) - (23) in
!  Zieger et al (2013), which are originally from Bowden (1950) and Babanin (2011)
!
!
!     SWAN (Simulating WAves Nearshore); a third generation wave model
!     Copyright (C) 1993-2017  Delft University of Technology
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation; either version 2 of
!     the License, or (at your option) any later version.
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     A copy of the GNU General Public License is available at
!     http://www.gnu.org/copyleft/gpl.html#SEC3
!     or by writing to the Free Software Foundation, Inc.,
!     59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
!
!
!  0. Authors
!     Yalin Fan
!
!  1. Updates
!
!  2. PURPOSE
!     Compute swell dissipation / attenuation for a third generation wave model
!
!  3. METHOD
!
!     Sswell(sigma,d) = -2.0*b1*sigma*sqrt(Bn)*E(sigma,d)
!
!  4. Argument variables
!
! i   SPCSIG: Relative frequencies in computational domain in sigma-space
!
!        IS          Counter of relative frequency band
!        ID          Counter of directional distribution
!        MSC         Maximum counter of relative frequency
!        MDC         Maximum counter of directional distribution
!
!        one and more dimensional arrays:
!        ---------------------------------
!        KWAVE     2D    Wavenumber
!        IMATDA    2D    Coefficients of main diagonal of matrix
!        IDCMIN    1D    Frequency dependent counter
!        IDCMAX    1D    Frequency dependent counter
!
!     5. SUBROUTINES CALLING
!        SOURCE
!-------------------------------------------------------------------------------------

    USE SWCOMM3, ONLY : KCGRD, DDIR, MSC, B1Z

    IMPLICIT NONE
!
    ! Subroutine arguments:
    LOGICAL, INTENT(IN)   :: TESTFL
    INTEGER, INTENT(IN)   :: ISSTOP,MDC,IPTST
    INTEGER, INTENT(IN)   :: IDCMIN(:) ! IDCMIN(MSC)
    INTEGER, INTENT(IN)   :: IDCMAX(:) ! IDCMAX(MSC)
    REAL  , INTENT(IN)    :: SPCSIG(:)  ! SPCSIG(MSC)
    REAL  , INTENT(IN)    :: CG(:,:)    ! CG(MSC,ICMAX)
    REAL  , INTENT(IN)    :: KWAVE(:,:) ! KWAVE(MSC,MICMAX)
    REAL  , INTENT(IN)    :: AC2(:,:,:) ! AC2(MDC,MSC,MCGRD)
    REAL  , INTENT(INOUT) :: IMATDA(:,:) ! IMATDA(MDC,MSC)
    REAL  , INTENT(INOUT) :: DISSC1(:,:,:) ! DISSC1(1:MDC,1:MSC,1:MDISP)
    REAL  , INTENT(OUT)   :: PLSWEL(:,:,:) ! PLSWEL(MDC,MSC,NPTST)

    ! Local variables:
    INTEGER              :: IDDUM ,ID    ,IS
    REAL                 :: DMAX,AINV,BN
    ! Ktheta is like D(theta), except max value at each freq is unity
    REAL                 :: KTHETA(MSC,MDC)
    REAL                 :: ANAR(MSC),SIGDENS(MSC),SQRTBN(MSC),CINV(MSC)
    REAL                 :: SWDIS(ISSTOP) ! this is S_SWELL / E(f,theta)

!   Calculate 1d spectrum E(sigma)
    DO  IS = 1, MSC
       SIGDENS(IS) = 0.
       DO  ID = 1, MDC
          SIGDENS(IS) = SIGDENS(IS) + SPCSIG(IS) * AC2(ID,IS,KCGRD(1))
          KTHETA(IS,ID)=AC2(ID,IS,KCGRD(1))
       END DO
       SIGDENS(IS)=SIGDENS(IS)*DDIR  ! units m^2/(radHz)
    END DO

!   Calculate K(sigma,theta)
    DO  IS = 1, MSC
       DMAX=-1.0
       DO  ID = 1, MDC
          IF(KTHETA(IS,ID).GT.DMAX)DMAX=KTHETA(IS,ID)
       END DO
       IF(DMAX.EQ.0.0)THEN
          ! a fix for freq bins (usually first two or so) that are empty
          DMAX=1.0
          DO  ID = 1, MDC
             KTHETA(IS,ID)=1.0
          END DO
       ENDIF

!      Optional: Calculate circular RMS directional spread (for comparison with
!      AINV) :
!      rmsdir(is)=kuik(ktheta(is,:),ddir,spcdir(:,2),spcdir(:,3),.true./.false)

       DO  ID = 1, MDC
          KTHETA(IS,ID)=KTHETA(IS,ID)/DMAX
       END DO
    END DO

!   Calculate SQRT(Bn(IS))
    DO IS=1,MSC
       AINV=0.0
       DO  ID = 1, MDC
          AINV=AINV+KTHETA(IS,ID)*DDIR
       END DO

!   Optional: use RMSDIR in place of AINV: to see how to do this, read notes in
!      rev 622 and prior

       ANAR(IS)=1.0/AINV
       BN=ANAR(IS)*SIGDENS(IS)*(KWAVE(IS,1)**3)*CG(IS,1)
       SQRTBN(IS)=SQRT(BN)
!       SQRTBN(IS)=ANAR(IS)*(SIGDENS(IS)*(KWAVE(IS,1)**3)*CG(IS,1))
    END DO

!   Calculate Swell dissipation Sswl divided by E(sigma,theta) (F(sigma,theta)
!   This is from WW3v4 and Zieger et al. (2015) eq. 23.
!   This does *not* include the steepness-dependent B1 that is introduced in
!   WW3v5 and Zieger et al. (2015) eq. 28.
    DO IS=1, ISSTOP
          SWDIS(IS) = 2.0*B1Z*SPCSIG(IS)*SQRTBN(IS)/3.0
    END DO

!   Re: usage of DISCC1(_,_,4): In case of Zieger Sswell,
!   this should go to tau_wave_to_ocean, not tau_wave_to_atm

    DO IS=1, ISSTOP
       !         Only fill the values for the current sweep
       DO IDDUM = IDCMIN(IS), IDCMAX(IS)
          ID = MOD(IDDUM - 1 + MDC, MDC) + 1
          IMATDA(ID,IS)   = IMATDA(ID,IS)   + SWDIS(IS)
          DISSC1(ID,IS,4) = DISSC1(ID,IS,4) + SWDIS(IS)
          IF (TESTFL) PLSWEL(ID,IS,IPTST) = -1.*(SWDIS(IS))
       END DO
    END DO

    RETURN
  END SUBROUTINE SSWELL_ZIEGER


!----------------------------------------------------------------------------
  SUBROUTINE integrate(ansNum, x,y,np)
!-------------------------------------------------------------------------

    IMPLICIT NONE
!
!
!     SWAN (Simulating WAves Nearshore); a third generation wave model
!     Copyright (C) 1993-2017  Delft University of Technology
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation; either version 2 of
!     the License, or (at your option) any later version.
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     A copy of the GNU General Public License is available at
!     http://www.gnu.org/copyleft/gpl.html#SEC3
!     or by writing to the Free Software Foundation, Inc.,
!     59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
!
    !
    !    Use Trapezoidal Rule to Integrate the area under the curve
    !    specified by the data points in x and y
    !
    !    John Mahaffy 4/9/95
    !
    !
    !   Local Variables
    !
    !    esterr   -   estimated error for Trapizoidal integration
    !    sum2     -   Trapizoidal integration using every other available point
    !

    REAL          , INTENT(IN)  ::  X(:)
    REAL          , INTENT(IN)  ::  Y(:)
    REAL          , INTENT(OUT) :: ANSNUM
    INTEGER       , INTENT(IN)  :: NP
    INTEGER                     :: I

    ANSNUM=0.0
    !  np=size(x) ! we cannot do this because array may be partially filled

    DO  I=1,NP-1
       ANSNUM=ANSNUM+.5*(Y(I)+Y(I+1))*(X(I+1)-X(I))
    END DO

    RETURN
  END  SUBROUTINE INTEGRATE
!-------------------------------------------------------------------------


      SUBROUTINE SURF_ROUGH_FAN (FPI,U,UST,CD)
!
!     Parameter list
!     ----------------------------------------------------------------
!       FPI     Real   I   Peak-input frequency.
!       U       Real   I   True Wind speed.
!       UST     Real   O   Friction velocity.
!       CD      Real   O   Drag coefficient.
!     ----------------------------------------------------------------
!/
      IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
      REAL, INTENT(IN)        :: FPI, U
      REAL, INTENT(OUT)       :: UST, CD
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
      INTEGER                 :: I1, ITT, I
      REAL                    :: SQRTH, SIX, R1, WNP, CP, UNZ, ALPHA, &
                                 RDCH, AFP, Z0,                       &
                                 G,CAPA,PII,CW,Z0FTN,CDFTN,CPEAK,     &
                                 ZCH,X,WAGE,A1,B1,A2,B2,FX,xx1,slope
!/
!/ ------------------------------------------------------------------- /
!/
!
! 0. Definition of constant
      G=9.8
      CAPA=0.4
      PII=3.141517
!
!----      A1=0.57   produce too high Hs  ----------
!---- First try, use 0.4521, still too large
!---- Second try, use 0.3931, Hs is still 1.5m higher
!---- Third try, use 0.354
      A1=0.354
      B1=-1.5
      B2=1.5
!
!---- a transition point of wave age
      CW=11

!
! 1. Defining CD function ---------------------------------------------- *
        IF(U.LE.0.) THEN
          Z0FTN=0.000001
          CDFTN=0.000001
        ELSE IF(U.LE.12.5) THEN
          Z0FTN=0.0185/G*(0.001*U*U+0.028*U)**2
          CDFTN=0.4*0.4*(LOG(10./Z0FTN)**(-2))
        ELSE
          Z0FTN=(0.085*U-0.58)*0.001
          CDFTN=0.4*0.4*(LOG(10./Z0FTN)**(-2))
        END IF
!
!---- UST
      IF(UST.LE.0.) UST=SQRT(CDFTN)*U
!
!---- Phase velocity at peak
      CPEAK=SQRT(G**2/(2.*PII*FPI)**2)
!
!---- First check----------------------------------
      IF(U.LE.0.) THEN
          Z0  =  0.000001
          ZCH =  0.000001
          CD  =  0.000001
          UST =  0.000001
        GO TO 900
      ELSE
!
!----- Initial guess of UST using a function
       X=SQRT(CDFTN)*U
!
!----- Iteration to get UST------------------------------------------- *
!
      DO I=1,5
!
!----- Wave age
         WAGE=CPEAK/X
!
!----- Relationship between wave age and Charnock Coeff.
          xx1=100.0
          slope= 0.15*(U/12.5)
          ZCH=0.023/xx1**slope*WAGE**slope
!
!----- Calculation of Charnock Coefficient(Zch-Wage relationship)

!
         FX=X/CAPA*LOG(10.*G/(X**2*ZCH))-U
!
        IF(ABS(FX).GE. 0.001) THEN
          X=X-(X/CAPA*LOG(10*G/(X**2*ZCH))-U)/                        &
                 ((1/CAPA)*(LOG(10*G/(X**2*ZCH))-2.))
        ELSE
          UST=X
        END IF

       END DO
!------------------ End of Iteration --------------------------------- *
!
        UST=X
        Z0=UST**2*ZCH/G
        IF (Z0.LT.0.00001) Z0 = 0.00001
        CD=(1/CAPA*LOG(10/Z0))**(-2)

!
!----- Second Check--------
!       IF((WAGE.GT.3.).AND.(WAGE.LT.32.)) THEN
!
!------- Roughness length and Drag coefficient
!         Z0=UST**2*ZCH/G
!         CD=(1/CAPA*LOG(10/Z0))**(-2)
!       ELSE
!         Z0=Z0FTN
!         CD=CDFTN
!         UST=SQRT(CD)*U
!         ZCH=Z0*G/UST**2
!       END IF
!---------------------------
!
      END IF
!--------------------------- First Check Ending
!
!----- Final check Cd
!       IF((ABS(CDFTN-CD).GT.0.001).OR.(CD.GT.0.005)) THEN
!          Z0=Z0FTN
!          CD=CDFTN
!          UST=SQRT(CD)*U
!          UST=SQRT(CD)*U/ASF(ISEA)
!          ZCH=Z0*G/UST**2
!       END IF
!
!---- Direction of Stress
!
900   continue

      RETURN

      END SUBROUTINE SURF_ROUGH_FAN


      SUBROUTINE SURF_ROUGH_ECMWF(U, UST, S_IN, SPCSIG, KWAVE, CD)
!------------------------------------------------------------------------------
!     variable list
!     ----------------------------------------------------------------
!       S_IN    Real   I   Wind input source function
!       SPCSIG  Real   I   Relative frequencies in computational domain in sigma-space
!       KWAVE   Real   I   2D    Wavenumber
!       U       Real   I   True Wind speed.
!       UST     Real   O   Friction velocity.
!       CD      Real   O   Drag coefficient.
!     ----------------------------------------------------------------
!/
      USE SWCOMM3
      USE SWCOMM4
      USE OCPCOMM4

      IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
      REAL, INTENT(IN)        :: S_IN(MDC,MSC,MGENR)
      REAL, INTENT(IN)        :: SPCSIG(MSC)
      REAL, INTENT(IN)        :: KWAVE(MSC,ICMAX)
      REAL, INTENT(IN)        :: U
      REAL, INTENT(INOUT)     :: UST
      REAL, INTENT(OUT)       :: CD
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
      INTEGER  :: I1, ITT, I, IS, IDIR, NDIR, NF_OLD, NF_NEW, K, II
      REAL     :: G, CAPA, PII, afa, TAU_wav, RHOW, RHOA, TAU_total
      REAL     :: FREQ_TMP, FRAT, Z0FTN, CDFTN, z0, ze, a, b, X, CDRAG

      REAL, ALLOCATABLE  ::  S_IN_L(:,:,:)
      REAL, ALLOCATABLE  ::  CINV_L(:), SIGMA_L(:), DF_L(:)

!      REAL                    :: SQRTH, SIX, R1, WNP, CP, UNZ, ALPHA, &
!                                 RDCH, AFP, Z0, afa,                  &
!                                 G,CAPA,PII,CW,Z0FTN,CDFTN,CPEAK,     &
!                                 ZCH,X,WAGE,A1,B1,A2,B2,FX,xx1,slope
!/
!/ ------------------------------------------------------------------- /
!/
!
! 0. Definition of constant
      G=9.81
      CAPA=0.4
      PII=3.141517
      AFA = 0.01
      RHOW = 1025.0
      RHOA = 1.28
!
! 1. Difining CD function ---------------------------------------------- *
        IF(U.LE.0.) THEN
          Z0FTN=0.000001
          CDFTN=0.000001
        ELSE IF(U.LE.12.5) THEN
          Z0FTN=0.0185/G*(0.001*U*U+0.028*U)**2
          CDFTN=0.4*0.4*(LOG(10./Z0FTN)**(-2))
        ELSE
          Z0FTN=(0.085*U-0.58)*0.001
          CDFTN=0.4*0.4*(LOG(10./Z0FTN)**(-2))
        END IF
!
!---- UST
      IF(UST.LE.1E-10) THEN
         CDRAG = CDFTN
         UST = SQRT ( CDRAG ) * U
         GOTO 110
      ENDIF

!...S_in(theta,freq)
      NDIR=SIZE(S_IN,1)
      NF_OLD=SIZE(S_IN,2)

!...find nf_new
      FRAT=FRINTF+1.0 ! =freq(nf_old)/freq(nf_old-1)
      FREQ_TMP=SPCSIG(NF_OLD)/(2.0*PII)
      NF_NEW=-99
      DO IS=(NF_OLD+1),(NF_OLD+100)
         FREQ_TMP=FREQ_TMP*FRAT
         IF(FREQ_TMP > 10.0)THEN ! 10.0 here is f=10 Hz
            NF_NEW=IS
            EXIT
         ENDIF
      ENDDO

      ALLOCATE(CINV_L(NF_NEW))
      ALLOCATE(SIGMA_L(NF_NEW))
      ALLOCATE(DF_L(NF_NEW))
      ALLOCATE(S_IN_L(NDIR,NF_NEW,MGENR))

      CINV_L = 0.0
      SIGMA_L = 0.0
      DF_L = 0.0
      S_IN_L = 0.0

      SIGMA_L(1:NF_OLD)=SPCSIG(1:NF_OLD)

      DO IS=1,NF_OLD
         CINV_L(IS)=KWAVE(IS,1)/SPCSIG(IS)
      ENDDO

      DO IS = 1, NF_OLD
         DO IDIR = 1, NDIR
            DO II = 1, MGENR
               S_IN_L(IDIR,IS,II)=S_IN(IDIR,IS,II)
            ENDDO
         ENDDO
      ENDDO

      DO IS=(NF_OLD+1),NF_NEW
         SIGMA_L(IS)=SIGMA_L(IS-1)*FRAT
         CINV_L(IS)=SIGMA_L(IS)*0.102 ! deep water assumption for hf tail is OK
!........For Sin, we use an f^-2 approximation
!.... ...Sin=Beta*E~sigma*(U/C)^2*E = f * f^2 * f^-5 = f^-2
!.... ...since GENC0 (S_IN here) = Sin / f, we use therefore f^-3 as tail
!...   Also, because LFACTOR is applied to the tail part of S_IN before calculate
!...   tau_sin_x and tau_sin_y, and the reduction is significant (on the order of 10^-4)
!...   we should apply the same reduction to GENC0 as well. This requires to bring array
!...   LFACTOR to this subroutine.
!      *** Note that LFACTOR_SRC is renamed here as RDFSIN ***

!Yalin: S_IN(NF_OLD) / GENC0(NF_OLD) was already drawn down by LFACTOR(NF_OLD)
! before pass into routine ADDDIS. But, in CALC_LFACTOR, the tail is added using
! the original S_IN(NF_OLD) before add tail then calculate LFACTOR.
! Thus, we need to restore S_IN(NF_OLD) to its original value before add
! tail here

         DO IDIR=1,NDIR
            DO II = 1, MGENR
               S_IN_L(IDIR,IS,II) = (S_IN_L(IDIR,NF_OLD,II)/RDFSIN(NF_OLD)) &
                                  * (SIGMA_L(NF_OLD)/SIGMA_L(IS))**3 * RDFSIN(IS)
            ENDDO
         ENDDO
      ENDDO

      TAU_wav=0.0
      DO IDIR=1,NDIR
         DO IS=1,NF_NEW
            DF_L(IS)=SIGMA_L(IS)*FRINTF
            TAU_wav=TAU_wav+RHOW*G*S_IN_L(IDIR,IS,1)*CINV_L(IS)*SIGMA_L(IS) &
                         *DF_L(IS)*DDIR
         END DO
      ENDDO

! now we got tau_w, need to cal z0 use tau_w, then iteration
      Z0 = AFA * UST**2 / G
      TAU_total = RHOA * UST**2
      ZE = Z0/SQRT(1-TAU_wav/TAU_total)
      X = SQRT(G*Z0/AFA)*LOG((10.0+ZE-Z0)/Z0)/CAPA - U

      IF (ABS(X) .LT. 1E-6) GO TO 100

      IF (X.GT.0.0) THEN
          A = 0
          B = Z0
      ELSE
          A = Z0
          B = 0.1
      ENDIF

      K = 0
80    Z0 = (A+B)/2
      UST = SQRT(Z0*G/AFA)
      TAU_total = RHOA * UST**2
      ZE = Z0/SQRT(1-TAU_wav/TAU_total)
      X = SQRT(G*Z0/AFA)*LOG((10.0+ZE-Z0)/Z0)/CAPA - U
      IF (X .GT. 0.0) THEN
          B = (A+B)/2
      ELSE
          A = (A+B)/2
      ENDIF
      K = K + 1
      IF (ABS(X) > u*1E-6 .or. X < 0) GOTO 80

100   CONTINUE

      UST = SQRT(Z0*G/AFA)

      DEALLOCATE(CINV_L)
      DEALLOCATE(SIGMA_L)
      DEALLOCATE(DF_L)
      DEALLOCATE(S_IN_L)

110   CONTINUE

      RETURN

      END SUBROUTINE SURF_ROUGH_ECMWF



!------------------------------------------------------------------------------
  REAL FUNCTION KUIK(E_AT_FREQ,DTHETA_RAD,CTH,STH,DIAGNOSTICS)
!------------------------------------------------------------------------------
!
  IMPLICIT NONE
!
!
!     SWAN (Simulating WAves Nearshore); a third generation wave model
!     Copyright (C) 1993-2017  Delft University of Technology
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation; either version 2 of
!     the License, or (at your option) any later version.
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     A copy of the GNU General Public License is available at
!     http://www.gnu.org/copyleft/gpl.html#SEC3
!     or by writing to the Free Software Foundation, Inc.,
!     59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
!
!
!  0. Update history
!
!     Jan 3 2013 Initial version by Erick Rogers
!
!  1. Purpose:
!
!     Compute circular rms directional spreading (Kuik et al.)
!
!  2. Method
!
!     Kiuk et al. (JPO 1988), see also Rogers and Wang (2007) and Mardia book
!     referenced therein.
!
!  3. Parameter list:
!
! Type    I/O        Name    Description
!--------------------------------------------------------------------
  REAL, INTENT(IN) ::  E_AT_FREQ(:) !  E(f,theta) for one frequency f
  REAL, INTENT(IN) ::  DTHETA_RAD  !  delta_theta
  REAL, INTENT(IN) ::  CTH(:)  !  cos(theta)
  REAL, INTENT(IN) ::  STH(:)  !  sin(theta)
  LOGICAL, INTENT(IN) ::  DIAGNOSTICS !  Flag for diagnostic output
!
!  4. Error messages
!
!  5. Called by:
!
!  6. Subroutines used:
!
!  7. Remarks
!
!  8. Structure
!
!  9. Switches
!
! 10. Source code:
!------------------------------------------------------------------------------
! local variables
!
  INTEGER NDIR,IDIR            ! number of directions and counter for directions
  REAL EF                              ! E(f) at this frequency
  REAL INTEG3,INTEG4,a1,b1,m1,DSPR     ! temporary variables for integration
!---------------------------------------------------------------------

  NDIR=SIZE(E_AT_FREQ)

  INTEG3=0.0
  INTEG4=0.0
  EF=0.0
  DO IDIR=1,NDIR
     INTEG3=INTEG3+E_AT_FREQ(IDIR)*CTH(IDIR)*DTHETA_RAD
     INTEG4=INTEG4+E_AT_FREQ(IDIR)*STH(IDIR)*DTHETA_RAD
     EF=EF+E_AT_FREQ(IDIR)*DTHETA_RAD
  END DO
! note that if dtheta_rad is constant, it cancels out, so it can be omitted from this routine
  A1=INTEG3/EF
  B1=INTEG4/EF
  M1=SQRT(A1**2+B1**2)
  DSPR=SQRT(2*(1-M1))

!NRL  IF(DIAGNOSTICS)THEN
!NRL     WRITE(415,*)'% calc_Kuik: using dtheta_rad = ',dtheta_rad,' radians'
!NRL     WRITE(415,*)'% calc_Kuik: DSPR = ',DSPR,' radians'
!NRL  END IF

  KUIK=DSPR

  RETURN
END  FUNCTION KUIK

END MODULE SDSBABANIN