w3ref1md.F90

Go to the documentation of this file.

 
 
!/ ------------------------------------------------------------------- /
MODULE w3ref1md
  !/
  !/                  +-----------------------------------+
  !/                  | WAVEWATCH III                     |
  !/                  |            F. Ardhuin             |
  !/                  |                        FORTRAN 90 |
  !/                  | Last update :         27-Jun-2014 |
  !/                  +-----------------------------------+
  !/
  !/    31-Mar-2010 : Origination.                        ( version 3.14.IFREMER )
  !/    03-Sep-2010 : Clean up                            ( version 3.14.IFREMER )
  !/    31-May-2011 : Adding variable reflections         ( version 4.01 )
  !/    02-Nov-2011 : Compatibility with unst. grids      ( version 4.04 )
  !/    24-Fev-2012 : Correction of angle in fluxes       ( version 4.05 )
  !/    27-Jul-2013 : Adding free infragravity waves      ( version 4.11 )
  !/    11-Nov-2013 : Extends IG energy into main band    ( version 4.13 )
  !/    11-Jun-2014 : Put reflection by subgrids back     ( version 5.01 )
  !/    27-Jun-2014 : Modifies subgrid reflection of IG   ( version 5.01 )
  !/
  !  1. Purpose :
  !
  !     This module computes :
  !        - shoreline reflection
  !        - unresolved islands and iceberg reflections
  !
  !  2. Variables and types :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !     ----------------------------------------------------------------
  !
  !  3. Subroutines and functions :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !      W3SREF    Subr. Public   Reflection of waves (shorline, islands...)
  !     ----------------------------------------------------------------
  !
  !  4. Subroutines and functions used :
  !
  !      Name      Type  Module   Description
  !     ----------------------------------------------------------------
  !      STRACE    Subr. W3SERVMD Subroutine tracing.
  !     ----------------------------------------------------------------
  !
  !  5. Remarks :
  !
  !
  !  6. Switches :
  !
  !     !/S  Enable subroutine tracing.
  !
  !  7. Source code :
  !/
  !/ ------------------------------------------------------------------- /
  !/
  !
  PUBLIC
  !/
  !/ Public variables
  !/
  !
  !/
CONTAINS
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3sref(A, CG, WN, EMEAN, FMEAN, DEPTH, CX1, CY1, REFLC, REFLD,     &
       TRNX, TRNY, BERG, DT, IX, IY, JSEA, S)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |            F. Ardhuin             |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         11-Jun-2014 |
    !/                  +-----------------------------------+
    !/
    !/    06-May-2010 : Origination.                          ( version 3.14-Ifremer )
    !/    31-May-2011 : Introduction of amplitude-dependent R ( version 4.05 )
    !/    27-Jul-2013 : Adding free infragravity waves        ( version 4.11 )
    !/    11-Nov-2013 : Expands IG energy frequency range     ( version 4.13 )
    !/    05-Mar-2014 : Fixing bug with ICALC = 1 and IG1     ( version 4.18 )
    !/    11-Jun-2014 : Put reflection by subgrids back       ( version 5.01 )
    !/
    !  1. Purpose :
    !
    !     Computes coastal and iceberg/island reflections and adds free IG energy
    !
    !  2. Method :
    !
    !     Adds the reflected components from 2 types of sources:
    !        shoreline reflection, subgrid obstruction and icebergs
    !
    !     In the case where the IG switch is present, there are two passes:
    !        - ICALC = 1, only the wind sea and swell are reflected (no IG added)
    !        - ICALC = 2, IG energy is added into all frequency bands
    !
    !
    !     When IG energy is put in the entire spectrum ( NINT(IGPARS(4)).EQ.2 )
    !        two passes are done: the first for the reflection of windsea and swell
    !                             the second for the addition of IG and IG reflection alone
    !
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       A         R.A.  I   Action density spectrum (1-D)
    !       CG        R.A.  I   Group velocities.
    !       WN        R.A.  I   Wavenumbers.
    !       DEPTH     Real  I   Mean water depth.
    !       S         R.A.  O   Source term (1-D version).
    !       D         R.A.  O   Diagonal term of derivative (1-D version).
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      STRACE    Subr. W3SERVMD Subroutine tracing.
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3SRCE    Subr. W3SRCEMD Source term integration.
    !      W3EXPO    Subr.   N/A    Point output post-processor.
    !      GXEXPO    Subr.   N/A    GrADS point output post-processor.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !       None.
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/S  Enable subroutine tracing.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants
    USE w3gdatmd, ONLY: nk, nth, nspec, sig, dth, dden,  smctype, &
         refpars, ecos, esin, ec2, mapth, mapwn, flagll, &
         sig2, dsii, iobpd, gtype, ungtype, mapfs, clgtype, rlgtype
    USE w3gdatmd, ONLY : clats, hpfac, hqfac, sx, sy, si
