SUBROUTINE SINPUT_ARD_LOCAL (IPP,F,FL,THWNEW,USNEW,Z0NEW,&
& ROAIRN,WSTAR,SL,XLLWS,SSIN,DSSIN)
! ----------------------------------------------------------------------
!**** *SINPUT* - COMPUTATION OF INPUT SOURCE FUNCTION.
! P.A.E.M. JANSSEN KNMI AUGUST 1990
! OPTIMIZED BY : H. GUENTHER
! MODIFIED BY :
! J-R BIDLOT NOVEMBER 1995
! J-R BIDLOT FEBRUARY 1996-97
! J-R BIDLOT FEBRUARY 1999 : INTRODUCE ICALL AND NCALL
! P.A.E.M. JANSSEN MAY 2000 : INTRODUCE GUSTINESS
! J-R BIDLOT FEBRUARY 2001 : MAKE IT FULLY IMPLICIT BY ONLY
! USING NEW STRESS AND ROUGHNESS.
! S. ABDALLA OCTOBER 2001: INTRODUCTION OF VARIABLE AIR
! DENSITY AND STABILITY-DEPENDENT
! WIND GUSTINESS
! P.A.E.M. JANSSEN OCTOBER 2008: INTRODUCE DAMPING WHEN WAVES ARE
! RUNNING FASTER THAN THE WIND.
! J-R BIDLOT JANUARY 2013: SHALLOW WATER FORMULATION.
! L. AOUF March 2011 : USE OF NEW DISSIPATION DEVELOPED BY ARDHUIN ET AL.2010
!* PURPOSE.
! ---------
! COMPUTE INPUT SOURCE FUNCTION AND STORE ADDITIVELY INTO NET
! SOURCE FUNCTION ARRAY, ALSO COMPUTE FUNCTIONAL DERIVATIVE OF
! INPUT SOURCE FUNCTION.
!
! GUSTINESS IS INTRODUCED FOLL0WING THE APPROACH OF JANSSEN(1986),
! USING A GAUSS-HERMITE APPROXIMATION SUGGESTED BY MILES(1997).
! IN THE PRESENT VERSION ONLY TWO HERMITE POLYNOMIALS ARE UTILISED
! IN THE EVALUATION OF THE PROBABILITY INTEGRAL. EXPLICITELY ONE THEN
! FINDS:
!
! <GAMMA(X)> = 0.5*( GAMMA(X(1+SIG)) + GAMMA(X(1-SIG)) )
!
! WHERE X IS THE FRICTION VELOCITY AND SIG IS THE RELATIVE GUSTINESS
! LEVEL.
!** INTERFACE.
! ----------
! *CALL* *SINPUT (F, FL, IJS, IJL, THWNEW, USNEW, Z0NEW,
! & ROAIRN,WSTAR, SL, XLLWS)
! *F* - SPECTRUM.
! *FL* - DIAGONAL MATRIX OF FUNCTIONAL DERIVATIVE.
! *IJS* - INDEX OF FIRST GRIDPOINT.
! *IJL* - INDEX OF LAST GRIDPOINT.
! *THWNEW* - WIND DIRECTION IN RADIANS IN OCEANOGRAPHIC
! NOTATION (POINTING ANGLE OF WIND VECTOR,
! CLOCKWISE FROM NORTH).
! *USNEW* - NEW FRICTION VELOCITY IN M/S.
! *Z0NEW* - ROUGHNESS LENGTH IN M.
! *ROAIRN* - AIR DENSITY IN KG/M3
! *WSTAR* - FREE CONVECTION VELOCITY SCALE (M/S).
! *SL* - TOTAL SOURCE FUNCTION ARRAY.
! *XLLWS*- 1 WHERE SINPUT IS POSITIVE
! METHOD.
! -------
! SEE REFERENCE.
! EXTERNALS.
! ----------
! WSIGSTAR.
! MODIFICATIONS
! -------------
! - REMOVAL OF CALL TO CRAY SPECIFIC FUNCTIONS EXPHF AND ALOGHF
! BY THEIR STANDARD FORTRAN EQUIVALENT EXP and ALOGHF
! - MODIFIED TO MAKE INTEGRATION SCHEME FULLY IMPLICIT
! - INTRODUCTION OF VARIABLE AIR DENSITY
! - INTRODUCTION OF WIND GUSTINESS
! REFERENCE.
! ----------
! P. JANSSEN, J.P.O., 1989.
! P. JANSSEN, J.P.O., 1991
! ----------------------------------------------------------------------
! USE YOWCOUP , ONLY : BETAMAX ,ZALP ,TAUWSHELTER, XKAPPA
! USE YOWFRED , ONLY : FR ,TH ,DFIM ,COSTH ,SINTH
! USE YOWMPP , ONLY : NINF ,NSUP
! USE YOWPARAM , ONLY : NANG ,NFRE ,NBLO
! USE YOWPCONS , ONLY : G ,ZPI ,ROWATER ,YEPS
! USE YOWSHAL , ONLY : TFAK ,INDEP, DEPTH
! USE YOWSTAT , ONLY : ISHALLO
! USE YOWTABL , ONLY : IAB ,SWELLFT
! 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, &
& DELUST, DELALP, BETAMAX, XKAPPA, IDAMPING, TAUWSHELTER, &
& SINTH, COSTH, ZERO, &
& ROWATER => RHOW, &
& ROAIR => RHOA, &
& TH => SPDIR, &
& DELTH => DDIR, &
& G => G9, &
& ZPI => PI2, &
& EPSMIN => SMALL, &
& NANG => MDC, &
& NFRE => MSC, &
& INDEP => DEP
! ----------------------------------------------------------------------
IMPLICIT NONE
INTEGER, INTENT(IN) :: IPP
INTEGER :: K,M,IND
REAL, PARAMETER :: ABMIN = 0.3
REAL, PARAMETER :: ABMAX = 8.
REAL(rkind) :: CONST1
REAL(rkind) :: X1,X2,ZLOG,ZLOG1,ZLOG2,ZLOG2X,XV1,XV2,ZBETA1,ZBETA2
REAL(rkind) :: XI,X,DELI1,DELI2
REAL(rkind) :: FU,FUD,FW,NU_AIR,SWELLFPAR,SWELLF,SWELLF2,SWELLF3,SWELLF4,SWELLF5
REAL(rkind) :: SWELLF7, SMOOTH, YEPS
REAL(rkind) :: ARG, DELAB, CONST11, CONST22
!REAL(KIND=JPRB) :: ZHOOK_HANDLE
REAL(rkind), DIMENSION(NANG) :: PP
REAL(rkind), DIMENSION(NFRE) :: FAC, CONST, SIG, CONSTF
REAL(rkind) :: THWNEW, USNEW, Z0NEW, ROAIRN, WSTAR
REAL(rkind) :: TAUX, TAUY, TAUPX,TAUPY,USTP,USDIRP
REAL(rkind) :: Z0VIS, Z0NOZ, FWW
REAL(rkind) :: PVISC, PTURB
REAL(rkind) :: UCN1, UCN2, ZCN, CM
REAL(rkind) :: SIG_N, UORBT, AORB, TEMP, RE, ZORB
REAL(rkind) :: CNSN
REAL(rkind) :: XSTRESS, YSTRESS
REAL(rkind), DIMENSION(NANG) :: TEMP1, UFAC2
REAL(rkind), DIMENSION(NFRE) :: XK
REAL(rkind), DIMENSION(NANG) :: DSTAB1, DSTAB2
REAL(rkind), DIMENSION(NANG,NFRE) :: F,FL,SL
REAL(rkind), DIMENSION(NANG,NFRE) :: DSTAB
REAL(rkind), DIMENSION(NANG,NFRE) :: XLLWS
REAL(rkind), DIMENSION(NANG,NFRE) :: SSIN, DSSIN
!IF (LHOOK) CALL DR_HOOK('SINPUT',0,ZHOOK_HANDLE)
! ----------------------------------------------------------------------
CONST1 = BETAMAX/XKAPPA**2 /ROWATER
NU_AIR = 1.4E-5
SWELLFPAR = 3.
SWELLF = 0.8
SWELLF2 = -0.018
SWELLF3 = 0.015
SWELLF4 = 1.E5
SWELLF5 = 1.2
SWELLF7 = 2.3E5
FU=ABS(SWELLF3)
FUD=SWELLF2
FW=MAX(ABS(SWELLF3),0.)
DELAB= (ABMAX-ABMIN)/REAL(IAB)
!
!* 1. PRECALCULATED ANGULAR DEPENDENCE.
! ---------------------------------
DO K=1,NANG
TEMP1(K) = COS(TH(K)-THWNEW)
ENDDO
! ESTIMATE THE STANDARD DEVIATION OF GUSTINESS.
CALL WSIGSTAR (IPP, IPP, USNEW, Z0NEW, WSTAR, SIG_N)
! ----------------------------------------------------------------------
! computation of Uorb and Aorb
! ----------------------------------------------------------------------
UORBT = ZERO
AORB = ZERO
DO M=1,NFRE
K=1
SIG(M) = ZPI*FR(M)
TEMP = F(K,M)
DO K=2,NANG
TEMP = TEMP+F(K,M)
ENDDO
UORBT = UORBT+DFIM(M)*(SIG(M)**2)*TEMP
AORB = AORB+DFIM(M)*TEMP
ENDDO
UORBT = 2*SQRT(UORBT) ! this is the significant orbital amplitude
AORB = 2*SQRT(AORB)
Z0VIS = 0.1*NU_AIR/MAX(USNEW,0.0001)
Z0NOZ = max(Z0VIS,0.04*Z0NEW)
ZORB = AORB/Z0NOZ
IF (UORBT.NE.ZERO) THEN
RE = 4*UORBT*AORB / NU_AIR ! this is the Reynolds number
ELSE
RE = ZERO
ENDIF
! calcul fww
FU=ABS(SWELLF3)
FUD=SWELLF2
XI=(LOG10(MAX(ZORB,3.)) -ABMIN)/DELAB
IND = MIN (IAB-1, INT(XI))
DELI1= MIN (1. ,XI-FLOAT(IND))
DELI2= 1. - DELI1
FWW =SWELLFT(IND)*DELI2+SWELLFT(IND+1)*DELI1
XSTRESS=ZERO
YSTRESS=ZERO
TAUX=(USNEW**2)*SIN(THWNEW)
TAUY=(USNEW**2)*COS(THWNEW)
!* 2. LOOP OVER FREQUENCIES.
! ----------------------
DO M=1,NFRE
TAUPX=TAUX-ABS(TAUWSHELTER)*XSTRESS
TAUPY=TAUY-ABS(TAUWSHELTER)*YSTRESS
USTP=(TAUPX**2+TAUPY**2)**0.25
USDIRP=ATAN2(TAUPX,TAUPY)
CONSTF(M) =ZPI*ROWATER*FR(M)*DFIM(M)
FAC(M) = ZPI*FR(M)
CONST(M)=FAC(M)*CONST1
DO K=1,NANG
TEMP1(K) = COS(TH(K)-USDIRP)
ENDDO
!* INVERSE OF PHASE VELOCITIES AND WAVE NUMBER.
! -------------------------------------------
IF (ISHALLO.EQ.1) THEN
CM = FAC(M)/G
XK(M) = ((ZPI*FR(M))**2)/G
ELSE
CM = WK(M,IPP)/FAC(M)
XK(M) = WK(M,IPP)
ENDIF
!* PRECALCULATE FREQUENCY DEPENDENCE.
! ----------------------------------
UCN1 = USTP*(1.+SIG_N)*CM+ ZALP
UCN2 = USTP*(1.-SIG_N)*CM+ ZALP
ZCN = LOG(G*Z0NEW*CM**2)
CNSN = CONST(M) * ROAIRN
!* 2.1 LOOP OVER DIRECTIONS.
! ---------------------
DO K=1,NANG
IF (TEMP1(K).GT.0.01) THEN
X = TEMP1(K)*UCN1
ZLOG = ZCN + XKAPPA/X
IF (ZLOG.LT.0.) THEN
ZLOG2X=ZLOG*ZLOG*X
UFAC2(K) = EXP(ZLOG)*ZLOG2X*ZLOG2X
XLLWS(K,M)= 1.
ELSE
UFAC2(K) = 0.
XLLWS(K,M)= 0.
ENDIF
X = TEMP1(K)*UCN2
ZLOG = ZCN + XKAPPA/X
IF (ZLOG.LT.0.) THEN
ZLOG2X=ZLOG*ZLOG*X
UFAC2(K) = UFAC2(K)+EXP(ZLOG)*ZLOG2X*ZLOG2X
XLLWS(K,M)= 1.
ENDIF
ELSE
UFAC2(K) = 0.
XLLWS(K,M)=0.
ENDIF
ENDDO
! SWELL DAMPING:
DO K=1,NANG
PP(K) = 1.
YEPS = ROAIR/ROWATER
DSTAB1(K) = -SWELLF5*YEPS*2*XK(M)*SQRT(2*NU_AIR*SIG(M))*PP(K)
FW = 0.04*ZORB**(-0.25)
DSTAB2(K) = -YEPS*SWELLF*(FWW*UORBT+(FU+FUD*TEMP1(K))*USTP) &
*16*SIG(M)**2/G
END DO
IF (SWELLF7.GT.0.) THEN
SMOOTH=0.5*TANH((RE-SWELLF4)/SWELLF7)
PTURB=0.5+SMOOTH
PVISC=0.5-SMOOTH
DO K=1,NANG
DSTAB(K,M) = PVISC*DSTAB1(K)+PTURB*DSTAB2(K)
END DO
ELSE
PTURB=0.5
PVISC=0.5
IF (RE.LE.SWELLF4) THEN
DO K=1,NANG
DSTAB(K,M) = PVISC*DSTAB1(K)
END DO
ELSE
DO K=1,NANG
DSTAB(K,M) = PTURB*DSTAB2(K)
END DO
ENDIF
ENDIF
!* 2.2 ADDING INPUT SOURCE TERM TO NET SOURCE FUNCTION.
! AND UPDATE THE SHELTERING STRESS.
! ------------------------------------------------
DO K=1,NANG
CONST11=CONSTF(M)*SINTH(K)
CONST22=CONSTF(M)*COSTH(K)
FL(K,M) = 0.5*CNSN*UFAC2(K)+DSTAB(K,M)
SL(K,M) = FL(K,M)*F(K,M)
SSIN(K,M) = FL(K,M)*F(K,M)
DSSIN(K,M) = FL(K,M)
XSTRESS=XSTRESS+SL(K,M)*CONST11/MAX(ROAIRN,1.)
YSTRESS=YSTRESS+SL(K,M)*CONST22/MAX(ROAIRN,1.)
ENDDO
ENDDO
!IF (LHOOK) CALL DR_HOOK('SINPUT_ARD',1,ZHOOK_HANDLE)
RETURN
END SUBROUTINE SINPUT_ARD_LOCAL