#include "wwm_functions.h" !********************************************************************** !* * !********************************************************************** SUBROUTINE SET_WIND(IP,WIND10,WINDTH) ! USE DATAPOOL IMPLICIT NONE INTEGER, INTENT(IN) :: IP REAL(rkind), INTENT(OUT) :: WIND10, WINDTH REAL(rkind) :: WINDX, WINDY, VEC2RAD WINDX = WINDXY(IP,1) WINDY = WINDXY(IP,2) WIND10 = MAX(TWO,SQRT(WINDX**2+WINDY**2)) * WINDFAC WINDTH = VEC2RAD(WINDX,WINDY) END SUBROUTINE SET_WIND !********************************************************************** !* * !********************************************************************** SUBROUTINE SET_FRICTION(IP,ACLOC,WIND10,WINDTH,FPM) ! ! Friction Velocities different formulations .... ! USE DATAPOOL IMPLICIT NONE INTEGER, INTENT(IN) :: IP INTEGER :: I REAL(rkind) , INTENT(IN) :: ACLOC(MSC,MDC) REAL(rkind) , INTENT(IN) :: WIND10, WINDTH REAL(rkind) , INTENT(OUT) :: FPM REAL(rkind) :: WINDX, WINDY REAL(rkind) :: CDRAG REAL(rkind) :: VEC2RAD REAL(rkind) :: EPS_D REAL(rkind) :: z00, z0_t, fU10, CD10 REAL(rkind) :: ULur, Ur , Lur REAL(rkind) :: UFRIC1, UFRIC2 REAL(rkind) :: TM_W, CGP_W, CP_W, TP_W, LP_W, HS_W, KP_W INTEGER, PARAMETER :: MAXITER_WIND = 100 REAL(rkind), PARAMETER :: KAPPA = 0.4_rkind REAL(rkind), PARAMETER :: GAMMA = 0.8_rkind REAL(rkind), PARAMETER :: CZ0T = 0.0075_rkind REAL(rkind), PARAMETER :: EPS_B = 0.5_rkind REAL(rkind), PARAMETER :: EPS_T = 0.24_rkind REAL(rkind) :: VISK SELECT CASE (IFRIC) CASE (1) IF (WIND10 >= 7.5_rkind) THEN CDRAG = (0.8_rkind+0.065_rkind*WIND10)*0.001_rkind ELSE CDRAG = 0.0012873_rkind ENDIF UFRIC(IP) = SQRT(CDRAG)*WIND10 UFRIC(IP) = MAX(1.0E-15_rkind,UFRIC(IP)) CD(IP) = CDRAG FPM = G9 / ( 28.0_rkind * UFRIC(IP) ) Z0(IP) = 10.0_rkind/EXP(KAPPA*WIND10 /UFRIC(IP)) ALPHA_CH(IP)=G9 * Z0(IP) /(UFRIC(IP)**2) TAUTOT(IP)=(UFRIC(IP)**2)*rhoa !write(*,'(i10,3F15.6)') ip, wind10, cd(ip), ufric(ip) CASE (2) UFRIC(IP) = WIND10 * 1.0_rkind / ( 40.0_rkind - 6.0_rkind * LOG(WIND10) ) UFRIC(IP) = MAX(1.0E-15_rkind,UFRIC(IP)) FPM = G9 / ( 28.0_rkind * UFRIC(IP) ) CD(IP) = (UFRIC(IP)/WIND10)**2 Z0(IP) = 10.0_rkind/EXP(KAPPA*WIND10 /UFRIC(IP)) ALPHA_CH(IP)=G9 * Z0(IP) /(UFRIC(IP)**2) TAUTOT(IP)=(UFRIC(IP)**2)*rhoa CASE (3) IF (WIND10 .GT. 1.0_rkind) THEN CD10 = (0.8_rkind + 0.065_rkind * WIND10) * 10E-3_rkind ELSE CD10 = 0.0_rkind END IF UFRIC(IP) = SQRT(CD10) * WIND10 UFRIC(IP) = MAX(1.0E-15_rkind,UFRIC(IP)) CD(IP) = CD10 FPM = G9 / ( 28.0_rkind * UFRIC(IP) ) Z0(IP) = 10.0_rkind/EXP(KAPPA*WIND10 /UFRIC(IP)) ALPHA_CH(IP)=G9 * Z0(IP) /(UFRIC(IP)**2) TAUTOT(IP)=(UFRIC(IP)**2)*rhoa CASE (4) UFRIC(IP) = WIND10 / 28.0_rkind ! First Guess VISK = 1.5E-5_rkind CALL WINDSEASWELLSEP( IP, ACLOC, TM_W, CGP_W, CP_W, TP_W, LP_W, HS_W, KP_W ) IF (HS_W .LT. THR) GOTO 101 EPS_D = MAX(THR,0.5_rkind*HS_W*KP_W) Lur = LP_W/(4.0_rkind*PI) UFRIC2 = ZERO UFRIC1 = UFRIC(IP) ! DO I = 1, MAXITER_WIND IF (I .GT. 1) UFRIC1 = UFRIC2 fU10 = 0.02_rkind * MAX(ZERO, MyTANH(0.075_rkind*WIND10 - 0.75_rkind)) ! Eq. 6 z0_t = MAX( THR, 0.1_rkind*(VISK/MAX(THR,UFRIC1)) + ( CZ0T + fU10 ) * UFRIC1**2/G9 ) ! Eq. 7 TAUHF(IP) = (KAPPA**2*WIND10**2) / LOG(TEN/z0_t)**2 ! Eq. 8 z00 = TEN * EXP( -( KAPPA*WIND10 / MAX(THR,UFRIC1) ) ) ! Estimate the roughness length according the log. profile ULur = UFRIC1/KAPPA * LOG ((EPS_B/KP_w)/MAX(THR,z00)) ! Estimate the velocitiy in the height of reference Ur = MAX (ZERO, ULur - CP_W) ! Estimate the effective velocity TAUW(IP) = EPS_B*GAMMA/PI2 * Ur**2 * EXP(-EPS_T**2/EPS_D**2) ! Stress due to AFS of dominant waves in the wind sea Eq. 9 ! WRITE(BG%FHNDL,*) EPS_D**2, HS_W, KP_W, EXP(-EPS_T**2/EPS_D**2) UFRIC2 = SQRT(TAUW(IP) + TAUHF(IP)) ! New friction velocity TAUTOT(IP) = TAUW(IP) + TAUHF(IP) ! Check for convergence ! WRITE(DBG%FHNDL,*) 'ITERATION =', I ! WRITE(DBG%FHNDL,*) 'wind10 =', wind10 ! WRITE(DBG%FHNDL,*) 'fU10 =', fU10 ! WRITE(DBG%FHNDL,*) 'z0_t =', z0_t ! WRITE(DBG%FHNDL,*) 'z0 =', z0(ip) ! WRITE(DBG%FHNDL,*) 'Hs =', HS_w, 'L =', Lur, 'TP =', TP_w ! WRITE(DBG%FHNDL,*) 'Ulur =' ,ULur, 'Ur =', Ur, 'EPS_D =', EPS_D ! WRITE(DBG%FHNDL,*) 'T_ds & T_t =', TAUW(ip), TAUHF(ip) ! WRITE(DBG%FHNDL,*) 'UFRIC =', UFRIC1, UFRIC2 !stop 'wwm_windinput.F90 l.141' IF ( (ABS(UFRIC2-UFRIC1))/UFRIC1 .LT. SMALL) THEN UFRIC(IP) = UFRIC2 UFRIC(IP) = MAX(THR,UFRIC(IP)) EXIT END IF END DO 101 CONTINUE !WRITE(DBG%FHNDL,*) 'ITERATION =', I !WRITE(DBG%FHNDL,*) 'wind10 =', wind10 !WRITE(DBG%FHNDL,*) 'fU10 =', fU10 !WRITE(DBG%FHNDL,*) 'z0_t =', z0_t !WRITE(DBG%FHNDL,*) 'z_0 =', z0(ip) !WRITE(DBG%FHNDL,*) 'Hs =', HS_w, 'L =', Lur, 'TP =', TP_w !WRITE(DBG%FHNDL,*) 'Ulur =' ,ULur, 'Ur =', Ur, 'EPS_D =', EPS_D !WRITE(DBG%FHNDL,*) 'T_ds & T_t =', TAUW(ip), TAUHF(ip) !WRITE(DBG%FHNDL,*) 'UFRIC =', UFRIC1, UFRIC2 FPM = G9 / ( 28.0_rkind * UFRIC(IP) ) CD(IP) = (UFRIC(IP)/WIND10)**2 Z0(IP) = 10.0_rkind/EXP(KAPPA*WIND10 /UFRIC(IP)) TAUTOT(IP)=(UFRIC(IP)**2)*rhoa IF (UFRIC2 .LT. VERYSMALL) THEN ALPHA_CH(IP) = 0._rkind ELSE ALPHA_CH(IP) = g9 * z0(ip) / UFRIC2 ENDIF CASE DEFAULT END SELECT RETURN END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE SIN_LIN_CAV( IP, WIND10, WINDTH, FPM, IMATRA, SSINL ) ! ! Linear growth term according to Cavaleri & Melanotte Rizolli ... ! USE DATAPOOL IMPLICIT NONE INTEGER, INTENT(IN) :: IP REAL(rkind) , INTENT(OUT) :: IMATRA(MSC,MDC) REAL(rkind) , INTENT(OUT) :: SSINL(MSC,MDC) REAL(rkind) , INTENT(IN) :: WINDTH, WIND10 REAL(rkind) , INTENT(IN) :: FPM INTEGER :: IS, ID REAL(rkind) :: AUX1, AUX2, AUX3, AUXH, DIFFWND REAL(rkind) :: SWINA, TOLFIL, COSWIND, SINWIND AUX1 = 0.0015 / ( PI2*G9**2 ) DO IS = 1, MSC TOLFIL = EXP(-(MIN(2., FPM / SPSIG(IS))**4)) AUX2 = AUX1 / SPSIG(IS) DO ID = 1, MDC IF (SPSIG(IS) .GE. (0.7 * FPM) ) THEN AUX3 = ( WIND10/28. * MAX( 0. , (COSTH(ID)*COS(WINDTH) + SINTH(ID)*SIN(WINDTH))))**4 IMATRA(IS,ID) = MAX( 0. , AUX2 * AUX3 * TOLFIL ) !WRITE(*,'(5F20.10)') AUX1, AUX2, AUX3, DIFFWND, TOLFIL, IMATRA(IS,ID) SSINL(IS,ID) = IMATRA(IS,ID) ENDIF ENDDO ENDDO END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE SIN_EXP_KOMEN( IP, WINDTH, ACLOC, IMATRA, IMATDA, SSINE ) USE DATAPOOL IMPLICIT NONE ! ! *** the exponential growth term by Komen et al. (1984) *** ! INTEGER, INTENT(IN) :: IP REAL(rkind) , INTENT(IN) :: WINDTH REAL(rkind) , INTENT(IN) :: ACLOC(MSC,MDC) REAL(rkind) , INTENT(OUT) :: SSINE(MSC,MDC) REAL(rkind) , INTENT(INOUT):: IMATRA(MSC,MDC), IMATDA(MSC,MDC) INTEGER :: IS, ID REAL(rkind) :: AUX1, AUX2, AUX3 REAL(rkind) :: SWINB, CINV, COSDIF, SFIE(MSC,MDC) AUX1 = 0.25_rkind * RHOAW AUX2 = 28._rkind * UFRIC(IP) DO IS = 1, MSC CINV = WK(IS,IP)/SPSIG(IS) AUX3 = AUX2 * CINV DO ID = 1, MDC COSDIF = MyCOS(SPDIR(ID)-WINDTH) SWINB = AUX1 * ( AUX3 * COSDIF - 1.0_rkind ) SWINB = MAX( 0.0_rkind, SWINB * SPSIG(IS) ) SSINE(IS,ID) = SWINB * ACLOC(IS,ID) !WRITE(DBG%FHNDL,'(2I10,4F15.8)') IS, ID, SSINE(IS,ID), AUX3, AUX2, AUX1 IMATRA(IS,ID) = IMATRA(IS,ID) + SSINE(IS,ID) END DO END DO RETURN END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE SIN_MAKIN(IP, WIND10, WINDTH, KMESPC, ETOT, ACLOC, IMATRA, IMATDA, SSINE) USE DATAPOOL IMPLICIT NONE INTEGER, INTENT(IN) :: IP REAL(rkind) , INTENT(IN) :: WIND10, WINDTH REAL(rkind) , INTENT(IN) :: ACLOC(MSC,MDC) REAL(rkind) , INTENT(OUT) :: SSINE(MSC,MDC) REAL(rkind) , INTENT(INOUT) :: IMATRA(MSC,MDC), IMATDA(MSC,MDC) INTEGER :: IS, ID REAL(rkind) :: AUX1, AUX2 REAL(rkind) :: SWINB, CINV, SIGMA REAL(rkind) :: COSWW, THETA REAL(rkind) :: NC_MK, MC_MK, MBETA, RMK REAL(rkind) :: KMESPC, ETOT, DS, ELOCAL REAL(rkind) :: STEEPLOCAL REAL(rkind) :: ALOCAL, CPHASE, CYS, ATTC LOGICAL :: LATT, LOPP ! ! PARAMETER FROM MAKIN 1999 ! MC_MK = 0.3_rkind NC_MK = 5.0_rkind LOPP = .FALSE. CYS = -25._rkind ! Opposing wind attenuation. LATT = .FALSE. ATTC = -10.0_rkind ! Attenuation coefficient MBETA = 32.0_rkind ! See orignial Paper A GU OCEANS VOL. 104, No.: C4, April 1999 and see Makin & Stam 2003 (KNMI) DO IS = 1, MSC CINV = WK(IS,IP) / SPSIG(IS) SIGMA = SPSIG(IS) IF (WIND10 .LE. THR) THEN AUX1 = 0.0_rkind ELSE AUX1 = 1._rkind/CINV/WIND10 END IF AUX2 = UFRIC(IP)*CINV RMK = 1 - MC_MK * AUX1 ** NC_MK DO ID = 1, MDC THETA = SPDIR(ID) COSWW = MyCOS(THETA-WINDTH) IF (LATT) THEN IF (RMK .GE. 0.0_rkind) THEN SWINB = MBETA * RMK * RHOAW * AUX2**2 * COSWW * ABS(COSWW) * SIGMA END IF IF (COSWW * ABS(COSWW) .GE. 0.0_rkind .AND. RMK .LT. 0.0_rkind) THEN SWINB = MAX(ATTC,MBETA*RMK) * RHOAW * AUX2**2 * COSWW * ABS(COSWW) * SIGMA END IF ELSE IF (RMK .GT. ZERO) THEN SWINB = MBETA * RMK * RHOAW * AUX2**2 * COSWW * ABS(COSWW) * SIGMA ELSE SWINB = ZERO END IF END IF IF (COSWW * ABS(COSWW) .LE. ZERO) THEN IF (LOPP) THEN CPHASE = 0.0_rkind/CINV IF (IS .EQ. 1) DS = SPSIG(IS) IF (IS .GT. 1) DS = SPSIG(IS) - SPSIG(IS-1) IF (IS .EQ. 1) ELOCAL = ONEHALF * ACLOC(IS,ID) * SPSIG(IS) * SPSIG(IS) ! Simpson IF (IS .GT. 1) ELOCAL = ONEHALF * ( ACLOC(IS,ID) * SPSIG(IS) + ACLOC(IS-1,ID) * SPSIG(IS-1) ) * DS ALOCAL = SQRT(8.0_rkind*ELOCAL) STEEPLOCAL = ALOCAL * WK(IS,IP) SWINB = CYS * RHOAW * STEEPLOCAL * STEEPLOCAL * (ONE - ((WIND10 * COSWW)/CPHASE) ) **2 * SIGMA ELSE SWINB = ZERO END IF END IF SSINE(IS,ID) = SWINB * ACLOC(IS,ID) IF (ICOMP .GE. 2) THEN IF (SWINB .LT. 0) THEN IMATDA(IS,ID) = - SWINB ELSE IMATRA(IS,ID) = IMATRA(IS,ID) + SSINE(IS,ID) END IF ELSE IF (ICOMP .LT. 2) THEN IF (SWINB .LT. 0) THEN IMATDA(IS,ID) = SWINB IMATRA(IS,ID) = IMATRA(IS,ID) + SSINE(IS,ID) ELSE IMATRA(IS,ID) = IMATRA(IS,ID) + SSINE(IS,ID) END IF END IF END DO END DO RETURN END SUBROUTINE !********************************************************************** !* * !**********************************************************************