#ifdef W3_PDLIB
    USE yownodepool, ONLY: pdlib_si
    USE w3gdatmd, ONLY: iobp_loc, iobpd_loc, iobpa_loc, iobdp_loc
#endif
#ifdef W3_IG1
    USE w3gdatmd, ONLY : igpars
    USE w3gig1md
    USE w3canomd, ONLY : w3add2ndorder
    USE w3dispmd, ONLY: nar1d, dfac, n1max, ecg1, ewn1, dsie
#endif
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    !/
    !
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    REAL, INTENT(IN)        :: CG(NK), WN(NK), DEPTH, EMEAN, FMEAN
    REAL, INTENT(INOUT)     :: A(NSPEC)
    REAL, INTENT(IN)        :: CX1, CY1, DT
    INTEGER, INTENT(IN)     :: REFLD(6), IX, IY, JSEA
    REAL, INTENT(IN)        :: REFLC(4), TRNX, &
         TRNY, BERG
    REAL, INTENT(OUT)       :: S(NSPEC)
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER         :: ISPECI, ISPEC, IK, ITH, ITH2, ITH3, ITH2X, ITH2Y, &
         NRS, IK1
    INTEGER         :: ISEA, ICALC, IOBPDIP(NTH)
    LOGICAL         :: IGBCOVERWRITE, IGSWELLMAX
    REAL            :: R1, R2, R3, R4, R2X, R2Y, DEPTHIG
    REAL            :: DELA, DELX, DELY, FACX
    REAL            :: FAC1, FAC2, FAC3, FAC4, RAMP0, RAMP, &
         RAMP1, RAMP2, RAMP4, MICHEFAC, SLOPE
    REAL             :: HS, HIG, HIG1, HIG2, EB, SB, EMEANA, FMEAN2,   &
         FMEANA, FREQIG, EFIG, EFIG1, SQRTH, SMEANA
#ifdef W3_IG1
    INTEGER        :: NKIG,NSPECIG,NSPECIGSTART, I1, I2
    REAL           :: ATMP(NSPEC),ATMP2(NSPEC), STMP1(NSPEC),      &
         STMP2(NSPEC), WNB(NK), CGB(NK), SIX, IGFAC1, IGFAC2
#endif
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
    CALL strace (ient, 'W3SREF')
#endif
    !
    ! 0.  Initializations ------------------------------------------------ *
    !
#ifdef W3_IG1
    igbcoverwrite =(mod( nint(igpars(4)),2).EQ.1)
    igswellmax =( nint(igpars(4)).GE.2)
    ! This following line is a quick fix before the bug is understood ....
    IF (gtype.EQ.ungtype) igswellmax =.false.
    igfac1 = 0.25
    igfac2 = 0.25
#endif
    emeana  = 0.
    fmeana  = 0.
    fmean2  = 0.
 
    delx=1.
    dely=1.
    ! set FACx for all grid types
    IF (flagll) THEN
      facx   =  1./(dera * radius)
    ELSE
      facx   =  1.
    END IF
 
    isea = mapfs(iy,ix)
    !!Li  SMCTYPE shares info with RLGTYPE.  JGLi12Oct2020
    IF (gtype.EQ.rlgtype .OR. gtype.EQ.smctype) THEN
      delx=sx*clats(isea)/facx
      dely=sy/facx
    END IF
 
    IF (gtype.EQ.clgtype) THEN
      ! Maybe what follows works also for RLGTYPE ... to be verified
      delx=hpfac(iy,ix)/ facx
      dely=hqfac(iy,ix)/ facx
    END IF
 
    IF (gtype.EQ.ungtype) THEN
      IF (lpdlib) THEN
#ifdef W3_PDLIB
        delx=5.*sqrt(pdlib_si(jsea))*(dera * radius)    ! first approximation ...
        dely=5.*sqrt(pdlib_si(jsea))*(dera * radius)    ! first approximation ...
