#include "wwm_functions.h" !* * !********************************************************************** SUBROUTINE COMPUTE_DIRECTION_CNTG_A ! ! *** Crank - Nicolson Method ( Implicit ) ! USE DATAPOOL IMPLICIT NONE LOGICAL :: ISEQ0 REAL(rkind) :: GF(MDC) REAL(rkind) :: AMAT(3,MDC) REAL(rkind) :: TMPAC(MDC) REAL(rkind) :: CAD(MSC,MDC) INTEGER :: ID1, ID2 INTEGER :: IP, IS, ID, IT, ISTEP REAL(rkind) :: AUX1, AUX2, REST, DT4DI REAL(rkind) :: AM1 ! Lower Left Element REAL(rkind) :: A1M ! Upper Right Element REAL(rkind) :: CFLCAD TMPAC(:) = 0.0 DO IP = 1, MNP IF ((ABS(IOBP(IP)) .EQ. 1 .OR. ABS(IOBP(IP)) .EQ. 3) .AND. .NOT. LTHBOUND) CYCLE IF (DEP(IP) .LT. DMIN) CYCLE IF (IOBP(IP) .EQ. 2) CYCLE CALL PROPTHETA(IP,CAD) DO IS = 1, MSC GF(:) = 0.0 AMAT(:,:) = 0.0 TMPAC(:) = AC2(IS,:,IP) CFLCAD = MAXVAL(ABS(CAD(IS,:))/DDIR * DT4D) IF (ISEQ0(CFLCAD)) CYCLE REST = ABS(MOD(CFLCAD,1.0_rkind)) IF (REST .GT. THR .AND. REST .LT. 0.5) THEN ISTEP = ABS(NINT(CFLCAD)) + 1 ELSE ISTEP = ABS(NINT(CFLCAD)) END IF DT4DI = DT4D/MyREAL(ISTEP) DO IT = 1, ISTEP DO ID = 1, MDC ! Start Iteration ID1 = ID - 1 ID2 = ID + 1 IF (ID == 1) THEN ID1 = MDC AUX1 = (DT4DI/DDIR/2.0*CAD(IS,ID1)) A1M = -1.0*AUX1*RTHETA ELSE AUX1 = (DT4DI/DDIR/2.0*CAD(IS,ID1)) AMAT(1,ID) = -1.0*AUX1*RTHETA END IF AMAT(2,ID) = 1.0 IF (ID == MDC) THEN ID2 = 1 AUX2 = (DT4DI/DDIR/2.0*CAD(IS,ID2)) AM1 = AUX2*RTHETA ELSE AUX2 = (DT4DI/DDIR/2.0*CAD(IS,ID2)) AMAT(3,ID) = AUX2*RTHETA END IF GF(ID) = TMPAC(ID)-(1.0-RTHETA)*(AUX2*TMPAC(ID2)-AUX1*TMPAC(ID1) ) END DO CALL GAUS1D (MDC, AMAT, GF, AM1, A1M) TMPAC(:) = GF(:) END DO ! End Iteration AC2(IS,:,IP) = TMPAC(:) END DO END DO END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE COMPUTE_DIRECTION_QUICKEST_A USE DATAPOOL IMPLICIT NONE INTEGER :: IP, IS, ID INTEGER :: IT, ITER REAL(rkind) :: CAD(MSC,MDC) REAL(rkind) :: ACQ(0:MDC+1) REAL(rkind) :: CADS(0:MDC+1) REAL(rkind) :: DT4DI REAL(rkind) :: REST REAL(rkind) :: CFLCAD !$OMP PARALLEL DEFAULT(NONE) & !$OMP& SHARED(AC2,SI,DAC_THE,DT4D,DDIR,LCIRD,LTHBOUND,MNP, & !$OMP& MSC,MDC,WK,DEP,DMIN,CURTXY,ICOMP,IOBP,DAC_ADV,DAC_SIG,DAC_SOU) & !$OMP& PRIVATE(IP,IS,DT4DI,ITER,CFLCAD,REST,CADS,CAD,ACQ) !$OMP DO DO IP = 1, MNP ! IF (LQSTEA .AND. IP_IS_STEADY(IP) .EQ. 1) CYCLE IF ((ABS(IOBP(IP)) .EQ. 1 .OR. IOBP(IP) .EQ. 3) .AND. .NOT. LTHBOUND) CYCLE ! skip boudary points if set so ... IF (DEP(IP) .LT. DMIN) CYCLE ! skip dry nodes ... IF (IOBP(IP) .EQ. 2) CYCLE ! skip active boundary points ... CALL PROPTHETA(IP,CAD) DO IS = 1, MSC ACQ(1:MDC) = AC2(IS,:,IP) CADS(1:MDC) = CAD(IS,:) CFLCAD = MAXVAL(ABS(CAD(IS,:)))*DT4D/DDIR IF (CFLCAD .LT. THR) CYCLE REST = ABS(MOD(CFLCAD,1.0_rkind)) IF (REST .GT. THR .AND. REST .LT. 0.5) THEN ITER = ABS(NINT(CFLCAD)) + 1 ELSE ITER = ABS(NINT(CFLCAD)) END IF ITER = MAX(1,ITER) DT4DI = DT4D / MyREAL(ITER) DO IT = 1, ITER ! Iteration CALL QUICKEST_DIR(MDC,LCIRD,ACQ,CADS,DT4DI,DDIR) END DO ! end Interation AC2(IS,:,IP) = MAX(ZERO,ACQ(1:MDC)) END DO END DO !$OMP END DO NOWAIT !$OMP END PARALLEL END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE COMPUTE_DIRECTION_UPWIND_A USE DATAPOOL IMPLICIT NONE INTEGER :: IP, IS, ID, IT, ISTEP, ID1, ID2 REAL(rkind) :: CAD(MSC,MDC) REAL(rkind) :: TMP(MDC), REST, CFLCAD, DT4DI REAL(rkind) :: LP(MDC), LM(MDC), CP(MDC), CM(MDC), U0(MDC) DO IP = 1, MNP IF ((ABS(IOBP(IP)) .EQ. 1 .OR. ABS(IOBP(IP)) .EQ. 3) .AND. .NOT. LTHBOUND) CYCLE IF (DEP(IP) .LT. DMIN) CYCLE IF (IOBP(IP) .EQ. 2) CYCLE CALL PROPTHETA(IP,CAD) DO IS = 1, MSC U0(:) = AC2(IS,:,IP) CP(:) = MAX(ZERO,CAD(IS,:)) CM(:) = MIN(ZERO,CAD(IS,:)) CFLCAD = MAXVAL ( ABS(CAD(IS,:))*DT4D/DDIR ) REST = ABS(MOD(CFLCAD,1._rkind)) IF (CFLCAD .LT. THR) CYCLE REST = ABS(MOD(CFLCAD,1._rkind)) IF (REST .GT. THR .AND. REST .LT. 0.5_rkind) THEN ISTEP = ABS(NINT(CFLCAD)) + 1 ELSE ISTEP = ABS(NINT(CFLCAD)) END IF DT4DI = DT4D / MyREAL(ISTEP) DO IT = 1, ISTEP DO ID = 1, MDC ID1 = ID - 1 ID2 = ID + 1 IF (ID .EQ. 1) ID1 = MDC IF (ID .EQ. MDC) ID2 = 1 LP(ID) = - 1./DDIR * ( CP(ID1) * U0(ID1) - CP(ID) * U0(ID) ) LM(ID) = - 1./DDIR * ( CM(ID2) * U0(ID2) - CM(ID) * U0(ID) ) TMP(ID) = U0(ID) - DT4DI * ( -LM(ID) + LP(ID) ) END DO TMP(ID) = TMP(1) ! That lines looks dangerous memorywise ! and useless. U0(:) = TMP(:) END DO AC2(IS,:,IP) = MAX(0._rkind,U0(:)) END DO END DO END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE COMPUTE_DIRECTION_UPWIND_IMPLICIT USE DATAPOOL IMPLICIT NONE INTEGER :: IP, IS, ID, IT, ISTEP, ID1, ID2 REAL(rkind) :: CAD(MSC,MDC) REAL(rkind) :: TMP(MDC) REAL(rkind) :: LP(MDC), LM(MDC), CP(MDC), CM(MDC), U0(MDC) REAL(rkind) :: EMAT(MDC,MDC) DO IP = 1, MNP IF ((ABS(IOBP(IP)) .EQ. 1 .OR. ABS(IOBP(IP)) .EQ. 3) .AND. .NOT. LTHBOUND) CYCLE IF (DEP(IP) .LT. DMIN) CYCLE IF (IOBP(IP) .EQ. 2) CYCLE CALL PROPTHETA(IP,CAD) DO IS = 1, MSC U0 = AC2(IS,:,IP) CP = MAX(ZERO,CAD(IS,:)) CM = MIN(ZERO,CAD(IS,:)) EMAT=ZERO DO ID=1,MDC ID1 = ID - 1 ID2 = ID + 1 IF (ID .EQ. 1) ID1 = MDC IF (ID .EQ. MDC) ID2 = 1 EMAT(ID,ID ) = 1 + (DT4D/DDIR) * (CP(ID) - CM(ID)) EMAT(ID,ID1) = - (DT4D/DDIR) * CP(ID1) EMAT(ID,ID2) = (DT4D/DDIR) * CM(ID2) END DO CALL GAUSS_SOLVER(MDC, EMAT, TMP, U0) AC2(IS,:,IP) = MAX(0._rkind,TMP) END DO END DO END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE COMPUTE_DIRECTION_WENO_A USE DATAPOOL IMPLICIT NONE INTEGER :: ISTEP INTEGER :: IP, IS, ID, IT INTEGER :: ID1, ID11, ID2, ID22, ID23 REAL(rkind) :: L(MDC), FLP(MDC) REAL(rkind) :: FLM(MDC), NFLP(MDC), NFLM(MDC) REAL(rkind) :: CP(MDC), CM(MDC) REAL(rkind) :: CAD(MSC,MDC) REAL(rkind) :: WP1, WP2, WP3, WM1, WM2, WM3 REAL(rkind) :: WI_P1, WI_P2, WI_P3, WI_M1, WI_M2, WI_M3 REAL(rkind) :: FO_FLP1, FO_FLP2, FO_FLP3, FO_FLM1, FO_FLM2, FO_FLM3, EPS, REST REAL(rkind) :: TMP, TMP1, TMP2, TMP3 REAL(rkind) :: BETAM1, BETAM2,BETAM3, BETAP1, BETAP2,BETAP3 REAL(rkind) :: U0(MDC), U1(MDC), U2(MDC), U3(MDC) REAL(rkind) :: FLPID, FLPID1, FLPID11, FLPID2, FLPID22 REAL(rkind) :: FLMID, FLMID1, FLMID2, FLMID22, FLMID23 REAL(rkind) :: AP11, AP12, AP13, AP21, AP22, AP23, AP31, AP32, AP33 REAL(rkind) :: BP1, BP2, GAMMA1, GAMMA2, GAMMA3, CFLCAD, DT4DI EPS = 10E-10_rkind GAMMA1 = 0.1_rkind GAMMA2 = 0.6_rkind GAMMA3 = 0.3_rkind AP11 = (2.0_rkind)/(6.0_rkind) AP12 = -(7.0_rkind)/(6.0_rkind) AP13 = (11.0_rkind)/(6.0_rkind) AP21 = -(1.0_rkind)/(6.0_rkind) AP22 = (5.0_rkind)/(6.0_rkind) AP23 = (2.0_rkind)/(6.0_rkind) AP31 = (2.0_rkind)/(6.0_rkind) AP32 = (5.0_rkind)/(6.0_rkind) AP33 = -(1.0_rkind)/(6.0_rkind) BP1 = (13.0_rkind)/(12.0_rkind) BP2 = (1.0_rkind)/(4.0_rkind) DO IP = 1, MNP IF ((ABS(IOBP(IP)) .EQ. 1 .OR. ABS(IOBP(IP)) .EQ. 3) .AND. .NOT. LTHBOUND) CYCLE IF (DEP(IP) .LT. DMIN) CYCLE IF (IOBP(IP) .EQ. 2) CYCLE CALL PROPTHETA(IP,CAD) DO IS = 1, MSC U0(:) = AC2(IS,:,IP) CP(:) = MAX(ZERO,CAD(IS,:)) CM(:) = MIN(ZERO,CAD(IS,:)) CFLCAD = MAXVAL(ABS(CAD(IS,:)))*DT4D/DDIR REST = ABS(MOD(CFLCAD,ONE)) ISTEP = ABS(NINT(CFLCAD)) + 1 DT4DI = DT4D / MyREAL(ISTEP) DO IT = 1, ISTEP FLP(:) = CP(:) * U0(:) FLM(:) = CM(:) * U0(:) DO ID = 1, MDC ID1 = ID - 1 ID11 = ID - 2 ID2 = ID + 1 ID22 = ID + 2 ID23 = ID + 3 IF (ID .EQ. 1) THEN ID1 = MDC ID11 = MDC - 1 ELSE IF (ID .EQ. 2) THEN ID11 = MDC ELSE IF (ID .EQ. MDC - 2) THEN ID23 = 1 ELSE IF (ID .EQ. MDC - 1) THEN ID22 = 1 ID23 = 2 ELSE IF (ID .EQ. MDC) THEN ID2 = 1 ID22 = 2 ID23 = 3 END IF FLPID = FLP(ID) FLPID1 = FLP(ID1) FLPID11 = FLP(ID11) FLPID2 = FLP(ID2) FLPID22 = FLP(ID22) FLMID = FLM(ID) FLMID1 = FLM(ID1) FLMID2 = FLM(ID2) FLMID22 = FLM(ID22) FLMID23 = FLM(ID23) FO_FLP1 = AP11 * FLPID11 + AP12 * FLPID1 + AP13 * FLPID FO_FLP2 = AP21 * FLPID1 + AP22 * FLPID + AP23 * FLPID2 FO_FLP3 = AP31 * FLPID + AP32 * FLPID2 + AP33 * FLPID22 FO_FLM1 = AP33 * FLMID1 + AP32 * FLMID + AP31 * FLMID2 FO_FLM2 = AP23 * FLMID + AP22 * FLMID2 + AP21 * FLMID22 FO_FLM3 = AP13 * FLMID2 + AP12 * FLMID22 + AP11 * FLMID23 BETAP1= BP1 * (FLPID11-2.*FLPID1+FLPID )**2+BP2*( FLPID11-4.*FLPID1+3.* FLPID )**2 BETAP2= BP1 * (FLPID1 -2.*FLPID+FLPID2 )**2+BP2*( FLPID1 - FLPID2 )**2 BETAP3= BP1 * (FLPID -2.*FLPID2+FLPID22)**2+BP2*(3.* FLPID -4.*FLPID2+ FLPID22)**2 BETAM1= BP1 * (FLMID2-2.*FLMID22+FLMID23)**2+BP2*( FLMID2-4.*FLMID22+3.* FLMID23)**2 BETAM2= BP1 * (FLMID -2.*FLMID2 +FLMID22)**2+BP2*( FLMID - FLMID22)**2 BETAM3= BP1 * (FLMID1-2.*FLMID +FLMID2 )**2+BP2*(3.* FLMID1-4.*FLMID + FLMID2 )**2 TMP = 1./(EPS + BETAP1) WP1 = GAMMA1 * TMP * TMP TMP = 1./(EPS + BETAP2) WP2 = GAMMA2 * TMP * TMP TMP = 1./(EPS + BETAP3) WP3 = GAMMA3 * TMP * TMP TMP = 1./(EPS + BETAM1) WM1 = GAMMA1 * TMP * TMP TMP = 1./(EPS + BETAM2) WM2 = GAMMA2 * TMP * TMP TMP = 1./(EPS + BETAM3) WM3 = GAMMA3 * TMP * TMP TMP = 1./(WP1+WP2+WP3) WI_P1 = WP1*TMP WI_P2 = WP2*TMP WI_P3 = WP3*TMP TMP = 1./(WM1+WM2+WM3) WI_M1 = WM1*TMP WI_M2 = WM2*TMP WI_M3 = WM3*TMP NFLP(ID) = WI_P1 * FO_FLP1 + WI_P2 * FO_FLP2 + WI_P3 * FO_FLP3 NFLM(ID) = WI_M3 * FO_FLM1 + WI_M2 * FO_FLM2 + WI_M1 * FO_FLM3 END DO DO ID = 1, MDC ID1 = ID - 1 IF (ID .EQ. 1) ID1 = MDC L(ID) = 1./DDIR * ( (NFLP(ID)-NFLP(ID1)) + (NFLM(ID)-NFLM(ID1)) ) END DO DO ID = 1, MDC U1(ID) = U0(ID) - DT4DI * L(ID) END DO FLP(:) = CP(:) * U1(:) FLM(:) = CM(:) * U1(:) DO ID = 1, MDC ID1 = ID - 1 ID11 = ID - 2 ID2 = ID + 1 ID22 = ID + 2 ID23 = ID + 3 IF (ID .EQ. 1) THEN ID1 = MDC ID11 = MDC - 1 ELSE IF (ID .EQ. 2) THEN ID11 = MDC ELSE IF (ID .EQ. MDC - 2) THEN ID23 = 1 ELSE IF (ID .EQ. MDC - 1) THEN ID22 = 1 ID23 = 2 ELSE IF (ID .EQ. MDC) THEN ID2 = 1 ID22 = 2 ID23 = 3 END IF FLPID = FLP(ID) FLPID1 = FLP(ID1) FLPID11 = FLP(ID11) FLPID2 = FLP(ID2) FLPID22 = FLP(ID22) FLMID = FLM(ID) FLMID1 = FLM(ID1) FLMID2 = FLM(ID2) FLMID22 = FLM(ID22) FLMID23 = FLM(ID23) FO_FLP1 = AP11 * FLPID11 + AP12 * FLPID1 + AP13 * FLPID FO_FLP2 = AP21 * FLPID1 + AP22 * FLPID + AP23 * FLPID2 FO_FLP3 = AP31 * FLPID + AP32 * FLPID2 + AP33 * FLPID22 FO_FLM1 = AP33 * FLMID1 + AP32 * FLMID + AP31 * FLMID2 FO_FLM2 = AP23 * FLMID + AP22 * FLMID2 + AP21 * FLMID22 FO_FLM3 = AP13 * FLMID2 + AP12 * FLMID22 + AP11 * FLMID23 BETAP1= BP1 * (FLPID11-2*FLPID1+FLPID )**2+BP2*( FLPID11-4.*FLPID1+3.* FLPID )**2 BETAP2= BP1 * (FLPID1 -2*FLPID+FLPID2 )**2+BP2*( FLPID1 - FLPID2 )**2 BETAP3= BP1 * (FLPID -2*FLPID2+FLPID22)**2+BP2*(3.* FLPID -4.*FLPID2+ FLPID22)**2 BETAM1= BP1 * (FLMID2-2*FLMID22+FLMID23)**2+BP2*( FLMID2-4.*FLMID22+3.* FLMID23)**2 BETAM2= BP1 * (FLMID -2*FLMID2 +FLMID22)**2+BP2*( FLMID - FLMID22)**2 BETAM3= BP1 * (FLMID1-2*FLMID +FLMID2 )**2+BP2*(3.* FLMID1-4.*FLMID + FLMID2 )**2 TMP1 = 1./(EPS + BETAP1) WP1 = GAMMA1 * TMP1 * TMP1 TMP2 = 1./(EPS + BETAP2) WP2 = GAMMA2 * TMP2 * TMP2 TMP3 = 1./(EPS + BETAP3) WP3 = GAMMA3 * TMP3 * TMP3 TMP1 = 1./(EPS + BETAM1) WM1 = GAMMA1 * TMP1 * TMP1 TMP2 = 1./(EPS + BETAM2) WM2 = GAMMA2 * TMP2 * TMP2 TMP3 = 1./(EPS + BETAM3) WM3 = GAMMA3 * TMP3 * TMP3 TMP = 1./(WP1+WP2+WP3) WI_P1 = WP1*TMP WI_P2 = WP2*TMP WI_P3 = WP3*TMP TMP = 1./(WM1+WM2+WM3) WI_M1 = WM1*TMP WI_M2 = WM2*TMP WI_M3 = WM3*TMP NFLP(ID) = WI_P1 * FO_FLP1 + WI_P2 * FO_FLP2 + WI_P3 * FO_FLP3 NFLM(ID) = WI_M3 * FO_FLM1 + WI_M2 * FO_FLM2 + WI_M1 * FO_FLM3 END DO DO ID = 1, MDC ID1 = ID - 1 IF (ID .EQ. 1) ID1 = MDC L(ID) = 1./DDIR * ( (NFLP(ID)-NFLP(ID1)) + (NFLM(ID)-NFLM(ID1)) ) END DO DO ID = 1, MDC U2(ID) = 3./4. * U0(ID) + 1./4. * U1(ID) - 1./4. * DT4DI * L(ID) END DO FLP(:) = CP(:) * U2(:) FLM(:) = CM(:) * U2(:) DO ID = 1, MDC ID1 = ID - 1 ID11 = ID - 2 ID2 = ID + 1 ID22 = ID + 2 ID23 = ID + 3 IF (ID .EQ. 1) THEN ID1 = MDC ID11 = MDC - 1 ELSE IF (ID .EQ. 2) THEN ID11 = MDC ELSE IF (ID .EQ. MDC - 2) THEN ID23 = 1 ELSE IF (ID .EQ. MDC - 1) THEN ID22 = 1 ID23 = 2 ELSE IF (ID .EQ. MDC) THEN ID2 = 1 ID22 = 2 ID23 = 3 END IF FLPID = FLP(ID) FLPID1 = FLP(ID1) FLPID11 = FLP(ID11) FLPID2 = FLP(ID2) FLPID22 = FLP(ID22) FLMID = FLM(ID) FLMID1 = FLM(ID1) FLMID2 = FLM(ID2) FLMID22 = FLM(ID22) FLMID23 = FLM(ID23) FO_FLP1 = AP11 * FLPID11 + AP12 * FLPID1 + AP13 * FLPID FO_FLP2 = AP21 * FLPID1 + AP22 * FLPID + AP23 * FLPID2 FO_FLP3 = AP31 * FLPID + AP32 * FLPID2 + AP33 * FLPID22 FO_FLM1 = AP33 * FLMID1 + AP32 * FLMID + AP31 * FLMID2 FO_FLM2 = AP23 * FLMID + AP22 * FLMID2 + AP21 * FLMID22 FO_FLM3 = AP13 * FLMID2 + AP12 * FLMID22 + AP11 * FLMID23 BETAP1= BP1 * (FLPID11-2*FLPID1+FLPID )**2+BP2*( FLPID11-4.*FLPID1+3.* FLPID )**2 BETAP2= BP1 * (FLPID1 -2*FLPID+FLPID2 )**2+BP2*( FLPID1 - FLPID2 )**2 BETAP3= BP1 * (FLPID -2*FLPID2+FLPID22)**2+BP2*(3.* FLPID -4.*FLPID2+ FLPID22)**2 BETAM1= BP1 * (FLMID2-2*FLMID22+FLMID23)**2+BP2*( FLMID2-4.*FLMID22+3.* FLMID23)**2 BETAM2= BP1 * (FLMID -2*FLMID2 +FLMID22)**2+BP2*( FLMID - FLMID22)**2 BETAM3= BP1 * (FLMID1-2*FLMID +FLMID2 )**2+BP2*(3.* FLMID1-4.*FLMID + FLMID2 )**2 TMP = 1./(EPS + BETAP1) WP1 = GAMMA1 * TMP * TMP TMP = 1./(EPS + BETAP2) WP2 = GAMMA2 * TMP * TMP TMP = 1./(EPS + BETAP3) WP3 = GAMMA3 * TMP * TMP TMP = 1./(EPS + BETAM1) WM1 = GAMMA1 * TMP * TMP TMP = 1./(EPS + BETAM2) WM2 = GAMMA2 * TMP * TMP TMP = 1./(EPS + BETAM3) WM3 = GAMMA3 * TMP * TMP TMP = 1./(WP1+WP2+WP3) WI_P1 = WP1*TMP WI_P2 = WP2*TMP WI_P3 = WP3*TMP TMP = 1./(WM1+WM2+WM3) WI_M1 = WM1*TMP WI_M2 = WM2*TMP WI_M3 = WM3*TMP NFLP(ID) = WI_P1 * FO_FLP1 + WI_P2 * FO_FLP2 + WI_P3 * FO_FLP3 NFLM(ID) = WI_M3 * FO_FLM1 + WI_M2 * FO_FLM2 + WI_M1 * FO_FLM3 END DO DO ID = 1, MDC ID1 = ID - 1 IF (ID .EQ. 1) ID1 = MDC L(ID) = 1./DDIR * ( (NFLP(ID)-NFLP(ID1)) + (NFLM(ID)-NFLM(ID1)) ) END DO DO ID = 1, MDC U3(ID) = 1./3. * U0(ID) + 2./3. * U2(ID) - 2./3. * DT4DI * L(ID) END DO U0(:) = U3(:) END DO AC2(IS,:,IP) = MAX(ZERO,MyREAL(U3(:))) END DO END DO END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE PROPTHETA(IP,CAD) USE DATAPOOL IMPLICIT NONE INTEGER, INTENT(IN) :: IP REAL(rkind), INTENT(OUT) :: CAD(MSC,MDC) INTEGER :: IS, ID REAL(rkind) :: WKDEP, DWDH, CFL CAD = 0. SELECT CASE (DIMMODE) CASE (1) IF (DEP(IP) .GT. DMIN) THEN ! ar: obsolete check since done in calling routine ... DO IS = 1, MSC WKDEP = WK(IS,IP) * DEP(IP) IF (WKDEP .LT. 13.) THEN DWDH = SPSIG(IS)/SINH(MIN(KDMAX,TWO*WKDEP)) DO ID = 1, MDC CAD(IS,ID) = DWDH*(SINTH(ID)*DDEP(IP,1)) END DO END IF END DO IF (LSTCU .OR. LSECU) THEN DO IS = 1, MSC DO ID = 1, MDC CAD(IS,ID) = CAD(IS,ID) + SIN2TH(ID)*DCUY(IP,1) + COSTH(ID)*DCUX(IP,1) END DO END DO END IF END IF CASE (2) IF (DEP(IP) .GT. DMIN) THEN DO IS = 1, MSC WKDEP = WK(IS,IP) * DEP(IP) IF (WKDEP .LT. 13.) THEN DWDH = SPSIG(IS)/SINH(MIN(KDMAX,2.*WKDEP)) !AR: Joseph: There is a small very tiny different in terms of CPU !when we do refraction ... it has something to do with the depth !gradients ... DO ID = 1, MDC CAD(IS,ID) = DWDH * ( SINTH(ID)*DDEP(IP,1)-COSTH(ID)*DDEP(IP,2) ) END DO ENDIF END DO IF (LSTCU .OR. LSECU) THEN DO IS = 1, MSC DO ID = 1, MDC CAD(IS,ID) = CAD(IS,ID) + SIN2TH(ID)*DCUY(IP,1)-COS2TH(ID)*DCUX(IP,2)+SINCOSTH(ID)*( DCUX(IP,1)-DCUY(IP,2) ) END DO END DO END IF IF (LDIFR) THEN DO IS = 1, MSC DO ID = 1, MDC CAD(IS,ID) = DIFRM(IP)*CAD(IS,ID)-CG(IS,IP)*(DIFRX(IP)*SINTH(ID)-DIFRY(IP)*COSTH(ID)) END DO END DO END IF ELSE CAD = 0. END IF CASE DEFAULT CALL WWM_ABORT('PROPTHETA') END SELECT IF (LSPHE) THEN DO IS = 1, MSC DO ID = 1, MDC CAD(IS,ID) = CAD(IS,ID)-CG(IS,IP)*COSTH(ID)*TAN(YP(IP)*DEGRAD)/REARTH !WRITE(*,*) IS, ID, CG(IS,IP), COSTH(ID), TAN(YP(IP)*DEGRAD)/REARTH END DO END DO END IF IF (LFILTERTH) THEN DO IS = 1, MSC CFL = MAXVAL(ABS(CAD(IS,:)))*DT4D/DDIR IF (CFL .GT. MAXCFLTH) THEN DO ID = 1, MDC CAD(IS,ID) = MAXCFLTH/CFL*CAD(IS,ID) END DO END IF END DO END IF END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE COMPUTE_DIRECTION USE DATAPOOL IMPLICIT NONE IF (WRITESTATFLAG == 1) THEN WRITE(STAT%FHNDL,'("+TRACE...",A)') 'ENTERING COMPUTE_DIRECTION' FLUSH(STAT%FHNDL) END IF IF (DMETHOD > 0) THEN IF (DMETHOD == 1) THEN CALL COMPUTE_DIRECTION_CNTG_A ELSE IF (DMETHOD == 2) THEN CALL COMPUTE_DIRECTION_QUICKEST_A ELSE IF (DMETHOD == 3) THEN CALL COMPUTE_DIRECTION_WENO_A ELSE IF (DMETHOD == 4) THEN CALL COMPUTE_DIRECTION_UPWIND_A ELSE IF (DMETHOD == 5) THEN CALL COMPUTE_DIRECTION_UPWIND_IMPLICIT END IF END IF IF (WRITESTATFLAG == 1) THEN WRITE(STAT%FHNDL,'("+TRACE...",A)') 'FINISHED COMPUTE_DIRECTION' FLUSH(STAT%FHNDL) END IF IF ( DMETHOD == 1) CALL RESCALE_SPECTRUM END SUBROUTINE COMPUTE_DIRECTION !********************************************************************** !* * !**********************************************************************