SUBROUTINE SDISS_ARDH_VEC_LOCAL (IPP, F, FL, SL, F1MEAN, XKMEAN, &
& PHIEPS, TAUWD, M, SSDS, DSSDS)
! ----------------------------------------------------------------------
!**** *SDISSIP_ARDH_VEC* - COMPUTATION OF DISSIPATION SOURCE FUNCTION.
! LOTFI AOUF METEO FRANCE 2013
! FABRICE ARDHUIN IFREMER 2013
!* PURPOSE.
! --------
! COMPUTE DISSIPATION SOURCE FUNCTION AND STORE ADDITIVELY INTO
! NET SOURCE FUNCTION ARRAY. ALSO COMPUTE FUNCTIONAL DERIVATIVE
! OF DISSIPATION SOURCE FUNCTION.
!** INTERFACE.
! ----------
! *CALL* *SDISSIP (F, FL, IJS, IJL, SL, F1MEAN, XKMEAN,*
! PHIEPS, TAUWD, M LCFLX)
! *F* - SPECTRUM.
! *FL* - DIAGONAL MATRIX OF FUNCTIONAL DERIVATIVE
! *IJS* - INDEX OF FIRST GRIDPOINT
! *IJL* - INDEX OF LAST GRIDPOINT
! *SL* - TOTAL SOURCE FUNCTION ARRAY
! *F1MEAN* - MEAN FREQUENCY BASED ON 1st MOMENT.
! *XKMEAN* - MEAN WAVE NUMBER BASED ON 1st MOMENT.
! *PHIEPS* - ENERGY FLUX FROM WAVES TO OCEAN INTEGRATED OVER
! THE PROGNOSTIC RANGE.
! *TAUWD* - DISSIPATION STRESS INTEGRATED OVER
! THE PROGNOSTIC RANGE.
! *MIJ* - LAST FREQUENCY INDEX OF THE PROGNOSTIC RANGE.
! *LCFLX* - TRUE IF THE CALCULATION FOR THE FLUXES ARE
! PERFORMED.
! METHOD.
! -------
! SEE REFERENCES.
! EXTERNALS.
! ----------
! NONE.
! REFERENCE.
! ----------
! ARDHUIN et AL. JPO DOI:10.1175/20110JPO4324.1
! ----------------------------------------------------------------------
!USE YOWFRED , ONLY : FR, TH, FRATIO, DELTH, GOM, COSTH, SINTH, DFIM
!USE YOWPCONS , ONLY : RAD ,G ,ZPI ,ROWATER ,YEPS
!USE YOWMEAN , ONLY : EMEAN, FMEAN
!USE YOWMPP , ONLY : NINF ,NSUP
!USE YOWPARAM , ONLY : NANG ,NFRE
!USE YOWSHAL , ONLY : TFAK ,INDEP, TCGOND
!USE YOWSPEC , ONLY : U10NEW, THWNEW, USNEW
!USE YOWSTAT , ONLY : ISHALLO
!USE YOWTEST , ONLY : IU06 ,ITEST
!USE PARKIND1 ,ONLY : JPIM ,JPRB
!USE YOMHOOK ,ONLY : LHOOK ,DR_HOOK
USE DATAPOOL, ONLY : FR, WETAIL, FRTAIL, WP1TAIL, ISHALLO, RKIND, &
& DFIM, DFIMOFR, DFFR, DFFR2, WK, ZALP, IAB, SWELLFT, &
& IUSTAR, IALPHA, USTARM, TAUT, ONETHIRD, RKIND, ONE, &
& DELUST, DELALP, BETAMAX, XKAPPA, IDAMPING, TAUWSHELTER, &
& FRATIO, EMEAN, USNEW, THWNEW, DEGRAD, LCFLX, MSC, MDC, &
& SINTH, COSTH, &
& ROWATER => RHOW, &
& ROAIR => RHOA, &
& TH => SPDIR, &
& DELTH => DDIR, &
& G => G9, &
& ZPI => PI2, &
& EPSMIN => SMALL, &
& NANG => MDC, &
& NFRE => MSC, &
& INDEP => DEP, &
& CG0 => CG, &
& WK0 => WK
! ----------------------------------------------------------------------
IMPLICIT NONE
INTEGER :: K , M, IPP
INTEGER :: I, J, I1, J1
INTEGER :: IS, SDSNTH, DIKCUMUL
INTEGER :: I_INT, J_INT, M2, K2
INTEGER :: NSMOOTH(NFRE)
INTEGER :: ISDSDTH
INTEGER :: ISSDSBRFDF
INTEGER :: MIJ
INTEGER, ALLOCATABLE :: SATINDICES(:,:)
REAL(rkind) :: XK(NFRE), CG(NFRE)
REAL(rkind) :: ALFAMEAN
REAL(rkind) :: TPIINV, TMP00, TMP01, TMP02, TMP03, TMP04
REAL(rkind) :: COSWIND
REAL(rkind) :: DTURB
REAL(rkind) :: EPSR
REAL(rkind) :: W, P0, MICHE, DELTA1, DELTA2
REAL(rkind) :: SCDFM
REAL(rkind) :: SDISS, YEPS
REAL(rkind) :: SDSBR, SDSBR2
REAL(rkind) :: SATURATION2,FACSAT
REAL(rkind) :: SSDSC3, SSDSC4, SSDSC5, SSDSC6, SDSCOS
REAL(rkind) :: SSDSHF, SSDSLF, X, DTEMP, TEMP
REAL(rkind) :: SSDSBRF1, XFR, SXFR,SSDSBR,SDSC1,SSDSP,SSDSC2
REAL(rkind) :: DSIP_
REAL(rkind) :: ROG
REAL(rkind) :: SSDSC1, SSDSC , SSDSISO
!REAL(KIND=JPRB) :: ZHOOK_HANDLE
REAL(rkind) :: SIG(NFRE)
REAL(rkind) :: BSIGBAJ
REAL(rkind) :: F1MEAN, XKMEAN, WNMEAN2
REAL(rkind) :: PHIEPS, TAUWD, CM
REAL(rkind) :: FACTOR, FACTURB
REAL(rkind) :: RENEWALFREQ(NANG)
REAL(rkind), DIMENSION(NFRE) :: CONSTFM
REAL(rkind) :: BTH0(NFRE) !saturation spectrum
REAL(rkind) :: BTH0S(NFRE) !smoothed saturation spectrum
REAL(rkind) :: BTH(NANG,NFRE) !saturation spectrum
REAL(rkind) :: BTHS(NANG,NFRE) !smoothed saturation spectrum
REAL(rkind),DIMENSION(NANG,NFRE) :: F,FL,SL,A, D
REAL(rkind),DIMENSION(NANG,NFRE) :: SSDS, DSSDS
REAL(rkind), ALLOCATABLE, DIMENSION(:) :: C_,C_C,C2_,C2_C2,DSIP_05_C2
REAL(rkind) , ALLOCATABLE :: SATWEIGHTS(:,:)
REAL(rkind) , ALLOCATABLE :: CUMULW(:,:,:,:)
LOGICAL :: LLTEST
!IF (LHOOK) CALL DR_HOOK('SDISSIP_ARDH_VEC',0,ZHOOK_HANDLE)
! ----------------------------------------------------------------------
W=26.
TPIINV = 1./ZPI
MIJ = MSC
ROG = ROWATER*G
IF (LCFLX) THEN
PHIEPS = 0.
TAUWD = 0.
! !!!! CONSTFM is only defined up to M=MIJ
SCDFM = 0.5*DELTH*(1.-1./FRATIO)
DO M=1,MIJ-1
CONSTFM(M) = ROG*DFIM(M)
ENDDO
CONSTFM = ROG*SCDFM*FR(MIJ)
ENDIF
XFR = 1.1 !! ??? JEAN BIDLOT: what is XFR is it FRATIO ????
SSDSBRF1 = 0.5
SXFR = 0.5*(FRATIO-1/FRATIO)
! SDSBR = 1.40E-3
! SDSBR = 1.20E-3 ! from Babanin (personal communication)
SDSBR = 9.0E-4
SDSBR2 = 0.8
SSDSBR = SDSBR
! SDSC1 = -2.1 !! This is Bidlot et al. 2005, Otherwise WAM4 uses -4.5
SDSC1 = 0.
SSDSC1 = SDSC1
SSDSC5 = 0.
SSDSISO = 1.
SSDSP = 2.
! SSDSC2 = -3.0E-5 ! Retuned by Fabrice from VdW
SSDSC2 = -2.8E-5
SSDSC3 = 0.8
SSDSC4 = 1.
SSDSC6 = 0.3
ISDSDTH = 80
SDSCOS = 2.
SSDSHF = 1.
SSDSLF = 1.
EPSR=SQRT(SSDSBR)
TMP00 = SSDSC1*ZPI
YEPS = ROAIR/ROWATER
TMP01 = SSDSC5/G*YEPS
DO M=1, NFRE
SIG(M) = ZPI*FR(M)
END DO
ALFAMEAN = (XKMEAN**2)*EMEAN(IPP)
FACTOR = TMP00 *F1MEAN*(ALFAMEAN*ALFAMEAN)
FACTURB = TMP01*USNEW(IPP)*USNEW(IPP)
WNMEAN2 = MAX( 1.E-10 , XKMEAN )
IF (ISHALLO.EQ.0) THEN
DO M=1, NFRE
XK(M) = WK0(M,IPP)!TFAK(INDEP,M)
CG(M) = CG0(M,IPP)!TCGOND(INDEP,M)
END DO
ELSE
DO M=1, NFRE
XK(M) = (SIG(M)**2)/G
CG(M) = G/(2*ZPI*FR(M))
END DO
ENDIF
DO M=1, NFRE
!cdir outerunroll=4
DO K=1, NANG
A(K,M) = TPIINV*CG(M)*F(K,M)
ENDDO
END DO
IF (ISDSDTH.NE.180) THEN
SDSNTH = MIN(NINT(ISDSDTH*DEGRAD/(DELTH)),NANG/2-1)
ALLOCATE(SATINDICES(NANG,SDSNTH*2+1))
ALLOCATE(SATWEIGHTS(NANG,SDSNTH*2+1))
DO K=1,NANG
DO I_INT=K-SDSNTH, K+SDSNTH
J_INT=I_INT
IF (I_INT.LT.1) J_INT=I_INT+NANG
IF (I_INT.GT.NANG) J_INT=I_INT-NANG
SATINDICES(K,I_INT-(K-SDSNTH)+1)=J_INT
SATWEIGHTS(K,I_INT-(K-SDSNTH)+1)=COS(TH(K)-TH(J_INT))**SDSCOS
END DO
END DO
END IF
! calcul de cumulw
SSDSBRF1 = 0.5
SXFR = 0.5*(FRATIO-1/FRATIO)
! initialise CUMULW
ALLOCATE(CUMULW(NFRE,NANG,NFRE,NANG))
DO I=1,NFRE
DO J=1,NANG
DO I1=1,NFRE
DO J1=1,NANG
CUMULW(I,J,I1,J1)=0.
ENDDO
ENDDO
ENDDO
END DO
IF (SSDSC3.NE.0.) THEN
! DIKCUMUL is the integer difference in frequency bands
! between the "large breakers" and short "wiped-out waves"
DIKCUMUL = NINT(SSDSBRF1/(XFR-1.))
ALLOCATE(C_(NFRE))
ALLOCATE(C_C(NFRE))
ALLOCATE(C2_(NFRE-DIKCUMUL))
ALLOCATE(C2_C2(NFRE-DIKCUMUL))
ALLOCATE(DSIP_05_C2(NFRE-DIKCUMUL))
DO M=1,NFRE
C_(M)=G/SIG(M) ! Valid in deep water only
C_C(M)=C_(M)*C_(M)
END DO
DO M=1,NFRE
DO K=1,NANG
DO M2=1,M-DIKCUMUL
C2_(M2)=G/SIG(M2)
C2_C2(M2)=C2_(M2)*C2_(M2)
DSIP_ = SIG(M2)*SXFR
DSIP_05_C2(M2)=DSIP_/(0.5*C2_(M2))
DO K2=1,NANG
CUMULW(M,K,M2,K2)=SQRT(C_C(M)+C2_C2(M2) &
& -2*C_(M)*C2_(M2)*COSTH(1+ABS(K2-K)))*DSIP_05_C2(M2)
END DO
END DO
END DO
END DO
DEALLOCATE(C_)
DEALLOCATE(C_C)
DEALLOCATE(C2_)
DEALLOCATE(C2_C2)
DEALLOCATE(DSIP_05_C2)
! Multiplies by lambda(k,theta)=1/(2*pi**2) and
! and the coefficient that transforms SQRT(B) to Banner et al. (2000)'s epsilon
! 2.26 is equal to 5.55 (Banner & al. 2000) times 1.6**2 / 2pi where
! 1.6 is the ratio between Banner's epsilon and SQRT(B)
TMP02 = 2*TPIINV*2.26
DO I=1,NFRE
DO J=1,NANG
DO I1=1,NFRE
DO J1=1,NANG
CUMULW(I,J,I1,J1)=CUMULW(I,J,I1,J1)*TMP02
END DO
ENDDO
END DO
END DO
END IF
DO M=1, NFRE
FACSAT=(XK(M)**3)*DELTH
BTH0(M)=SUM(A(1:NANG,M))*FACSAT
END DO
! DO K=1,NANG
DO M=1, NFRE
DO K=1,NANG
FACSAT=(XK(M)**3)*DELTH
! integrates around full circle
BTH(K,M) = SUM(SATWEIGHTS(K,1:SDSNTH*2+1)*A(SATINDICES(K,1:SDSNTH*2+1),M))*FACSAT
END DO
BTH0(M) = MAXVAL(BTH(1:NANG,M))
END DO
!/ST3 SDSBR = 1.20E-3 ! Babanin (personnal communication)
ISSDSBRFDF = 22 ! test pour DC
ISSDSBRFDF = 0
!/ST3 SDSBRF1 = 0.5
!/ST3 SDSBRF2 = 0.
IF (ISSDSBRFDF.GT.0.AND.ISSDSBRFDF.LT.NFRE/2) THEN
!cdir collapse
BTH0S=BTH0
!cdir collapse
NSMOOTH=ONE
!cdir collapse
BTHS=BTH
!cdir outerunroll=4
DO M=1, ISSDSBRFDF
BTH0S (1+ISSDSBRFDF)=BTH0S (1+ISSDSBRFDF)+BTH0(M)
NSMOOTH(1+ISSDSBRFDF)=NSMOOTH(1+ISSDSBRFDF)+1
!cdir collapse
DO K=1,NANG
BTHS(K,M)=BTHS(K,M)+BTH(K,M)
END DO
ENDDO
DO M=2+ISSDSBRFDF,1+2*ISSDSBRFDF
!cdir nodep
BTH0S (1+ISSDSBRFDF)=BTH0S (1+ISSDSBRFDF)+BTH0(M)
NSMOOTH(1+ISSDSBRFDF)=NSMOOTH(1+ISSDSBRFDF)+1
!cdir collapse
DO K=1,NANG
BTHS(K,M)=BTHS(K,M)+BTH(K,M)
END DO
ENDDO
DO M=ISSDSBRFDF,1,-1
!cdir nodep
BTH0S (M)=BTH0S (M+1)-BTH0(M+ISSDSBRFDF+1)
NSMOOTH(M)=NSMOOTH(M+1)-1
!cdir collapse
DO K=1,NANG
BTHS(K,M)=BTHS(K,M)-BTH(K,M)
END DO
ENDDO
DO M=2+ISSDSBRFDF,NFRE-ISSDSBRFDF
!cdir nodep
BTH0S (M)=BTH0S (M-1)-BTH0(M-ISSDSBRFDF-1)+BTH0(M+ISSDSBRFDF)
NSMOOTH(M)=NSMOOTH(M-1)
!cdir collapse
DO K=1,NANG
BTHS(K,M)=BTHS(K,M)-BTH(K,M)+BTH(K,M)
END DO
ENDDO
!cdir novector
DO M=NFRE-ISSDSBRFDF+1,NFRE
!cdir nodep
BTH0S (M)=BTH0S (M-1)-BTH0(M-ISSDSBRFDF)
NSMOOTH(M)=NSMOOTH(M-1)-1
!cdir collapse
DO K=1,NANG
BTHS(K,M)=BTHS(K,M)-BTH(K,M)
END DO
END DO
! final division by NSMOOTH
!cdir collapse
DO M=1,NFRE
BTH0(M)=MAX(0.,BTH0S(M)/NSMOOTH(M))
END DO
DO M=1,NFRE
!cdir outerunroll=4
DO K=1,NANG
BTH(K,M)=MAX(0.,BTHS(K,M)/NSMOOTH(M))
END DO
END DO
END IF
! DELTA1 = 0.4
! DELTA2 = 0.6
DELTA1 = 0.
DELTA2 = 0.
MICHE = 1.0
TMP03 = 1.0/(SSDSBR*MICHE)
DO M=1, NFRE
LLTEST = (SSDSC3.NE.0.AND.M.GT.DIKCUMUL)
IF (XKMEAN.NE.0) THEN
X = WK0(M,IPP)/XKMEAN!TFAK(INDEP(IJ),M)/XKMEAN(IJ)
BSIGBAJ = FACTOR*( (1.-DELTA2)*X + DELTA2*X**2)
ELSE
BSIGBAJ = 0
ENDIF
IF (ISHALLO.EQ.0) THEN
CM=WK0(M,IPP)/SIG(M)!TFAK(INDEP(IJ),M)/SIG(M)
ELSE
!AR: below is a nice bug always the last index ... in the original, should be IJ
CM=SIG(M)/G
ENDIF
DO K=1,NANG
RENEWALFREQ(K)=0.
! Correction of saturation level for shallow-water kinematics
! Cumulative effect based on lambda (breaking probability is
! the expected rate of sweeping by larger breaking waves)
IF (LLTEST) THEN
DO M2=1,M-DIKCUMUL
!AR: below the index is wrong ...
DO K2=1,NANG
IF (BTH0(M2).GT.SSDSBR) THEN
! Integrates over frequencies M2 and directions K2 to
! Integration is performed from M2=1 to a frequency lower than M: IK-DIKCUMUL
RENEWALFREQ(K)=RENEWALFREQ(K)+ CUMULW(M,K,M2,K2) &
& *(MAX(SQRT(BTH(K2,M2))-EPSR,0.))**2
ENDIF
END DO
END DO
ENDIF
SATURATION2=TANH(10*(((BTH(K,M)/SSDSBR)**0.5)-SDSBR2))
COSWIND=(COSTH(K)*COS(THWNEW(IPP))+SINTH(K)*SIN(THWNEW(IPP))) ! �vérifier K ?
DTURB=-2.*SIG(M)*XK(M)*FACTURB*COSWIND ! Theory -> stress direction
P0=SSDSP ! -0.5*SSDSC3*(1-TANH(W*USTAR*XK(M)/SIG(M)-0.1)) ! for SDSC3=1 this is vdW et al.
TMP04 = SSDSC3*RENEWALFREQ(K)
! DTEMP=SSDSC2 * SIG(M) &
DTEMP=SSDSC2 * SIG(M) &
& * ( SSDSC6 *(MAX(0.,BTH0( M)*TMP03-SSDSC4))**P0 &
& + (1-SSDSC6)*(MAX(0.,BTH (K,M)*TMP03-SSDSC4))**P0)&
& - (TMP04+DTURB) !terme cumulatif
D(K,M) = DTEMP + BSIGBAJ*SSDSLF *0.5*(1-SATURATION2) &
& + BSIGBAJ*SSDSHF *0.5*(SATURATION2+1)
WRITE(111116,'(10F15.8)') D(K,M),DTEMP,SIG(M),RENEWALFREQ(K),BTH0(M),BTH(K,M)
END DO
!cdir outerunroll=4
DO K=1, NANG
SL(K,M) = SL(K,M)+D(K,M)*F(K,M)
FL(K,M) = FL(K,M)+D(K,M)
SSDS(K,M) = D(K,M)*F(K,M)
DSSDS(K,M) = D(K,M)
IF (LCFLX.AND.M.LE.MIJ) THEN
SDISS = D(K,M)*F(K,M)
PHIEPS = PHIEPS+SDISS*CONSTFM(M)
TAUWD = TAUWD+CM*SDISS*CONSTFM(M)
ENDIF
END DO
END DO
IF (ALLOCATED(CUMULW)) DEALLOCATE (CUMULW)
IF (ALLOCATED(SATWEIGHTS)) DEALLOCATE (SATWEIGHTS)
IF (ALLOCATED(SATINDICES)) DEALLOCATE (SATINDICES)
!IF (LHOOK) CALL DR_HOOK('SDISSIP_ARDH_VEC',1,ZHOOK_HANDLE)
RETURN
END SUBROUTINE SDISS_ARDH_VEC_LOCAL