#endif
      ELSE
        delx=5.*sqrt(si(ix))*(dera * radius)    ! first approximation ...
        dely=5.*sqrt(si(ix))*(dera * radius)    ! first approximation ...
      ENDIF
    END IF
 
    ik1=1
#ifdef W3_IG1
    ik1=nint(igpars(5))+1
    nspecigstart = nint(igpars(5))*nth
#endif
    DO ik=ik1, nk
      eb  = 0.
      DO ith=1, nth
        eb = eb + a(ith+(ik-1)*nth)
      END DO
      eb   = eb * dden(ik) / cg(ik)
      emeana    = emeana  + eb
      fmean2    = fmean2  + eb /sig(ik)**2
      fmeana    = fmeana  + eb /sig(ik)
    END DO
    fmeana  = tpiinv * (emeana / max( 1.e-7 , fmeana ))
    fmean2  = tpiinv * sqrt(emeana / max( 1.e-7 , fmean2 ))
    fmeana  = max(fmeana,sig(1))
    !
    ! 1.  Sets reflection term to zero
    !
    icalc=1
#ifdef W3_IG1
    stmp1 = 0.
    stmp2 = 0.
#endif
    hs=4.*sqrt(emeana)
#ifdef W3_IG1
    atmp(:)=a(:)      ! the IG energy will be added to this ATMP
    atmp2(:)=a(:)     ! this is really to keep in memory the original spectrum
    IF (igbcoverwrite.AND.reflc(1).GT.0) THEN
      igfac1 = 1.
      atmp2(1:nspecigstart) = 0.
    END IF
    !
    ! resets IG band energy to zero
    !
    DO icalc=1,2
#endif
      s = 0.
#ifdef W3_IG1
      IF (igbcoverwrite) a(1:nspecigstart)=0.
      IF (icalc.EQ.1) a=atmp2
      IF (icalc.EQ.2) THEN
        !
        ! 1.1 Replaces IG part by forced IG
        !
        ! determines highest IG frequency
        !
        IF (igswellmax) THEN
          nkig=nk
        ELSE
          nkig=nint(igpars(5))
        ENDIF
        freqig=sig(nint(igpars(5)))*tpiinv
        !
        nspecig=nkig*nth
        atmp(1:nspecigstart)=0.               ! flat bottom approximation (Hasselmann 1962)
        !                 is not valid for long waves
        IF (nint(igpars(3)).EQ.1) THEN        ! IGPARS(3) = IGSOURCE
          IF (nint(igpars(8)).EQ.1) THEN      ! in this case, uses depth at break point
            depthig=max(1.,hs/0.3)            ! to be modified later with a proper gamma
          ELSE
            depthig=depth
          END IF
          IF (igpars(10).GT.0.) depthig = igpars(10) ! fixed depth for 2nd order calculation
          !
          !/ --- INLINED WAVNU1 (START) ---------------------------------------- /
          !
          DO ik=1, nk
            sqrth  = sqrt(depthig)
            six    = sig(ik) * sqrth
            i1     = int(six/dsie)
            IF (i1.LE.n1max) THEN
              i2 = i1 + 1
              r1 = six/dsie - real(i1)
              r2 = 1. - r1
              wnb(ik) = ( r2*ewn1(i1) + r1*ewn1(i2) ) / depth
              cgb(ik) = ( r2*ecg1(i1) + r1*ecg1(i2) ) * sqrth
            ELSE
              wnb(ik) = sig(ik)*sig(ik)/grav
              cgb(ik) = 0.5 * grav / sig(ik)
            END IF
          END DO
          !
          !/ --- INLINED WAVNU1 (END) ------------------------------------------ /
          !
          IF (nint(igpars(1)).EQ.1) THEN         ! IGPARS(1) = IGMETHOD
            CALL w3addig(atmp,depthig,wnb,cgb,1)
          ELSE
            CALL w3add2ndorder(atmp,depthig,wnb,cgb,1)
          END IF
          ! Transforms energy back to proper depth
          DO ik=1,nkig
            atmp(1+(ik-1)*nth:ik*nth)=atmp(1+(ik-1)*nth:ik*nth)*(cgb(ik)*wn(ik))/(cg(ik)*wnb(ik))
          END DO
          a(1:nspecig)=atmp(1:nspecig)
          IF (igswellmax) THEN
            DO ispec=1,nspecig
              a(ispec)=max(atmp(ispec)-a(ispec),0.)
            END DO
          ELSE
            a(1:nspecig)=atmp(1:nspecig)
          ENDIF
          !
        ELSEIF (nint(igpars(3)).EQ.2) THEN   ! Empirical source of IG energy
          !
          ! This empirical source was adjusted to Waimea and Duck data
          ! When applied to deep water the 1/Depth must be replaced with k/Cg
          ! Hence the proper coefficient is WN(IK)/CG(IK)*GRAV**2/SIG(IK)
          !
          ! The empirical form is HIG = IGEMPIRICAL * ...
          ! this is not quite yet a wave height: multiplied below by MIN(0.0036, ...
          !
          hig= hs/(max(fmean2,freqig)**2)
          efig=(hig*0.25)**2/0.0279  ! this is (HsIG/4)^2 / df
          hig2 = 0.
          DO ik=1,nkig
            !
            ! First approximation: constant IG spectrum  with frequency
            !
            efig1=efig*min(0.0036,igpars(11)*wn(ik)/cg(ik)*grav**2/sig(ik))
            !
            ! Correction: gives a frequency shape ...
            !
            IF (ik.LT.ik1) THEN
              ! The 1.5 exponent of Ardhuin et al. 2014 (see figure 8.a) was probably too high  ... now reduced to 1.0
              efig1=efig1*igfac1*1.2*min(1.,0.013/(sig(ik)*tpiinv))**1.0
            ELSE
              efig1=efig1*igfac2*1.2*min(1.,0.013/(sig(ik)*tpiinv))**1.0
            END IF
            !
            ! Conversion to action spectral density A(k,theta), assuming isotropic dir.
            !
            a(1+(ik-1)*nth:ik*nth)=efig1*cg(ik)/((sig(ik)*tpi)*tpi)
            hig2 = hig2 + efig1*dsii(ik)*tpiinv
          END DO
        ELSE
          nspecig=0
        END IF
        !
      END IF ! ICALC EQ 2
#endif
      !
      nrs=nint(refpars(8))
      IF (refpars(6).GT.0) THEN
        !
        ! This is the Miche parameter for a beach slope of REFLC(3)
        !
        IF(reflc(3)/=reflc(3)) THEN ! isnan test
          slope=0.001
        ELSE
          slope=max(0.001,reflc(3))
        END IF
        michefac=0.0001*grav**2*(slope**5)  &
             /(max(emeana,1e-4)*max(fmeana,0.001)**4)
        ramp0=max(0.07*(alog10(michefac)+4.5)+1.5*michefac,0.005)   ! IF REFLC(1)=1, 0.07 should be 0.007
        ! NB: these constants are adjusted for REFLC(1) = 0.1. If  REFLC(1)=1, 0.07 should be 0.007
      ELSE
        ramp0=1.
      ENDIF
 
      !
      ! 2.  Shoreline reflection =============================================== *
      !
      IF (reflc(1).GT.0) THEN
        fac1=1/(0.5*real(nth))
        fac2=1.57/(0.5*real(nth))
        !           FAC3=2.6/(0.5*REAL(NTH))  ! this is for NRS=4
        fac3=2./sum(abs(ecos(1:nth))**nrs)
        fac4=1.
        !
        DO ik=1, nk
          !
          ! Includes frequency dependence (see Elgar et al. JPO 1994)
          !
          IF (refpars(6).GT.0) THEN
            ramp=(max((0.75*tpi*fmeana/sig(ik)),1.)**refpars(10))*ramp0
            ramp1=min(refpars(9),reflc(1)*ramp)
            ramp2=min(refpars(9),reflc(2)*ramp)
          ELSE
            ramp1=ramp0*reflc(1)
            ramp2=ramp0*reflc(2)
          END IF
          !
          ! Special treatment for unstructured grids when not using source term
          !
          IF (gtype.EQ.ungtype.AND.refpars(3).LT.0.5) THEN
            IF (lpdlib) THEN
#ifdef W3_PDLIB
              iobpdip = iobpd_loc(:,jsea)
#endif
            ELSE
              iobpdip = iobpd(:,ix)
            ENDIF
            DO ith=1, nth
              ispeci=ith+(ik-1)*nth
              a(ispeci)=a(ispeci)*iobpdip(ith) !puts to zero the energy not going to coast
            END DO
            !
            DO ith=1, nth
              r1=ecos(1+mod(abs(ith-refld(1)),nth))
              r1=iobpdip(ith)
              ispeci=ith+(ik-1)*nth
              r2=ramp1*a(ispeci)
              IF (r1.GT.0.AND.r2.GT.0) THEN
                !
                ! Determines direction of specular reflection: th3=pi+2*n-th1
                !
                ith3=1+mod(nth/2+nth+2*refld(1)-ith-1,nth)
                DO ith2=1,nth
                  !
                  !  Adds energy into reflected directions (ITH2)
                  !
                  ispec=ith2+(ik-1)*nth
                  r3=ecos(1+mod(abs(ith2-refld(1)),nth))
                  IF (r3.LT.0) THEN
                    r4=ecos(1+mod(abs(ith2-ith3),nth))*(1-iobpdip(ith2))
                    IF (r4.GT.0.) THEN
                      !
                      ! Tests the type of shoreline geometry
                      !
                      SELECT CASE (refld(2))
                      CASE (0)
                        ! Sharp corner: broad reflection
                        s(ispec)=s(ispec)+r2*fac1/dt
 
                        ! FA: analog to following lines to be swapped in if reflection method changed
                        ! RECT CASE:
                        !  S(ISPEC)=S(ISPEC)+    &
                        !    REAL(REFLD(3))*R2*CG(IK)*ABS(ECOS(ITH2X)/DELX)*FAC1   &
                        !   +REAL(REFLD(4))*R2*CG(IK)*ABS(ESIN(ITH2Y)/DELY)*FAC1
 
 
                      CASE (1)
                        ! mild corner: average reflection
                        s(ispec)=s(ispec)+r2*abs(r4)*fac2/dt
                      CASE (2)
                        ! straight coast: narrow reflection
                        !     IF(ITH3.EQ.ITH2) S(ISPEC)=S(ISPEC)+R2/DT  ! THIS IS FOR SPECULAR REF.
                        s(ispec)=s(ispec)+r2*(r4**nrs) *fac3/dt
                      END SELECT
                    END IF  ! (R4.GT.0.)
                  END IF  ! (R3.LT.0)
                END DO ! ITH2=1,NTH
              END IF  ! (R1.GT.0.AND.R2.GT.0)
            END DO  ! ITH=1, NTH
          ELSE ! (GTYPE.NE.UNGTYPE)
            !
            ! This is for structured grids ....
            !
            !
            !  Loop on  incident wave direction (ITH)
            !
            DO ith=1, nth
              r1=ecos(1+mod(abs(ith-refld(1)),nth))
              r2=ramp1*a(ith+(ik-1)*nth)
              IF (r1.GT.0.AND.r2.GT.0) THEN
                DO ith2=1,nth
                  !
                  !  Adds energy into reflected directions (ITH2)
                  !
                  ispec=ith2+(ik-1)*nth
                  ith2x=1+mod(nth+ith2-refld(5)-1,nth)
                  ith2y=1+mod(nth+ith2-refld(6)-1,nth)
                  r3=ecos(1+mod(abs(ith2-refld(1)),nth))
                  IF (r3.LT.0) THEN
                    !
                    ! Determines direction of specular reflection: th3=pi+2*n-th1
                    !
                    ith3=1+mod(nth/2+nth+2*refld(1)-ith-1,nth)
                    r4=ecos(1+mod(abs(ith2-ith3),nth))
                    IF (r4.GT.0.) THEN
                      !
                      ! Tests the type of shoreline geometry
                      ! NB: REFLD(3) or REFLD(4) is equal to 1 if the reflection is applied (real land neighbor)
                      SELECT CASE (refld(2))
                      CASE (0)
                        ! Sharp corner: broad reflection
                        s(ispec)=s(ispec)+    &
                             REAL(REFLD(3))*R2*CG(IK)*ABS(ECOS(ITH2X)/DELX)*FAC1   &
                             +REAL(REFLD(4))*R2*CG(IK)*ABS(ESIN(ITH2Y)/DELY)*FAC1
                      CASE (1)
                        ! mild corner: average reflection
                        !
                        s(ispec)=s(ispec)+    &
                             REAL(REFLD(3))*R2*CG(IK)*ABS(ECOS(ITH2X)/DELX)*ABS(R4)*FAC2 &
                             + REAL(REFLD(4))*R2*CG(IK)*ABS(ESIN(ITH2Y)/DELY)*ABS(R4)*FAC2
                      CASE (2)
                        ! straight coast: narrow reflection
                        ! Specular for tests
                        !                   S(ISPEC)=S(ISPEC)+REAL(REFLD(3))*R2*CG(IK)*ABS(ECOS(ITH2)/DELX)  &
                        !                                    +REAL(REFLD(4))*R2*CG(IK)*ABS(ESIN(ITH2)/DELY)
                        !
                        s(ispec)=s(ispec)+real(refld(3))*r2*cg(ik)*abs(ecos(ith2x)/delx) &
                             *(r4**nrs) *fac3                          &
                             +real(refld(4))*r2*cg(ik)*abs(esin(ith2y)/dely) &
                             *(r4**nrs) *fac3
                      END SELECT
                    END IF ! (R4.GT.0.)
                  END IF  ! (R3.LT.0)
                END DO  ! ITH2=1,NTH
              END IF  ! (R1.GT.0.AND.R2.GT.0)
            END DO  ! ITH=1,NTH
          END IF  ! UNGTYPE
 
        END DO ! loop on IK
      END IF   ! end of test on REFLC(1)
      !
      !  Add diffuse reflection due to subgrid islands and icebergs
      !  At present this feature is not supported for unstructured grids.
      !
      IF (    ((refpars(2).GT.0.).AND.((trnx+trny).LT.2))     &
           .OR.((refpars(4).GT.0.).AND.(berg.GT.0)       )   ) THEN
        !
        ! Includes frequency dependence (see Elgar et al. JPO 1994)
        !
        IF (refpars(6).GT.0) THEN
          ramp=(max((0.75*tpi*fmeana/sig(ik)),1.)**refpars(10))*ramp0
          ramp2=min(refpars(9),reflc(2)*ramp)
          !
          ! recomputes coefficients for iceberg slope given by REFLC(4)
          !
          slope=max(0.001,reflc(4))
          michefac=0.0001*grav**2*(slope**5)  &
               /(max(emeana,1e-4)*max(fmeana,0.001)**4)
          ramp0=max(0.007*(alog10(michefac)+4.5)+1.5*michefac,0.005)   ! IF REFLC(1)=1, 0.07 should be 0.007
          ramp=(max((0.75*tpi*fmeana/sig(ik)),1.)**refpars(10))*ramp0
          ramp4=min(refpars(9),ramp)
        ELSE
          ramp2=ramp0*reflc(2)
          ramp4=ramp0*reflc(4)
        END IF
        !
        !
        r2x=  ramp2*refpars(2)*max(0.,min(1.,(1-trnx)))  &
             + ramp4*refpars(4)*max(0.,min(1.,(1-exp(-berg*delx*0.0001))))
        r2y=  ramp2*refpars(2)*max(0.,min(1.,(1-trny)))  &
             + ramp4*refpars(4)*max(0.,min(1.,(1-exp(-berg*dely*0.0001))))
        fac1=1/(0.5*real(nth))
        DO ik=1, nk
          DO ith=1, nth
            r2=a(ith+(ik-1)*nth)
            IF (r2.GT.0.) THEN
 
              DO ith2=1,nth
                ispec=ith2+(ik-1)*nth
                r3=ecos(1+mod(nth+ith2-ith,nth))
                IF (r3.LT.0) THEN
                  s(ispec)=s(ispec)+ &
                       cg(ik)*r2x*r2*abs(ecos(ith2)/delx)*fac1 &
                       + cg(ik)*r2y*r2*abs(esin(ith2)/dely)*fac1
                END IF
              END DO
            END IF
          END DO
        END DO
      END IF
 
#ifdef W3_IG1
      IF (icalc.EQ.1) THEN
        stmp1(nspecigstart+1:nspec) = s(nspecigstart+1:nspec)
      ELSE
        stmp2 = s
        DO ispec = 1, nspec
          s(ispec) = max(stmp2(ispec),stmp1(ispec))
        END DO
      END IF
    ENDDO ! ICALC = 1,2
    a(1:nspecig)=atmp2(1:nspecig)   ! removes bound IG components ...
#endif
    !/
    !/ End of W3SREF ----------------------------------------------------- /
    !/
  END SUBROUTINE w3sref!/ ------------------------------------------------------------------- /
 
  !/
  !/ End of module W3REF1MD -------------------------------------------- /
  !/
END MODULE w3ref1md