#include "wwm_functions.h" #undef positivity #undef DEBUG_COHERENCY_FLUCT !********************************************************************** !* * !********************************************************************** ! for ICOMP == 0 Fully Explicit SUBROUTINE FLUCT_EXPLICIT USE DATAPOOL IMPLICIT NONE INTEGER :: IS, ID, IP REAL(rkind) :: DTMAX !$OMP PARALLEL IF (AMETHOD == 1) THEN !$OMP DO PRIVATE (ID,IS) DO ID = 1, MDC DO IS = 1, MSC CALL EXPLICIT_N_SCHEME(IS,ID) END DO END DO ELSE IF (AMETHOD == 2) THEN !$OMP DO PRIVATE (ID,IS) DO ID = 1, MDC DO IS = 1, MSC CALL EXPLICIT_PSI_SCHEME(IS,ID) END DO END DO ELSE IF (AMETHOD == 3) THEN !$OMP DO PRIVATE (ID,IS) DO ID = 1, MDC DO IS = 1, MSC CALL EXPLICIT_LFPSI_SCHEME(IS,ID) END DO END DO END IF !$OMP END PARALLEL END SUBROUTINE !********************************************************************** !* * !********************************************************************** ! for ICOMP == 1 SUBROUTINE FLUCT_IMP_EXP_SOURCES USE DATAPOOL #ifdef PETSC use PETSC_CONTROLLER, only : EIMPS_PETSC use petsc_block, only: EIMPS_PETSC_BLOCK #endif IMPLICIT NONE INTEGER :: IS, ID, IP REAL(rkind) :: DTMAX #ifdef PETSC ! petsc block has its own loop over MSC MDC IF(AMETHOD == 5) THEN call EIMPS_PETSC_BLOCK RETURN END IF #endif !$OMP PARALLEL IF (AMETHOD == 1) THEN !$OMP DO PRIVATE (ID,IS) DO ID = 1, MDC DO IS = 1, MSC CALL EIMPS_V1( IS, ID) END DO END DO ELSE IF (AMETHOD == 2) THEN !$OMP DO PRIVATE (ID,IS) DO ID = 1, MDC DO IS = 1, MSC CALL CNIMPS( IS, ID) END DO END DO ELSE IF (AMETHOD == 3) THEN !$OMP DO PRIVATE (ID,IS) DO ID = 1, MDC DO IS = 1, MSC CALL CNEIMPS( IS, ID, DTMAX) END DO END DO ELSE IF (AMETHOD == 4) THEN #ifdef PETSC !$OMP DO PRIVATE (ID,IS) DO ID = 1, MDC DO IS = 1, MSC CALL EIMPS_PETSC(IS, ID) END DO END DO #endif ELSE IF (AMETHOD == 6) THEN #ifdef WWM_SOLVER # ifdef MPI_PARALL_GRID CALL WWM_SOLVER_EIMPS(MainLocalColor, SolDat) # endif #endif END IF !$OMP END PARALLEL END SUBROUTINE !********************************************************************** !* * !********************************************************************** ! for ICOMP == 2 SUBROUTINE FLUCT_IMP_SOURCES USE DATAPOOL #ifdef PETSC use PETSC_CONTROLLER, only : EIMPS_PETSC use petsc_block, only: EIMPS_PETSC_BLOCK #endif IMPLICIT NONE INTEGER :: IS, ID, IP REAL(rkind) :: DTMAX !2DO MATHIEU: Please clean this ... and please check this #ifdef PETSC ! petsc block has its own loop over MSC MDC IF(AMETHOD == 5) THEN call EIMPS_PETSC_BLOCK RETURN ENDIF #endif #ifdef WWM_SOLVER # ifdef MPI_PARALL_GRID IF (AMETHOD == 6) THEN CALL WWM_SOLVER_EIMPS(MainLocalColor, SolDat) RETURN ELSE IF (AMETHOD == 7) THEN CALL EIMPS_TOTAL_JACOBI_ITERATION RETURN ENDIF # endif #endif ! !$OMP PARALLEL IF (AMETHOD == 1) THEN !$OMP DO PRIVATE (ID,IS) DO ID = 1, MDC DO IS = 1, MSC CALL EIMPS_V1( IS, ID) END DO END DO ELSE IF (AMETHOD == 2) THEN !$OMP DO PRIVATE (ID,IS) DO ID = 1, MDC DO IS = 1, MSC CALL CNIMPS( IS, ID) END DO END DO ELSE IF (AMETHOD == 3) THEN !$OMP DO PRIVATE (ID,IS) DO ID = 1, MDC DO IS = 1, MSC CALL CNEIMPS( IS, ID, DTMAX) END DO END DO ELSE IF (AMETHOD == 4) THEN #ifdef PETSC !$OMP DO PRIVATE (ID,IS) DO ID = 1, MDC DO IS = 1, MSC CALL EIMPS_PETSC(IS, ID) END DO END DO #endif ENDIF !$OMP END PARALLEL END SUBROUTINE !********************************************************************** !* * !********************************************************************** ! for ICOMP == 3 SUBROUTINE FLUCT_IMP_ALL USE DATAPOOL #ifdef PETSC use PETSC_CONTROLLER, only : EIMPS_PETSC use petsc_block, only: EIMPS_PETSC_BLOCK #endif IMPLICIT NONE INTEGER :: IS, ID, IP REAL(rkind) :: DTMAX IF (AMETHOD .eq.5) THEN #ifdef PETSC CALL EIMPS_PETSC_BLOCK #endif ELSE IF (AMETHOD .eq. 7) THEN #ifdef WWM_SOLVER CALL EIMPS_TOTAL_JACOBI_ITERATION #endif END IF END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE FLUCTCFL(IS, ID, DTMAX) USE DATAPOOL IMPLICIT NONE REAL(rkind) :: K(3,MNE) REAL(rkind), INTENT(OUT) :: DTMAX INTEGER :: IS, ID INTEGER :: I, J, I1, I2, I3 INTEGER :: IP, IE, POS REAL(rkind) :: KSUM, KMAX, LAMBDA(2) REAL(rkind) :: DTMAX_EXP, DTMAX_GLOBAL_EXP REAL(rkind) :: REST, C(2,MNP) DTMAX_GLOBAL_EXP = 10.D14 CALL CADVXY(IS,ID,C) DO IE = 1, MNE I1 = INE(1,IE) I2 = INE(2,IE) I3 = INE(3,IE) LAMBDA(1) = ONESIXTH * (C(1,I1)+C(1,I2)+C(1,I3)) LAMBDA(2) = ONESIXTH * (C(2,I1)+C(2,I2)+C(2,I3)) K(1,IE) = LAMBDA(1) * IEN(1,IE) + LAMBDA(2) * IEN(2,IE) K(2,IE) = LAMBDA(1) * IEN(3,IE) + LAMBDA(2) * IEN(4,IE) K(3,IE) = LAMBDA(1) * IEN(5,IE) + LAMBDA(2) * IEN(6,IE) END DO J = 0 DO IP = 1, MNP KSUM = ZERO KMAX = ZERO DO I = 1, CCON(IP) J = J + 1 IE = IE_CELL(J) POS = POS_CELL(J) KSUM = KSUM + MAX(K(POS,IE),ZERO) IF ( ABS(K(POS,IE)) > KMAX ) KMAX = ABS(K(POS,IE)) END DO IF (KSUM > ZERO) THEN DTMAX_EXP = SI(IP)/KSUM ELSE DTMAX_EXP = 10.d14 END IF ! IF (KMAX > ZERO) THEN ! DTMAX_EXP = SI(IP)/KMAX ! Somewhat smaller due to the CRD approach ... ! ELSE ! DTMAX_EXP = 10E14 ! END IF IF (DTMAX_GLOBAL_EXP > DTMAX_EXP) DTMAX_GLOBAL_EXP = DTMAX_EXP END DO REST = ABS(MOD(DT4A/DTMAX_GLOBAL_EXP,ONE)) IF (REST > THR .AND. REST < ONEHALF) THEN ITER_EXP(IS,ID) = ABS(NINT(DT4A/DTMAX_GLOBAL_EXP)) + 1 ELSE ITER_EXP(IS,ID) = ABS(NINT(DT4A/DTMAX_GLOBAL_EXP)) END IF DTMAX = DTMAX_GLOBAL_EXP END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE EXPLICIT_N_SCHEME ( IS, ID ) USE DATAPOOL IMPLICIT NONE INTEGER, INTENT(IN) :: IS,ID ! ! local integer ! INTEGER :: IP, IE, IT, IP_TEST INTEGER :: I1, I2, I3, I, J, IMETHOD, IPOS INTEGER :: NI(3), POS ! ! local double ! REAL(rkind) :: UTILDE REAL(rkind) :: DTMAX_GLOBAL_EXP, DTMAX_EXP #ifdef MPI_PARALL_GRID REAL(rkind) :: DTMAX_GLOBAL_EXP_LOC #endif REAL(rkind) :: REST, TESTMIN REAL(rkind) :: LAMBDA(2), DT4AI REAL(rkind) :: FL11,FL12,FL21,FL22,FL31,FL32 REAL(rkind) :: KTMP(3) REAL(rkind) :: KKSUM(MNP), KKMAX(MNP), ST(MNP), N(MNE), U3(3) REAL(rkind) :: C(2,MNP), U(MNP), DTSI(MNP), CFLXY REAL(rkind) :: FLALL(3,MNE), UTILDEE(MNE) REAL(rkind) :: FL111, FL112, FL211, FL212, FL311, FL312 REAL(rkind) :: KELEM(3,MNE) REAL(rkind) :: CBAR_1_1(2), CBAR_1_2(2) REAL(rkind) :: CBAR_2_1(2), CBAR_2_2(2) REAL(rkind) :: CBAR_3_1(2), CBAR_3_2(2) REAL(rkind) :: Ftest(MNP) ! ! local parameter ! REAL(rkind) :: TMP ! ! Calculate phase speeds for the certain spectral component ... ! CALL CADVXY(IS,ID,C) IP_TEST = 180 ! ! Calculate K-Values and contour based quantities ... ! !!$OMP PARALLEL DO DEFAULT(NONE) & !!$OMP& PRIVATE(IE,I1,I2,I3,LAMBDA,KTMP,TMP,FL11,FL12,FL21,FL22,FL31,FL32,FL111,FL112,FL211,FL212,FL311,FL312) & !!$OMP& SHARED(MNE,INE,IEN,C,KELEM,FLALL,N) DO IE = 1, MNE ! IF (IE_IS_STEADY(IE) .GT. 2) THEN ! WRITE(DBG%FHNDL,*) '1st IE LOOP CYCLE', IE, IE_IS_STEADY(IE) ! CYCLE ! ENDIF I1 = INE(1,IE) I2 = INE(2,IE) I3 = INE(3,IE) LAMBDA(1) = ONESIXTH *(C(1,I1)+C(1,I2)+C(1,I3)) LAMBDA(2) = ONESIXTH *(C(2,I1)+C(2,I2)+C(2,I3)) KELEM(1,IE) = LAMBDA(1) * IEN(1,IE) + LAMBDA(2) * IEN(2,IE) KELEM(2,IE) = LAMBDA(1) * IEN(3,IE) + LAMBDA(2) * IEN(4,IE) KELEM(3,IE) = LAMBDA(1) * IEN(5,IE) + LAMBDA(2) * IEN(6,IE) #ifdef positivity CBAR_1_1 = 0.5 * (C(:,I1) + C(:,I3)) CBAR_1_2 = 0.5 * (C(:,I1) + C(:,I2)) CBAR_2_1 = 0.5 * (C(:,I2) + C(:,I3)) CBAR_2_2 = 0.5 * (C(:,I2) + C(:,I1)) CBAR_3_1 = 0.5 * (C(:,I2) + C(:,I3)) CBAR_3_2 = 0.5 * (C(:,I1) + C(:,I3)) KELEM(1,IE) = -DOT_PRODUCT(CBAR_1_1,IEN(3:4,IE)) -DOT_PRODUCT(CBAR_1_2,IEN(5:6,IE)) KELEM(2,IE) = -DOT_PRODUCT(CBAR_2_1,IEN(1:2,IE)) -DOT_PRODUCT(CBAR_2_2,IEN(5:6,IE)) KELEM(3,IE) = -DOT_PRODUCT(CBAR_3_1,IEN(1:2,IE)) -DOT_PRODUCT(CBAR_3_2,IEN(3:4,IE)) #endif KTMP = KELEM(:,IE) TMP = SUM(MIN(ZERO,KTMP)) N(IE) = - ONE/MIN(-THR,TMP) KELEM(:,IE) = MAX(ZERO,KTMP) FL11 = C(1,I2) * IEN(1,IE) + C(2,I2) * IEN(2,IE) FL12 = C(1,I3) * IEN(1,IE) + C(2,I3) * IEN(2,IE) FL21 = C(1,I3) * IEN(3,IE) + C(2,I3) * IEN(4,IE) FL22 = C(1,I1) * IEN(3,IE) + C(2,I1) * IEN(4,IE) FL31 = C(1,I1) * IEN(5,IE) + C(2,I1) * IEN(6,IE) FL32 = C(1,I2) * IEN(5,IE) + C(2,I2) * IEN(6,IE) FL111 = TWO*FL11+FL12 FL112 = TWO*FL12+FL11 FL211 = TWO*FL21+FL22 FL212 = TWO*FL22+FL21 FL311 = TWO*FL31+FL32 FL312 = TWO*FL32+FL31 FLALL(1,IE) = (FL311 + FL212) * ONESIXTH + KELEM(1,IE) FLALL(2,IE) = (FL111 + FL312) * ONESIXTH + KELEM(2,IE) FLALL(3,IE) = (FL211 + FL112) * ONESIXTH + KELEM(3,IE) END DO ! If the current field or water level changes estimate the iteration ! number based on the new flow field and the CFL number of the scheme IF (LCALC) THEN ! KKSUM = ZERO ! DO IE = 1, MNE ! IF (IE_IS_STEADY(IE) .GT. 2) THEN ! CYCLE ! ENDIF ! NI = INE(:,IE) ! KKSUM(NI) = KKSUM(NI) + KELEM(:,IE) ! END DO !AR: Experimental ... improves speed by 20% but maybe unstable in !certain situations ... must be checked thoroughly KKMAX = ZERO KKSUM = ZERO J = 0 DO IP = 1, MNP DO I = 1, CCON(IP) J = J + 1 IE = IE_CELL(J) POS = POS_CELL(J) KKSUM(IP) = KKSUM(IP) + MAX(KELEM(POS,IE),ZERO) ! IF ( ABS(KELEM(POS,IE)) > KKMAX(IP) ) KKMAX(IP) = ABS(KELEM(POS,IE)) END DO END DO #ifdef MPI_PARALL_GRID DTMAX_GLOBAL_EXP = VERYLARGE DTMAX_GLOBAL_EXP_LOC = VERYLARGE DO IP = 1, NP_RES ! IF (IP_IS_STEADY(IP) .GT. 2) THEN ! CYCLE ! ENDIF DTMAX_EXP = SI(IP)/MAX(THR,KKSUM(IP)) !WRITE(DBG%FHNDL,'(I10,3F15.6)') IP, SI(IP), KKSUM(IP), DEPTH(IP), DTMAX_EXP IF (LCFL) THEN CFLCXY(1,IP) = MAX(CFLCXY(1,IP), C(1,IP)) CFLCXY(2,IP) = MAX(CFLCXY(2,IP), C(2,IP)) CFLCXY(3,IP) = DTMAX_EXP/DT4A END IF DTMAX_GLOBAL_EXP_LOC = MIN(DTMAX_GLOBAL_EXP_LOC,DTMAX_EXP) END DO CALL MPI_ALLREDUCE(DTMAX_GLOBAL_EXP_LOC,DTMAX_GLOBAL_EXP,1,rtype,MPI_MIN,COMM,IERR) !WRITE(STAT%FHNDL,'(2I10,2F15.4)') IS, ID, DTMAX_GLOBAL_EXP, DT4A/DTMAX_GLOBAL_EXP #else DTMAX_GLOBAL_EXP = VERYLARGE DO IP = 1, MNP IF (IOBP(IP) .NE. 0) CYCLE ! IF (IP_IS_STEADY(IP) .GT. 2) THEN ! CYCLE ! ENDIF DTMAX_EXP = SI(IP)/MAX(THR,KKSUM(IP)) !DTMAX_EXP = SI(IP)/MAX(THR,KKMAX(IP)) IF (LCFL) THEN CFLCXY(1,IP) = MAX(CFLCXY(1,IP), C(1,IP)) CFLCXY(2,IP) = MAX(CFLCXY(2,IP), C(2,IP)) CFLCXY(3,IP) = KKSUM(IP) END IF DTMAX_GLOBAL_EXP = MIN ( DTMAX_GLOBAL_EXP, DTMAX_EXP) !WRITE(22227,*) IP, CCON(IP), SI(IP) !IF (IP == 24227 .AND. IS == 1) WRITE(DBG%FHNDL,'(2I10,6F20.8)') IP, ID, XP(IP), YP(IP), SI(IP), KKSUM(IP), DEP(IP), CFLCXY(3,IP) END DO !WRITE(STAT%FHNDL,'(2I10,2F15.4)') IS, ID, DTMAX_GLOBAL_EXP, DT4A/DTMAX_GLOBAL_EXP !AR: Makes very strange error in the code ... #endif CFLXY = DT4A/DTMAX_GLOBAL_EXP REST = ABS(MOD(CFLXY,ONE)) IF (REST .LT. THR) THEN ITER_EXP(IS,ID) = ABS(NINT(CFLXY)) ELSE IF (REST .GT. THR .AND. REST .LT. ONEHALF) THEN ITER_EXP(IS,ID) = ABS(NINT(CFLXY)) + 1 ELSE ITER_EXP(IS,ID) = ABS(NINT(CFLXY)) END IF END IF DT4AI = DT4A/ITER_EXP(IS,ID) DTSI(:) = DT4AI/SI(:) U = AC2(IS,ID,:) #ifdef DEBUG_COHERENCY_FLUCT WRITE(STAT%FHNDL,*) 'IS=', IS, ' ID=', ID CALL Print_SumScalar(SI, "SI at start of EXPLICIT_N_SCHEME") CALL Print_SumScalar(DTSI, "DTSI at start of EXPLICIT_N_SCHEME") CALL Print_SumScalar(U, "U at start of EXPLICIT_N_SCHEME") Ftest=MyREAL(IOBWB) CALL Print_SumScalar(Ftest, "IOBWB at start of EXPLICIT_N_SCHEME") Ftest=MyREAL(IOBPD(ID,:)) CALL Print_SumScalar(Ftest, "IOBPD at start of EXPLICIT_N_SCHEME") Ftest=MyREAL(IOBDP) CALL Print_SumScalar(Ftest, "IOBDP at start of EXPLICIT_N_SCHEME") Ftest=C(1,:) CALL Print_SumScalar(Ftest, "C(1,:) at start of EXPLICIT_N_SCHEME") Ftest=C(2,:) CALL Print_SumScalar(Ftest, "C(2,:) at start of EXPLICIT_N_SCHEME") #endif IF (LADVTEST) THEN CALL CHECKCONS(U,SUMAC1) END IF ! ! Loop over all sub time steps, all quantities in this loop depend on the solution U itself !!! ! ! WRITE(STAT%FHNDL,'(3I10,4F15.4)') IS, ID, ITER_EXP(IS,ID), SQRT(MAXVAL(C(1,:))**2+MAXVAL(C(2,:))**2), & ! & SQRT(MAXVAL(C(1,:))**2+MAXVAL(C(2,:))**2)*DT4A/MINVAL(EDGELENGTH), MAXVAL(CG(IS,:)), SQRT(G9*MAXVAL(DEP)) IMETHOD = 1 IF (IMETHOD == 1) THEN DO IT = 1, ITER_EXP(IS,ID) #ifdef DEBUG_COHERENCY_FLUCT WRITE(STAT%FHNDL,*) 'IT=', IT #endif ST = ZERO ! Init. ... only used over the residual nodes see IP loop DO IE = 1, MNE ! IF (IE_IS_STEADY(IE) .GT. 2) THEN ! WRITE(DBG%FHNDL,*) '2nd IE LOOP CYCLE', IT, IE, IE_IS_STEADY(IE) ! CYCLE ! ENDIF NI = INE(:,IE) U3 = U(NI) UTILDE = N(IE) * (DOT_PRODUCT(FLALL(:,IE),U3)) !* IOBED(ID,IE) ST(NI) = ST(NI) + KELEM(:,IE) * (U3 - UTILDE) ! the 2nd term are the theta values of each node ... END DO #ifdef DEBUG_COHERENCY_FLUCT # ifdef MPI_PARALL_GRID CALL EXCHANGE_P2D(ST) ! Simply for debugging purposes # endif CALL Print_SumScalar(ST, "ST used in update of U") #endif DO IP = 1, MNP ! IF (IP_IS_STEADY(IP) .GT. 2) THEN ! WRITE(DBG%FHNDL,*) '1st IP LOOP CYCLE', IT, IP, IP_IS_STEADY(IP) ! CYCLE ! ENDIF ! IF (IOBP(IP) .NE. 0 .AND. CCON(IP) .GT. 3) THEN ! TESTMIN = U(IP)-DTSI(IP)*ST(IP) ! IF (TESTMIN .LT. ZERO) WRITE(99999,*) TESTMIN U(IP) = MAX(ZERO,U(IP)-DTSI(IP)*ST(IP)*IOBWB(IP))*IOBPD(ID,IP)*IOBDP(IP) ! ELSE IF (IOBP(IP) .EQ. 0) THEN ! U(IP) = MAX(ZERO,U(IP)-DTSI(IP)*ST(IP)*IOBWB(IP))*IOBPD(ID,IP)*IOBDP(IP) ! ENDIF ENDDO ! WRITE(*,'(2I10,F20.10,2I20,F20.10)') ID, IS, U(IP_TEST), IOBPD(ID,IP_TEST), IOBDP(IP_TEST), DEP(IP_TEST) #ifdef DEBUG_COHERENCY_FLUCT CALL Print_SumScalar(U, "U after the iteration") #endif #ifdef MPI_PARALL_GRID CALL EXCHANGE_P2D(U) ! Exchange after each update of the res. domain #endif #ifdef DEBUG_COHERENCY_FLUCT CALL Print_SumScalar(U, "U after the exchange") #endif END DO ! ----> End Iteration ELSE IF (IMETHOD == 2) THEN DO IT = 1, ITER_EXP(IS,ID) UTILDEE = N*(FLALL(1,:)*U(INE(1,:))+FLALL(2,:)*U(INE(2,:))+FLALL(3,:)*U(INE(3,:))) ST = ZERO J = 0 DO IP = 1, MNP DO I = 1, CCON(IP) IE = IE_CELL2(IP,I) IPOS = POS_CELL2(IP,I) ST(IP) = ST(IP) + KELEM(IPOS,IE) * (U(IP) - UTILDEE(IE)) END DO END DO U = MAX(ZERO,U-DTSI*ST*IOBWB)*IOBPD(ID,:)*IOBDP(:) #ifdef MPI_PARALL_GRID CALL EXCHANGE_P2D(U) ! Exchange after each update of the res. domain #endif END DO ELSE IF (IMETHOD == 3) THEN DO IT = 1, ITER_EXP(IS,ID) ST = 0.d0 ! Init. ... only used over the residual nodes see IP loop DO IE = 1, MNE NI = INE(:,IE) U3 = U(NI) UTILDE = N(IE)*(DOT_PRODUCT(FLALL(:,IE),U3*IOBPD(ID,NI))) !UTILDE = N(IE)*(DOT_PRODUCT(FLALL(:,IE),U3)) ST(NI) = ST(NI)*IOBPD(ID,NI)+KELEM(:,IE)*(U3 - UTILDE) !ST(NI) = ST(NI)+KELEM(:,IE)*(U3 - UTILDE) END DO U = MAX(0.d0,U-DTSI*ST*MyREAL(IOBWB))!*DBLE(IOBPD(ID,:)) #ifdef MPI_PARALL_GRID CALL EXCHANGE_P2D(U) #endif END DO ! ----> End Iteration END IF ! IMETHOD AC2(IS,ID,:) = U IF (LADVTEST) THEN WRITE(4001) SNGL(RTIME) WRITE(4001) (SNGL(C(1,IP)), SNGL(C(2,IP)), SNGL(AC2(1,1,IP)), & & IP = 1, MNP) CALL CHECKCONS(U,SUMAC2) IF (MINVAL(U) .LT. MINTEST) MINTEST = MINVAL(U) WRITE (DBG%FHNDL,*) 'VOLUMES AT T0, T1 and T2',SUMACt0, & & SUMAC1, SUMAC2, MINTEST WRITE (DBG%FHNDL,*) 'VOLUME ERROR: TOTAL and ACTUAL', & & 100.0_rkind-((SUMACt0-SUMAC2)/SUMACt0)*100.0_rkind, & & 100.0_rkind- ((SUMAC1-SUMAC2)/SUMAC1)*100.0_rkind END IF END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE EXPLICIT_PSI_SCHEME ( IS, ID ) USE DATAPOOL IMPLICIT NONE INTEGER, INTENT(IN) :: IS,ID ! ! local integer ! INTEGER :: IP, IE, IT INTEGER :: I1, I2, I3, K INTEGER :: NI(3) ! ! local double ! REAL(rkind) :: FT REAL(rkind) :: UTILDE REAL(rkind) :: DTMAX_GLOBAL_EXP, DTMAX_EXP #ifdef MPI_PARALL_GRID REAL(rkind) :: DTMAX_GLOBAL_EXP_LOC #endif REAL(rkind) :: REST REAL(rkind) :: LAMBDA(2), DT4AI REAL(rkind) :: FL11,FL12,FL21,FL22,FL31,FL32 REAL(rkind) :: THETA_L(3) REAL(rkind) :: KTMP(3) REAL(rkind) :: BET1(3), BETAHAT(3) REAL(rkind) :: KKSUM(MNP), ST(MNP), N(MNE) REAL(rkind) :: C(2,MNP), U(MNP), DTSI(MNP), CFLXY REAL(rkind) :: FLALL(3,MNE) REAL(rkind) :: FL111, FL112, FL211, FL212, FL311, FL312 REAL(rkind) :: KELEM(3,MNE) ! ! local parameter ! REAL(rkind) :: TMP U(:) = AC2(IS,ID,:) ! ! Calculate phase speeds for the certain spectral component ... ! CALL CADVXY(IS,ID,C) ! ! Calculate K-Values and contour based quantities ... ! DO IE = 1, MNE I1 = INE(1,IE) I2 = INE(2,IE) I3 = INE(3,IE) LAMBDA(1) = ONESIXTH *(C(1,I1)+C(1,I2)+C(1,I3)) LAMBDA(2) = ONESIXTH *(C(2,I1)+C(2,I2)+C(2,I3)) KELEM(1,IE) = LAMBDA(1) * IEN(1,IE) + LAMBDA(2) * IEN(2,IE) KELEM(2,IE) = LAMBDA(1) * IEN(3,IE) + LAMBDA(2) * IEN(4,IE) KELEM(3,IE) = LAMBDA(1) * IEN(5,IE) + LAMBDA(2) * IEN(6,IE) KTMP = KELEM(:,IE) TMP = SUM(MIN(ZERO,KTMP)) N(IE) = - ONE/MIN(-THR,TMP) KELEM(:,IE) = MAX(ZERO,KTMP) FL11 = C(1,I2) * IEN(1,IE) + C(2,I2) * IEN(2,IE) FL12 = C(1,I3) * IEN(1,IE) + C(2,I3) * IEN(2,IE) FL21 = C(1,I3) * IEN(3,IE) + C(2,I3) * IEN(4,IE) FL22 = C(1,I1) * IEN(3,IE) + C(2,I1) * IEN(4,IE) FL31 = C(1,I1) * IEN(5,IE) + C(2,I1) * IEN(6,IE) FL32 = C(1,I2) * IEN(5,IE) + C(2,I2) * IEN(6,IE) FL111 = TWO*FL11+FL12 FL112 = TWO*FL12+FL11 FL211 = TWO*FL21+FL22 FL212 = TWO*FL22+FL21 FL311 = TWO*FL31+FL32 FL312 = TWO*FL32+FL31 FLALL(1,IE) = FL311 + FL212 FLALL(2,IE) = FL111 + FL312 FLALL(3,IE) = FL211 + FL112 END DO ! If the current field or water level changes estimate the iteration ! number based on the new flow field and the CFL number of the scheme IF (LCALC) THEN KKSUM = ZERO DO IE = 1, MNE NI = INE(:,IE) KKSUM(NI) = KKSUM(NI) + KELEM(:,IE) END DO ! #ifdef MPI_PARALL_GRID DTMAX_GLOBAL_EXP = VERYLARGE DTMAX_GLOBAL_EXP_LOC = VERYLARGE DO IP = 1, NP_RES DTMAX_EXP = SI(IP)/MAX(THR,KKSUM(IP)) IF (LCFL) THEN CFLCXY(1,IP) = MAX(CFLCXY(1,IP), C(1,IP)) CFLCXY(2,IP) = MAX(CFLCXY(2,IP), C(2,IP)) CFLCXY(3,IP) = MAX(CFLCXY(3,IP), DTMAX_EXP) END IF DTMAX_GLOBAL_EXP_LOC=MIN(DTMAX_GLOBAL_EXP_LOC, DTMAX_EXP) END DO CALL MPI_ALLREDUCE(DTMAX_GLOBAL_EXP_LOC,DTMAX_GLOBAL_EXP, & & 1,rtype,MPI_MIN,comm,ierr) #else DTMAX_GLOBAL_EXP = VERYLARGE DO IP = 1, MNP DTMAX_EXP = SI(IP)/MAX(THR,KKSUM(IP)) IF (LCFL) THEN CFLCXY(1,IP) = MAX(CFLCXY(1,IP), C(1,IP)) CFLCXY(2,IP) = MAX(CFLCXY(2,IP), C(2,IP)) CFLCXY(3,IP) = MAX(CFLCXY(3,IP), DTMAX_EXP) END IF DTMAX_GLOBAL_EXP = MIN ( DTMAX_GLOBAL_EXP, DTMAX_EXP) END DO #endif ! ! ITER_EXP(IS,ID) is the number of sub time step in order to fullfill ! the CFL number .LT. 1 for the certain wave component ... ! CFLXY = DT4A/DTMAX_GLOBAL_EXP REST = ABS(MOD(CFLXY,ONE)) IF (REST .LT. THR) THEN ITER_EXP(IS,ID) = ABS(NINT(CFLXY)) ELSE IF (REST .GT. THR .AND. REST .LT. ONEHALF) THEN ITER_EXP(IS,ID) = ABS(NINT(CFLXY)) + 1 ELSE ITER_EXP(IS,ID) = ABS(NINT(CFLXY)) END IF ITER_EXP(IS,ID) = MAX(1,ITER_EXP(IS,ID)) END IF DT4AI = DT4A/ITER_EXP(IS,ID) DTSI(:) = DT4AI/SI(:) ! ! Loop over all sub time steps, all quantities in this loop depend ! on the solution U itself !!! ! U(:) = AC2(IS,ID,:) DO IT = 1, ITER_EXP(IS,ID) ST = ZERO DO IE = 1, MNE NI = INE(:,IE) FT = -ONESIXTH*DOT_PRODUCT(U(NI),FLALL(:,IE)) UTILDE = N(IE) * ( DOT_PRODUCT(KELEM(:,IE),U(NI)) - FT ) THETA_L(:) = KELEM(:,IE) * (U(NI) - UTILDE) IF (ABS(FT) .GT. ZERO) THEN BET1(:) = THETA_L(:)/FT IF (ANY( BET1 .LT. ZERO) ) THEN BETAHAT(1) = BET1(1) + ONEHALF * BET1(2) BETAHAT(2) = BET1(2) + ONEHALF * BET1(3) BETAHAT(3) = BET1(3) + ONEHALF * BET1(1) BET1(1) = MAX(ZERO,MIN(BETAHAT(1), & & ONE-BETAHAT(2),ONE)) BET1(2) = MAX(ZERO,MIN(BETAHAT(2), & & ONE-BETAHAT(3),ONE)) BET1(3) = MAX(ZERO,MIN(BETAHAT(3), & & ONE-BETAHAT(1),ONE)) THETA_L(:) = FT * BET1 END IF ELSE THETA_L(:) = ZERO END IF ! the 2nd term are the theta values of each node ... ST(NI) = ST(NI) + THETA_L END DO U = MAX(ZERO,U-DTSI*ST*IOBWB)*IOBPD(ID,:) DO IP = 1, MNP U(IP) = MAX(ZERO,U(IP)-DTSI(IP)*ST(IP)*IOBWB(IP))*IOBPD(ID,IP)*IOBDP(IP) END DO #ifdef MPI_PARALL_GRID CALL EXCHANGE_P2D(U) ! Exchange after each update of the res. domain #endif END DO ! ----> End Iteration AC2(IS,ID,:) = U(:) IF (LADVTEST) THEN WRITE(4001) SNGL(RTIME) WRITE(4001) (SNGL(C(1,IP)), SNGL(C(2,IP)), SNGL(AC2(1,1,IP)), & & IP = 1, MNP) CALL CHECKCONS(U,SUMAC2) IF (MINVAL(U) .LT. MINTEST) MINTEST = MINVAL(U) WRITE (DBG%FHNDL,*) 'VOLUMES AT T0, T1 and T2',SUMACt0, & & SUMAC1, SUMAC2, MINTEST WRITE (DBG%FHNDL,*) 'VOLUME ERROR: TOTAL and ACTUAL', & & 100.0_rkind-((SUMACt0-SUMAC2)/SUMACt0)*100.0_rkind, & & 100.0_rkind- ((SUMAC1-SUMAC2)/SUMAC1)*100.0_rkind END IF END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE EXPLICIT_LFPSI_SCHEME(IS,ID) USE DATAPOOL IMPLICIT NONE INTEGER, INTENT(IN) :: IS,ID ! ! local integer ! INTEGER :: IP, IE, IT, I, J INTEGER :: I1, I2, I3 INTEGER :: NI(3) ! ! local double ! REAL(rkind) :: FT REAL(rkind) :: UTILDE REAL(rkind) :: DTMAX_GLOBAL_EXP, DTMAX_EXP REAL(rkind) :: REST REAL(rkind) :: TMP(3), TMP1 REAL(rkind) :: LAMBDA(2), DT4AI REAL(rkind) :: BET1(3), BETAHAT(3), BL REAL(rkind) :: FL11,FL12,FL21,FL22,FL31,FL32 REAL(rkind) :: THETA_L(3,MNE), THETA_H(3), THETA_ACE(3,MNE) REAL(rkind) :: UTMP(3) REAL(rkind) :: WII(2,MNP), UL(MNP), USTARI(2,MNP) REAL(rkind) :: KKSUM(MNP), ST(MNP), PM(MNP), PP(MNP), UIM(MNE) REAL(rkind) :: UIP(MNE), UIPIP(MNP), UIMIP(MNP) REAL(rkind) :: C(2,MNP), U(MNP), DTSI(MNP), CFLXY, N(MNE) REAL(rkind) :: FL111, FL112, FL211, FL212, FL311, FL312 REAL(rkind) :: KELEM(3,MNE), FLALL(3,MNE) #ifdef MPI_PARALL_GRID REAL(rkind) :: DTMAX_GLOBAL_EXP_LOC #endif ! ! local parameter ! BL = ZERO ! ! Calculate phase speeds for the certain spectral component ... ! CALL CADVXY(IS,ID,C) ! ! Calculate K-Values and contour based quantities ... ! DO IE = 1, MNE I1 = INE(1,IE) I2 = INE(2,IE) I3 = INE(3,IE) LAMBDA(1) = ONESIXTH *(C(1,I1)+C(1,I2)+C(1,I3)) LAMBDA(2) = ONESIXTH *(C(2,I1)+C(2,I2)+C(2,I3)) KELEM(1,IE) = LAMBDA(1) * IEN(1,IE) + LAMBDA(2) * IEN(2,IE) KELEM(2,IE) = LAMBDA(1) * IEN(3,IE) + LAMBDA(2) * IEN(4,IE) KELEM(3,IE) = LAMBDA(1) * IEN(5,IE) + LAMBDA(2) * IEN(6,IE) N(IE) = - ONE/MIN(-THR,SUM(MIN(ZERO,KELEM(:,IE)))) FL11 = C(1,I2) * IEN(1,IE) + C(2,I2) * IEN(2,IE) FL12 = C(1,I3) * IEN(1,IE) + C(2,I3) * IEN(2,IE) FL21 = C(1,I3) * IEN(3,IE) + C(2,I3) * IEN(4,IE) FL22 = C(1,I1) * IEN(3,IE) + C(2,I1) * IEN(4,IE) FL31 = C(1,I1) * IEN(5,IE) + C(2,I1) * IEN(6,IE) FL32 = C(1,I2) * IEN(5,IE) + C(2,I2) * IEN(6,IE) FL111 = 2*FL11+FL12 FL112 = 2*FL12+FL11 FL211 = 2*FL21+FL22 FL212 = 2*FL22+FL21 FL311 = 2*FL31+FL32 FL312 = 2*FL32+FL31 FLALL(1,IE) = FL311 + FL212 FLALL(2,IE) = FL111 + FL312 FLALL(3,IE) = FL211 + FL112 END DO ! If the current field or water level changes estimate the iteration ! number based on the new flow field and the CFL number of the scheme IF (LCALC) THEN KKSUM = ZERO DO IE = 1, MNE NI = INE(:,IE) KKSUM(NI) = KKSUM(NI) + MAX(ZERO,KELEM(:,IE)) END DO #ifdef MPI_PARALL_GRID DTMAX_GLOBAL_EXP = VERYLARGE DTMAX_GLOBAL_EXP_LOC = VERYLARGE DO IP = 1, NP_RES DTMAX_EXP = SI(IP)/MAX(THR,KKSUM(IP)) ! DTMAX_EXP = MAX( ABS(IOBP(IP)*VERYLARGE), SI(IP)/MAX(THR,KKSUM(IP)*IOBPD(ID,IP)) ) IF (LCFL) THEN CFLCXY(1,IP) = MAX(CFLCXY(1,IP), C(1,IP)) CFLCXY(2,IP) = MAX(CFLCXY(2,IP), C(2,IP)) CFLCXY(3,IP) = MAX(CFLCXY(3,IP), DTMAX_EXP) END IF DTMAX_GLOBAL_EXP_LOC=MIN ( DTMAX_GLOBAL_EXP_LOC, DTMAX_EXP) END DO CALL MPI_ALLREDUCE(DTMAX_GLOBAL_EXP_LOC,DTMAX_GLOBAL_EXP, 1,rtype,MPI_MIN,comm,ierr) #else DTMAX_GLOBAL_EXP = VERYLARGE DO IP = 1, MNP DTMAX_EXP = SI(IP)/MAX(THR,KKSUM(IP)) IF (LCFL) THEN CFLCXY(1,IP) = MAX(CFLCXY(1,IP), C(1,IP)) CFLCXY(2,IP) = MAX(CFLCXY(2,IP), C(2,IP)) CFLCXY(3,IP) = MAX(CFLCXY(3,IP), DTMAX_EXP) END IF DTMAX_GLOBAL_EXP = MIN ( DTMAX_GLOBAL_EXP, DTMAX_EXP) END DO #endif CFLXY = DT4A/DTMAX_GLOBAL_EXP REST = ABS(MOD(CFLXY,1._rkind)) IF (REST .LT. THR) THEN ITER_EXP(IS,ID) = ABS(NINT(CFLXY)) ELSE IF (REST .GT. THR .AND. REST .LT. ONEHALF) THEN ITER_EXP(IS,ID) = ABS(NINT(CFLXY)) + 1 ELSE ITER_EXP(IS,ID) = ABS(NINT(CFLXY)) END IF END IF ! LCALC DT4AI = DT4A/ITER_EXP(IS,ID) DTSI(:) = DT4AI/SI(:) U(:) = AC2(IS,ID,:) UL = U #ifdef MPI_PARALL_GRID CALL EXCHANGE_P2D(U) CALL EXCHANGE_P2D(UL) #endif ! ! Loop over all sub time steps, all quantities in this loop depend ! on the solution U itself !!! ! DO IT = 1, ITER_EXP(IS,ID) ! ! Element loop ! ST = ZERO PM = ZERO PP = ZERO DO IE = 1, MNE NI = INE(:,IE) UTMP = U(NI) FT = -ONESIXTH*DOT_PRODUCT(UTMP,FLALL(:,IE)) TMP = MAX(ZERO,KELEM(:,IE)) UTILDE = N(IE) * ( DOT_PRODUCT(TMP,UTMP) - FT ) THETA_L(:,IE) = TMP * ( UTMP - UTILDE ) IF (ABS(FT) .GT. ZERO) THEN BET1(:) = THETA_L(:,IE)/FT IF (ANY( BET1 .LT. ZERO) ) THEN BETAHAT(1) = BET1(1) + ONEHALF * BET1(2) BETAHAT(2) = BET1(2) + ONEHALF * BET1(3) BETAHAT(3) = BET1(3) + ONEHALF * BET1(1) BET1(1) = MAX(ZERO,MIN(BETAHAT(1),ONE-BETAHAT(2),ONE)) BET1(2) = MAX(ZERO,MIN(BETAHAT(2),ONE-BETAHAT(3),ONE)) BET1(3) = MAX(ZERO,MIN(BETAHAT(3),ONE-BETAHAT(1),ONE)) THETA_L(:,IE) = FT * BET1 END IF ELSE THETA_L(:,IE) = ZERO END IF ST(NI) = ST(NI) + THETA_L(:,IE) THETA_H = (ONETHIRD+DT4AI/(TWO*TRIA(IE)) * KELEM(:,IE) ) * FT ! LAX ! THETA_H = (ONETHIRD+TWOTHIRD*KELEM(:,IE)/SUM(MAX(ZERO,KELEM(:,IE))))*FT ! CENTRAL THETA_ACE(:,IE) = THETA_H-THETA_L(:,IE) PP(NI) = PP(NI) + MAX(ZERO, -THETA_ACE(:,IE)) * DTSI(NI) PM(NI) = PM(NI) + MIN(ZERO, -THETA_ACE(:,IE)) * DTSI(NI) END DO DO IP = 1, MNP UL(IP) = MAX(ZERO,U(IP)-DTSI(IP)*ST(IP)*IOBWB(IP))*IOBPD(ID,IP)*IOBDP(IP) ENDDO #ifdef MPI_PARALL_GRID CALL EXCHANGE_P2D(UL) ! Exchange after each update of the res. domain CALL EXCHANGE_P2D(U) ! Exchange after each update of the res. domain #endif USTARI(1,:) = MAX(UL,U) * IOBPD(ID,:) USTARI(2,:) = MIN(UL,U) * IOBPD(ID,:) UIP = 0. UIM = 0. DO IE = 1, MNE NI = INE(:,IE) UIP(IE) = MAXVAL(USTARI(1,NI)) UIM(IE) = MINVAL(USTARI(2,NI)) END DO J = 0 ! Counter ... UIPIP = 0. UIMIP = 0. DO IP = 1, MNP DO I = 1, CCON(IP) J = J + 1 IE = IE_CELL(J) UIPIP(IP) = MAX(UIPIP(IP),UIP(IE)) UIMIP(IP) = MIN(UIMIP(IP),UIM(IE)) ENDDO END DO !I: loop over connected elemen DO IP = 1, MNP WII(1,IP) = MIN(ONE,(UIPIP(IP)-UL(IP))/MAX( THR,PP(IP))) WII(2,IP) = MIN(ONE,(UIMIP(IP)-UL(IP))/MIN(-THR,PM(IP))) IF (ABS(PP(IP)) .LT. THR) WII(1,IP) = ZERO IF (ABS(PM(IP)) .LT. THR) WII(2,IP) = ZERO END DO ST = ZERO DO IE = 1, MNE I1 = INE(1,IE) I2 = INE(2,IE) I3 = INE(3,IE) IF (THETA_ACE(1,IE) .LT. ZERO) THEN TMP(1) = WII(1,I1) ELSE TMP(1) = WII(2,I1) END IF IF (THETA_ACE(2,IE) .LT. ZERO) THEN TMP(2) = WII(1,I2) ELSE TMP(2) = WII(2,I2) END IF IF (THETA_ACE(3,IE) .LT. ZERO) THEN TMP(3) = WII(1,I3) ELSE TMP(3) = WII(2,I3) END IF TMP1 = MINVAL(TMP) ST(I1) = ST(I1) + THETA_ACE(1,IE) * TMP1! * (ONE - BL) + BL * THETA_L(1,IE) ST(I2) = ST(I2) + THETA_ACE(2,IE) * TMP1! * (ONE - BL) + BL * THETA_L(2,IE) ST(I3) = ST(I3) + THETA_ACE(3,IE) * TMP1! * (ONE - BL) + BL * THETA_L(3,IE) END DO !U = MAX(ZERO,UL-DTSI*ST*IOBWB)*IOBPD(ID,:) DO IP = 1, MNP U(IP) = MAX(ZERO,UL(IP)-DTSI(IP)*ST(IP)*IOBWB(IP))*IOBPD(ID,IP)*IOBDP(IP) ENDDO #ifdef MPI_PARALL_GRID CALL EXCHANGE_P2D(U) ! Exchange after each update of the res. domain #endif END DO ! ----> End Iteration AC2(IS,ID,:) = UL(:) IF (LADVTEST) THEN WRITE(4001) SNGL(RTIME) WRITE(4001) (SNGL(C(1,IP)), SNGL(C(2,IP)), SNGL(AC2(1,1,IP)), & & IP = 1, MNP) CALL CHECKCONS(U,SUMAC2) IF (MINVAL(U) .LT. MINTEST) MINTEST = MINVAL(U) WRITE (*,*) 'VOLUMES AT T0, T1 and T2',SUMACt0, & & SUMAC1, SUMAC2, MINTEST WRITE (*,*) 'VOLUME ERROR: TOTAL and ACTUAL', & & 100.0_rkind-((SUMACt0-SUMAC2)/SUMACt0)*100.0_rkind, & & 100.0_rkind- ((SUMAC1-SUMAC2)/SUMAC1)*100.0_rkind END IF END SUBROUTINE !********************************************************************** !* !********************************************************************** SUBROUTINE EIMPS_V1( IS, ID) USE DATAPOOL IMPLICIT NONE INTEGER, INTENT(IN) :: IS,ID INTEGER :: I, J INTEGER :: IP, IPGL1, IE, POS INTEGER :: I1, I2, I3 REAL(rkind) :: DTK, TMP3 REAL(rkind) :: LAMBDA(2), GTEMP1, GTEMP2, FLHAB REAL(rkind) :: FL11, FL12, FL21, FL22, FL31, FL32 REAL(rkind):: CRFS(3), K1, KM(3), K(3), TRIA03, DELFL, USFM REAL(rkind) :: DELTAL(3,MNE) REAL(rkind) :: KP(3,MNE), NM(MNE) REAL(rkind) :: U(MNP), C(2,MNP) REAL(rkind) :: X(MNP) REAL(rkind) :: B(MNP) INTEGER :: IPAR(16) INTEGER :: IERROR INTEGER :: IWKSP(20*MNP) INTEGER :: FLJU(MNP) INTEGER :: FLJAU(NNZ+1) REAL(rkind) :: FPAR(16) REAL(rkind) :: WKSP( 20*MNP ) REAL(rkind) :: AU(NNZ+1) REAL(rkind) :: INIU(MNP) REAL(rkind) :: ASPAR(NNZ), LIMFAC REAL(rkind) :: TIME1, TIME2, TIME3, TIME4 REAL(rkind) :: TEMP, DELT, XIMP, DELT5 INTEGER :: POS_TRICK(3,2) external bcgstab external gmres #ifdef TIMINGS ! CALL WAV_MY_WTIME(TIME1) #endif IWKSP = 0 WKSP = ZERO POS_TRICK(1,1) = 2 POS_TRICK(1,2) = 3 POS_TRICK(2,1) = 3 POS_TRICK(2,2) = 1 POS_TRICK(3,1) = 1 POS_TRICK(3,2) = 2 CALL CADVXY(IS,ID,C) ! ! Calculate countour integral quantities ... ! DO IE = 1, MNE I1 = INE(1,IE) I2 = INE(2,IE) I3 = INE(3,IE) LAMBDA(1) = ONESIXTH * (C(1,I1)+C(1,I2)+C(1,I3)) LAMBDA(2) = ONESIXTH * (C(2,I1)+C(2,I2)+C(2,I3)) K(1) = LAMBDA(1) * IEN(1,IE) + LAMBDA(2) * IEN(2,IE) K(2) = LAMBDA(1) * IEN(3,IE) + LAMBDA(2) * IEN(4,IE) K(3) = LAMBDA(1) * IEN(5,IE) + LAMBDA(2) * IEN(6,IE) KP(1,IE) = MAX(ZERO,K(1)) KP(2,IE) = MAX(ZERO,K(2)) KP(3,IE) = MAX(ZERO,K(3)) KM(1) = MIN(ZERO,K(1)) KM(2) = MIN(ZERO,K(2)) KM(3) = MIN(ZERO,K(3)) FL11 = C(1,I2)*IEN(1,IE)+C(2,I2)*IEN(2,IE) FL12 = C(1,I3)*IEN(1,IE)+C(2,I3)*IEN(2,IE) FL21 = C(1,I3)*IEN(3,IE)+C(2,I3)*IEN(4,IE) FL22 = C(1,I1)*IEN(3,IE)+C(2,I1)*IEN(4,IE) FL31 = C(1,I1)*IEN(5,IE)+C(2,I1)*IEN(6,IE) FL32 = C(1,I2)*IEN(5,IE)+C(2,I2)*IEN(6,IE) CRFS(1) = - ONESIXTH * (TWO *FL31 + FL32 + FL21 + TWO * FL22 ) CRFS(2) = - ONESIXTH * (TWO *FL32 + TWO * FL11 + FL12 + FL31 ) CRFS(3) = - ONESIXTH * (TWO *FL12 + TWO * FL21 + FL22 + FL11 ) DELTAL(:,IE) = CRFS(:)-KP(:,IE) NM(IE) = ONE/MIN(-THR,SUM(KM(:))) END DO U(:) = AC2(IS,ID,:) J = 0 ! Counter ... ASPAR = ZERO ! Mass matrix ... B = ZERO ! Right hand side ... ! ! ... assembling the linear equation system .... ! DO IP = 1, MNP DO I = 1, CCON(IP) J = J + 1 IE = IE_CELL(J) POS = POS_CELL(J) K1 = KP(POS,IE) ! Flux Jacobian TRIA03 = ONETHIRD * TRIA(IE) DTK = K1 * DT4A * IOBPD(ID,IP) * IOBWB(IP) * IOBDP(IP) TMP3 = DTK * NM(IE) I1 = POSI(1,J) ! Position of the recent entry in the ASPAR matrix ... ASPAR is shown in fig. 42, p.122 I2 = POSI(2,J) I3 = POSI(3,J) ASPAR(I1) = TRIA03 + DTK - TMP3 * DELTAL(POS ,IE) + ASPAR(I1) ! Diagonal entry ASPAR(I2) = - TMP3 * DELTAL(POS_TRICK(POS,1),IE) + ASPAR(I2) ! off diagonal entries ... ASPAR(I3) = - TMP3 * DELTAL(POS_TRICK(POS,2),IE) + ASPAR(I3) B(IP) = B(IP) + TRIA03 * U(IP) END DO !I: loop over connected elements ... END DO !IP IF (LBCWA .OR. LBCSP) THEN IF (LINHOM) THEN DO IP = 1, IWBMNP IPGL1 = IWBNDLC(IP) ASPAR(I_DIAG(IPGL1)) = SI(IPGL1) ! Set boundary on the diagonal B(IPGL1) = SI(IPGL1) * WBAC(IS,ID,IP) END DO ELSE DO IP = 1, IWBMNP IPGL1 = IWBNDLC(IP) ASPAR(I_DIAG(IPGL1)) = SI(IPGL1) ! Set boundary on the diagonal B(IPGL1) = SI(IPGL1) * WBAC(IS,ID,1) END DO ENDIF END IF IF (ICOMP .GE. 2 .AND. SMETHOD .GT. 0 .AND. LSOURCESWAM) THEN DO IP = 1, MNP IF (IOBWB(IP) .EQ. 1) THEN ! GTEMP1 = MAX((1.-DT4A*IMATDAA(IS,ID,IP)),1.) GTEMP2 = IMATRAA(IS,ID,IP)/MAX((ONE-DT4A*IMATDAA(IS,ID,IP)),ONE) DELFL = COFRM4(IS)*DT4S USFM = USNEW(IP)*MAX(FMEANWS(IP),FMEAN(IP)) FLHAB = ABS(GTEMP2*DT4S) FLHAB = MIN(FLHAB,USFM*DELFL)/DT4S B(IP) = B(IP)+SIGN(FLHAB,GTEMP2)*DT4S*SI(IP) !LIMFAC = MIN(ONE,ABS(SIGN(FLHAB,GTEMP2))/MAX(THR,ABS(IMATRAA(IP,IS,ID)))) !ASPAR(I_DIAG(IP)) = ASPAR(I_DIAG(IP))-DT4A*LIMFAC*IMATDAA(IP,IS,ID) ENDIF END DO ELSE IF (ICOMP .GE. 2 .AND. SMETHOD .GT. 0 .AND. .NOT. LSOURCESWAM) THEN DO IP = 1, MNP IF (IOBWB(IP) .EQ. 1) THEN ASPAR(I_DIAG(IP)) = ASPAR(I_DIAG(IP)) + IMATDAA(IS,ID,IP) * DT4A * SI(IP) ! Add source term to the diagonal B(IP) = B(IP) + IMATRAA(IS,ID,IP) * DT4A * SI(IP) ! Add source term to the right hand side ENDIF END DO ENDIF ! IPAR(1) = 0 ! always 0 to start an iterative solver IPAR(2) = 1 ! right preconditioning IPAR(3) = 1 ! use convergence test scheme 1 IPAR(4) = 200*MNP ! IPAR(5) = 15 IPAR(6) = 1000 ! use at most 1000 matvec's FPAR(1) = 1.0E-10 ! relative tolerance 1.0E-6 FPAR(2) = 1.0E-12 ! absolute tolerance 1.0E-10 FPAR(11) = ZERO ! clearing the FLOPS counter AU = 0. FLJAU = 0. FLJU = 0. CALL ILU0 (MNP, ASPAR, JA, IA, AU, FLJAU, FLJU, IWKSP, IERROR) ! WRITE(DBG%FHNDL,*) 'CALL SOLVER' ! WRITE(DBG%FHNDL,*) DT4A, MSC, MDC, MNE ! WRITE(DBG%FHNDL,*) 'WRITE CG', SUM(CG) ! WRITE(DBG%FHNDL,*) SUM(XP), SUM(YP) ! WRITE(DBG%FHNDL,*) SUM(IMATRAA), SUM(IMATDAA) ! WRITE(DBG%FHNDL,*) SUM(B), SUM(X) ! WRITE(DBG%FHNDL,*) SUM(IPAR), SUM(FPAR) ! WRITE(DBG%FHNDL,*) SUM(WKSP), SUM(INIU) ! WRITE(DBG%FHNDL,*) SUM(ASPAR), SUM(IA), SUM(JA) ! WRITE(DBG%FHNDL,*) SUM(AU), SUM(FLJAU), SUM(FLJU) INIU = AC2(IS,ID,:) * IOBPD(ID,:) X = ZERO CALL RUNRC (MNP, NNZ, B, X, IPAR, FPAR, WKSP, INIU, ASPAR, JA, IA, AU, FLJAU, FLJU, BCGSTAB) DO IP = 1, MNP AC2(IS,ID,IP) = MAX(ZERO,X(IP)) !* MyREAL(IOBPD(ID,IP)) END DO #ifdef TIMINGS CALL WAV_MY_WTIME(TIME4) #endif ! WRITE(DBG%FHNDL,*) 'SOLUTION' ! WRITE(DBG%FHNDL,*) SUM(B), SUM(X) ! WRITE(DBG%FHNDL,*) SUM(IPAR), SUM(FPAR) ! WRITE(DBG%FHNDL,*) SUM(WKSP), SUM(INIU) ! WRITE(DBG%FHNDL,*) SUM(ASPAR), SUM(JA), SUM(JA) ! WRITE(DBG%FHNDL,*) SUM(AU), SUM(FLJAU), SUM(FLJU) IF (LADVTEST) THEN WRITE(4001) SNGL(RTIME) WRITE(4001) (SNGL(C(1,IP)), SNGL(C(2,IP)), SNGL(AC2(1,1,IP)), & & IP = 1, MNP) CALL CHECKCONS(U,SUMAC2) IF (MINVAL(U) .LT. MINTEST) MINTEST = MINVAL(U) WRITE (*,*) 'VOLUMES AT T0, T1 and T2',SUMACt0, & & SUMAC1, SUMAC2, MINTEST WRITE (*,*) 'VOLUME ERROR: TOTAL and ACTUAL', & & 100.0_rkind-((SUMACt0-SUMAC2)/SUMACt0)*100.0_rkind, & & 100.0_rkind- ((SUMAC1-SUMAC2)/SUMAC1)*100.0_rkind END IF END SUBROUTINE !********************************************************************** !* !********************************************************************** SUBROUTINE EIMPS_ASPAR_B_SOURCES_LOCAL( IS, ID, ASPAR, B, U) USE DATAPOOL IMPLICIT NONE INTEGER, INTENT(IN) :: IS,ID REAL(rkind), intent(inout) :: ASPAR(NNZ) REAL(rkind), intent(inout) :: B(MNP) REAL(rkind), intent(inout) :: U(MNP) INTEGER :: POS_TRICK(3,2) REAL(rkind) :: FL11, FL12, FL21, FL22, FL31, FL32 REAL(rkind):: CRFS(3), K1, KM(3), K(3), TRIA03 REAL(rkind) :: C(2,MNP) REAL(rkind) :: DELTAL(3,MNE) REAL(rkind) :: GTEMP1, GTEMP2, DELT, XIMP, DELT5, USFM REAL(rkind) :: FLHAB, DELFL, TEMP INTEGER :: I1, I2, I3 INTEGER :: IP, IE, POS INTEGER :: I, J, IPGL1, IPrel REAL(rkind) :: KP(3,MNE), NM(MNE) REAL(rkind) :: DTK, TMP3 REAL(rkind) :: LAMBDA(2) POS_TRICK(1,1) = 2 POS_TRICK(1,2) = 3 POS_TRICK(2,1) = 3 POS_TRICK(2,2) = 1 POS_TRICK(3,1) = 1 POS_TRICK(3,2) = 2 CALL CADVXY(IS,ID,C) ! ! Calculate countour integral quantities ... ! DO IE = 1, MNE I1 = INE(1,IE) I2 = INE(2,IE) I3 = INE(3,IE) LAMBDA(1) = ONESIXTH * (C(1,I1)+C(1,I2)+C(1,I3)) LAMBDA(2) = ONESIXTH * (C(2,I1)+C(2,I2)+C(2,I3)) K(1) = LAMBDA(1) * IEN(1,IE) + LAMBDA(2) * IEN(2,IE) K(2) = LAMBDA(1) * IEN(3,IE) + LAMBDA(2) * IEN(4,IE) K(3) = LAMBDA(1) * IEN(5,IE) + LAMBDA(2) * IEN(6,IE) KP(1,IE) = MAX(ZERO,K(1)) KP(2,IE) = MAX(ZERO,K(2)) KP(3,IE) = MAX(ZERO,K(3)) KM(1) = MIN(ZERO,K(1)) KM(2) = MIN(ZERO,K(2)) KM(3) = MIN(ZERO,K(3)) FL11 = C(1,I2)*IEN(1,IE)+C(2,I2)*IEN(2,IE) FL12 = C(1,I3)*IEN(1,IE)+C(2,I3)*IEN(2,IE) FL21 = C(1,I3)*IEN(3,IE)+C(2,I3)*IEN(4,IE) FL22 = C(1,I1)*IEN(3,IE)+C(2,I1)*IEN(4,IE) FL31 = C(1,I1)*IEN(5,IE)+C(2,I1)*IEN(6,IE) FL32 = C(1,I2)*IEN(5,IE)+C(2,I2)*IEN(6,IE) CRFS(1) = - ONESIXTH * (TWO *FL31 + FL32 + FL21 + TWO * FL22 ) CRFS(2) = - ONESIXTH * (TWO *FL32 + TWO * FL11 + FL12 + FL31 ) CRFS(3) = - ONESIXTH * (TWO *FL12 + TWO * FL21 + FL22 + FL11 ) DELTAL(:,IE) = CRFS(:)- KP(:,IE) NM(IE) = ONE/MIN(-THR,SUM(KM(:))) END DO J = 0 ! Counter ... ASPAR = ZERO ! Mass matrix ... B = ZERO ! Right hand side ... ! ! ... assembling the linear equation system .... ! DO IP = 1, NP_RES IF (IOBPD(ID,IP) .EQ. 1 .AND. IOBWB(IP) .EQ. 1 .AND. DEP(IP) .GT. DMIN) THEN DO I = 1, CCON(IP) J = J + 1 IE = IE_CELL(J) POS = POS_CELL(J) K1 = KP(POS,IE) ! Flux Jacobian TRIA03 = ONETHIRD * TRIA(IE) DTK = K1 * DT4A * IOBPD(ID,IP) * IOBWB(IP) * IOBDP(IP) TMP3 = DTK * NM(IE) I1 = POSI(1,J) ! Position of the recent entry in the ASPAR matrix ... ASPAR is shown in fig. 42, p.122 I2 = POSI(2,J) I3 = POSI(3,J) ASPAR(I1) = TRIA03 + DTK - TMP3 * DELTAL(POS ,IE) + ASPAR(I1) ! Diagonal entry ASPAR(I2) = - TMP3 * DELTAL(POS_TRICK(POS,1),IE) + ASPAR(I2) ! off diagonal entries ... ASPAR(I3) = - TMP3 * DELTAL(POS_TRICK(POS,2),IE) + ASPAR(I3) B(IP) = B(IP) + TRIA03 * U(IP) END DO !I: loop over connected elements ... ELSE DO I = 1, CCON(IP) J = J + 1 IE = IE_CELL(J) TRIA03 = ONETHIRD * TRIA(IE) I1 = POSI(1,J) ! Position of the recent entry in the ASPAR matrix ... ASPAR is shown in fig. 42, p.122 ASPAR(I1) = TRIA03 + ASPAR(I1) ! Diagonal entry B(IP) = 0.!B(IP) + TRIA03 * 0. END DO !I: loop over connected elements ... END IF END DO !IP #if defined DEBUG WRITE(3000+myrank,*) 'IS, ID, sum=', IS, ID, sum(ASPAR) #endif IF (LBCWA .OR. LBCSP) THEN IF (LINHOM) THEN IPrel=IP ELSE IPrel=1 ENDIF DO IP = 1, IWBMNP IPGL1 = IWBNDLC(IP) ASPAR(I_DIAG(IPGL1)) = SI(IPGL1) ! Set boundary on the diagonal B(IPGL1) = SI(IPGL1) * WBAC(IS,ID,IPrel) END DO END IF IF (ICOMP .GE. 2 .AND. SMETHOD .GT. 0 .AND. LSOURCESWAM) THEN DO IP = 1, MNP IF (IOBWB(IP) .EQ. 1) THEN !GTEMP1 = MAX((1.-DT4A*FL(IP,ID,IS)),1.) !GTEMP2 = SL(IP,ID,IS)/GTEMP1/PI2/SPSIG(IS) GTEMP1 = MAX((ONE-DT4A*IMATDAA(IS,ID,IP)),ONE) GTEMP2 = IMATRAA(IP,IS,ID)/GTEMP1!/PI2/SPSIG(IS) DELT = DT4S XIMP = 1.0 DELT5 = XIMP*DELT DELFL = COFRM4(IS)*DELT USFM = USNEW(IP)*MAX(FMEANWS(IP),FMEAN(IP)) TEMP = USFM*DELFL!/PI2/SPSIG(IS) FLHAB = ABS(GTEMP2*DT4S) FLHAB = MIN(FLHAB,TEMP)/DT4S B(IP) = B(IP) + SIGN(FLHAB,GTEMP2) * DT4A * SI(IP) ! Add source term to the right hand side ASPAR(I_DIAG(IP)) = ASPAR(I_DIAG(IP)) - SIGN(FLHAB,GTEMP2) * SI(IP) !!B(IP) = B(IP) + GTEMP2 * DT4A * SI(IP) ! Add source term to the right hand side !ASPAR(I_DIAG(IP)) = ASPAR(I_DIAG(IP)) - GTEMP2 * SI(IP) !This is then for the shallow water physics take care about ISELECT !ASPAR(I_DIAG(IP)) = ASPAR(I_DIAG(IP)) + IMATDAA(IS,ID,IP) * DT4A * SI(IP) ! Add source term to the diagonal !B(IP) = B(IP) + IMATRAA(IP,IS,ID) * DT4A * SI(IP) ! Add source term to the right hand side ENDIF END DO ELSE IF (ICOMP .GE. 2 .AND. SMETHOD .GT. 0 .AND. .NOT. LSOURCESWAM) THEN DO IP = 1, MNP IF (IOBWB(IP) .EQ. 1) THEN ASPAR(I_DIAG(IP)) = ASPAR(I_DIAG(IP)) + IMATDAA(IS,ID,IP) * DT4A * SI(IP) ! Add source term to the diagonal B(IP) = B(IP) + IMATRAA(IS,ID,IP) * DT4A * SI(IP) ! Add source term to the right hand side ENDIF END DO ENDIF END SUBROUTINE !********************************************************************** !* !********************************************************************** SUBROUTINE EIMPS( IS, ID) USE DATAPOOL IMPLICIT NONE INTEGER, INTENT(IN) :: IS,ID REAL(rkind) :: ASPAR(NNZ) REAL(rkind) :: X(MNP) REAL(rkind) :: B(MNP) REAL(rkind) :: U(MNP) INTEGER :: IPAR(16) INTEGER :: IERROR INTEGER :: IWKSP(20*MNP) INTEGER :: FLJU(MNP) INTEGER :: FLJAU(NNZ+1) REAL(rkind) :: FPAR(16) REAL(rkind) :: WKSP( 20*MNP ) REAL(rkind) :: AU(NNZ+1) REAL(rkind) :: INIU(MNP) REAL(rkind) :: TIME1, TIME2, TIME3, TIME4 INTEGER :: IP external bcgstab external gmres IWKSP = 0 WKSP = ZERO U(:) = AC2(IS,ID,:) CALL EIMPS_ASPAR_B_SOURCES_LOCAL( IS, ID, ASPAR, B, U) ! IPAR(1) = 0 ! always 0 to start an iterative solver IPAR(2) = 1 ! right preconditioning IPAR(3) = 1 ! use convergence test scheme 1 IPAR(4) = 200*MNP ! IPAR(5) = 15 IPAR(6) = 1000 ! use at most 1000 matvec's FPAR(1) = 1.0E-10 ! relative tolerance 1.0E-6 FPAR(2) = 1.0E-12 ! absolute tolerance 1.0E-10 FPAR(11) = ZERO ! clearing the FLOPS counter AU = 0. FLJAU = 0. FLJU = 0. ! CALL ILU0(MNP, ASPAR, JA, IA, AU, FLJAU, FLJU, IWKSP, IERROR) CALL SOR(MNP, ASPAR, JA, IA, AU, FLJAU, FLJU, IWKSP, IERROR) ! WRITE(DBG%FHNDL,*) 'CALL SOLVER' ! WRITE(DBG%FHNDL,*) DT4A, MSC, MDC, MNE ! WRITE(DBG%FHNDL,*) 'WRITE CG', SUM(CG) ! WRITE(DBG%FHNDL,*) SUM(XP), SUM(YP) ! WRITE(DBG%FHNDL,*) SUM(IMATRAA), SUM(IMATDAA) ! WRITE(DBG%FHNDL,*) SUM(B), SUM(X) ! WRITE(DBG%FHNDL,*) SUM(IPAR), SUM(FPAR) ! WRITE(DBG%FHNDL,*) SUM(WKSP), SUM(INIU) ! WRITE(DBG%FHNDL,*) SUM(ASPAR), SUM(IA), SUM(JA) ! WRITE(DBG%FHNDL,*) SUM(AU), SUM(FLJAU), SUM(FLJU) INIU = AC2(IS,ID,:) * IOBPD(ID,:) X = ZERO CALL RUNRC (MNP, NNZ, B, X, IPAR, FPAR, WKSP, INIU, ASPAR, JA, IA, AU, FLJAU, FLJU, BCGSTAB) DO IP = 1, MNP AC2(IS,ID,IP) = MAX(ZERO,X(IP)) * MyREAL(IOBPD(ID,IP)) END DO #ifdef TIMINGS ! CALL WAV_MY_WTIME(TIME4) #endif ! WRITE(DBG%FHNDL,*) 'SOLUTION' ! WRITE(DBG%FHNDL,*) SUM(B), SUM(X) ! WRITE(DBG%FHNDL,*) SUM(IPAR), SUM(FPAR) ! WRITE(DBG%FHNDL,*) SUM(WKSP), SUM(INIU) ! WRITE(DBG%FHNDL,*) SUM(ASPAR), SUM(JA), SUM(JA) ! WRITE(DBG%FHNDL,*) SUM(AU), SUM(FLJAU), SUM(FLJU) IF (LADVTEST) THEN WRITE(4001) SNGL(RTIME) ! WRITE(4001) (SNGL(C(1,IP)), SNGL(C(2,IP)), SNGL(AC2(1,1,IP)), & ! & IP = 1, MNP) CALL CHECKCONS(U,SUMAC2) IF (MINVAL(U) .LT. MINTEST) MINTEST = MINVAL(U) WRITE (*,*) 'VOLUMES AT T0, T1 and T2',SUMACt0, & & SUMAC1, SUMAC2, MINTEST WRITE (*,*) 'VOLUME ERROR: TOTAL and ACTUAL', & & 100.0_rkind-((SUMACt0-SUMAC2)/SUMACt0)*100.0_rkind, & & 100.0_rkind- ((SUMAC1-SUMAC2)/SUMAC1)*100.0_rkind END IF END SUBROUTINE !********************************************************************** !* !********************************************************************** SUBROUTINE SQUARE_NORM(eV1, eV2, eScal) USE DATAPOOL, only : rkind, MNP, MDC, NP_RES #ifdef MPI_PARALL_GRID USE DATAPOOL, only : nwild_loc_res USE datapool, only : myrank, comm, ierr, nproc, istatus, rtype #endif implicit none real(rkind), intent(in) :: eV1(MNP), eV2(MNP) real(rkind), intent(out) :: eScal real(rkind) :: LScal(1), RScal(1) integer IP, iProc #ifndef MPI_PARALL_GRID eScal=0 DO IP=1,MNP eScal=eScal + (eV1(IP) - eV2(IP) )**2 END DO #else eScal=0 DO IP=1,NP_RES eScal=eScal + nwild_loc_res(IP)*(eV1(IP) - eV2(IP) )**2 END DO LScal(1)=eScal IF (myrank == 0) THEN DO iProc=2,nproc CALL MPI_RECV(RScal,1,rtype, iProc-1, 19, comm, istatus, ierr) LScal = LScal + RScal END DO DO iProc=2,nproc CALL MPI_SEND(LScal,1,rtype, iProc-1, 23, comm, ierr) END DO ELSE CALL MPI_SEND(LScal,1,rtype, 0, 19, comm, ierr) CALL MPI_RECV(LScal,1,rtype, 0, 23, comm, istatus, ierr) END IF eScal=LScal(1) #endif END SUBROUTINE !********************************************************************** !* !********************************************************************** SUBROUTINE EIMPS_JACOBI_ITERATION( IS, ID) USE DATAPOOL IMPLICIT NONE INTEGER, INTENT(IN) :: IS,ID REAL(rkind) :: ASPAR(NNZ) REAL(rkind) :: X(MNP), B(MNP), U(MNP) REAL(rkind) :: eSum, eSqrNorm INTEGER :: IP, idx, J, nbIter REAL(rkind) :: LOCAL_SOLVERTHR X(:) = AC2(IS,ID,:) CALL EIMPS_ASPAR_B_SOURCES_LOCAL( IS, ID, ASPAR, B, X) LOCAL_SOLVERTHR=10E-10 nbIter=0 DO DO IP=1,NP_RES eSum=B(IP) DO J=IA(IP),IA(IP+1)-1 IF (J .ne. I_DIAG(IP)) THEN idx=JA(J) eSum=eSum - ASPAR(J)*X(idx) END IF END DO eSum=eSum/ASPAR(I_DIAG(IP)) U(IP)=eSum END DO #ifdef MPI_PARALL_GRID CALL EXCHANGE_P2D(U) #endif X=U DO IP=1,NP_RES eSum=0 DO J=IA(IP),IA(IP+1)-1 idx=JA(J) eSum=eSum + ASPAR(J)*X(idx) END DO U(IP)=eSum END DO #ifdef MPI_PARALL_GRID CALL EXCHANGE_P2D(U) #endif CALL SQUARE_NORM(U, B, eSqrNorm) nbIter=nbIter+1 IF (eSqrNorm .lt. LOCAL_SOLVERTHR) THEN EXIT END IF END DO DO IP = 1, MNP AC2(IS,ID,IP) = MAX(ZERO,X(IP)) * MyREAL(IOBPD(ID,IP)) END DO END SUBROUTINE !********************************************************************** !* ZYL: Crank-Nicolson implicit !********************************************************************** SUBROUTINE CNIMPS_ASPAR_B(IS, ID, ASPAR, B, U) USE DATAPOOL implicit none INTEGER, INTENT(IN) :: IS,ID REAL(rkind), intent(out) :: ASPAR(NNZ) REAL(rkind), intent(out) :: B(MNP) REAL(rkind), intent(in) :: U(MNP) INTEGER :: IP, IE, POS INTEGER :: I1, I2, I3 INTEGER :: I, J INTEGER :: IP_J, IP_K !, LFIL INTEGER :: POS_TRICK(3,2) REAL(rkind) :: LAMBDA(2), FL11, FL12, FL21, FL22, FL31, FL32 REAL(rkind) :: CRFS(3), K1 REAL(rkind) :: IEN1(2), IEN2(2), IEN3(2), NM(MNE) REAL(rkind):: DELTAL(3,MNE), C(2,MNP), K(3), KP(3,MNE), KM(3,MNE) REAL(rkind) :: DT05 REAL(rkind) :: TMP1, TMP2, TMP3, TMP5, TMP6, TMP7, TMP8, BT1 POS_TRICK(1,1) = 2 POS_TRICK(1,2) = 3 POS_TRICK(2,1) = 3 POS_TRICK(2,2) = 1 POS_TRICK(3,1) = 1 POS_TRICK(3,2) = 2 CALL CADVXY(IS,ID,C) DT05 = 0.5 * DT4A DO IE = 1, MNE I1 = INE(1,IE) I2 = INE(2,IE) I3 = INE(3,IE) IEN1 = IEN(1:2,IE) IEN2 = IEN(3:4,IE) IEN3 = IEN(5:6,IE) LAMBDA(1) = ONESIXTH * (C(1,I1)+C(1,I2)+C(1,I3)) LAMBDA(2) = ONESIXTH * (C(2,I1)+C(2,I2)+C(2,I3)) K(1) = LAMBDA(1) * IEN(1,IE) + LAMBDA(2) * IEN(2,IE) K(2) = LAMBDA(1) * IEN(3,IE) + LAMBDA(2) * IEN(4,IE) K(3) = LAMBDA(1) * IEN(5,IE) + LAMBDA(2) * IEN(6,IE) KP(1,IE) = MAX(ZERO,K(1)) KP(2,IE) = MAX(ZERO,K(2)) KP(3,IE) = MAX(ZERO,K(3)) KM(1,IE) = MIN(ZERO,K(1)) KM(2,IE) = MIN(ZERO,K(2)) KM(3,IE) = MIN(ZERO,K(3)) FL11 = C(1,I2)*IEN1(1)+C(2,I2)*IEN1(2) FL12 = C(1,I3)*IEN1(1)+C(2,I3)*IEN1(2) FL21 = C(1,I3)*IEN2(1)+C(2,I3)*IEN2(2) FL22 = C(1,I1)*IEN2(1)+C(2,I1)*IEN2(2) FL31 = C(1,I1)*IEN3(1)+C(2,I1)*IEN3(2) FL32 = C(1,I2)*IEN3(1)+C(2,I2)*IEN3(2) CRFS(1) = - ONESIXTH * (2. *FL31 + FL32 + FL21 + 2. * FL22 ) CRFS(2) = - ONESIXTH * (2. *FL32 + 2. * FL11 + FL12 + FL31 ) CRFS(3) = - ONESIXTH * (2. *FL12 + 2. * FL21 + FL22 + FL11 ) DELTAL(:,IE) = CRFS(:) - MAX(K(:),ZERO) NM(IE) = ONE/MIN(-THR,SUM(KM(:,IE))) END DO B = ZERO ASPAR = ZERO J = 0 DO IP = 1, MNP DO I = 1, CCON(IP) J = J + 1 IE = IE_CELL(J) POS = POS_CELL(J) K1 = KP(POS,IE) IP_J = INE(POS_TRICK(POS,1),IE) IP_K = INE(POS_TRICK(POS,2),IE) TMP3 = NM(IE) TMP2 = ONETHIRD * TRIA(IE) TMP1 = DT05 * K1 * TMP3 TMP5 = DT05 * K1 TMP6 = TMP1 * DELTAL(POS,IE) TMP7 = TMP1 * DELTAL(POS_TRICK(POS,1),IE) TMP8 = TMP1 * DELTAL(POS_TRICK(POS,2),IE) BT1 = TMP2 - TMP5 + TMP6 B(IP) = ( BT1 * U(IP) + TMP7 * U(IP_J) + TMP8 * U(IP_K) ) + B(IP) I1 = POSI(1,J) I2 = POSI(2,J) I3 = POSI(3,J) ASPAR(I1) = TMP2 + TMP5 - TMP6 + ASPAR(I1) ASPAR(I2) = - TMP7 + ASPAR(I2) ASPAR(I3) = - TMP8 + ASPAR(I3) END DO END DO IF (ICOMP .GE. 2 .AND. SMETHOD .GT. 0) THEN DO IP = 1, MNP ASPAR(I_DIAG(IP)) = ASPAR(I_DIAG(IP)) + IMATDAA(IS,ID,IP) * DT4A * SI(IP) ! Add source term to the diagonal B(IP) = B(IP) + IMATRAA(IS,ID,IP) * DT4A * SI(IP) ! Add source term to the right hand side END DO END IF END SUBROUTINE !********************************************************************** !* ZYL: Crank-Nicolson implicit !********************************************************************** SUBROUTINE CNIMPS(IS, ID) USE DATAPOOL IMPLICIT NONE INTEGER, INTENT(IN) :: IS,ID INTEGER :: IPAR(16) INTEGER :: IERROR ! IWK ! ERROR Indicator and Work Array Size, INTEGER :: IP INTEGER :: IWKSP( 8*MNP ) ! Integer Workspace INTEGER :: JAU(NNZ+1), JU(MNP) REAL(rkind) :: FPAR(16) ! DROPTOL REAL(rkind) :: WKSP( 8 * MNP ) ! REAL WORKSPACES REAL(rkind) :: AU(NNZ+1), ASPAR(NNZ) REAL(rkind) :: U(MNP), X(MNP), B(MNP) REAL(rkind) :: INIU(MNP) external bcgstab external gmres U(:) = AC2(IS,ID,:) CALL CNIMPS_ASPAR_B(IS, ID, ASPAR, B, U) IPAR(1) = 0 ! always 0 to start an iterative solver IPAR(2) = 1 ! right preconditioning IPAR(3) = 1 ! use convergence test scheme 1 IPAR(4) = 8*MNP ! the 'w' has 10,000 elements IPAR(5) = 30 IPAR(6) = 100 ! use at most 100 matvec's FPAR(1) = 1.0d-8 ! relative tolerance 1.0E-6 FPAR(2) = 1.0d-10 ! absolute tolerance 1.0E-10BCGSTAB FPAR(11) = 0.0 ! clearing the FLOPS counter ! IWK = NNZ+1 ! DROPTOL = VERYSMALL ! LFIL = MNP - 20 JU = 0 IWKSP = 0 WKSP = ZERO INIU(:) = AC2(IS,ID,:) ! CALL MILU0 (MNP, ASPAR, JA, IA, AU, JAU, JU, IWKSP, IERROR) CALL ILU0 (MNP, ASPAR, JA, IA, AU, JAU, JU, IWKSP, IERROR) ! CALL ILUT (MNP, ASPAR, JA, IA, LFIL, DROPTOL, AU, JAU, JU, IWK, WKSP, IWKSP, IERROR) !... O.K ! CALL ILUK (MNP, ASPAR, JA, IA, LFIL, AU, JAU, JU, LEVS, IWK, WKSP, IWKSP, IERROR) CALL RUNRC(MNP, NNZ, B, X, IPAR, FPAR, WKSP, INIU, ASPAR, JA, IA, AU, JAU, JU, BCGSTAB) DO IP = 1, MNP AC2(IS,ID,IP) = MAX(ZERO,X(IP)) * IOBPD(ID,IP) END DO IF (LADVTEST) THEN WRITE(4001) SNGL(RTIME) ! WRITE(4001) (SNGL(C(1,IP)), SNGL(C(2,IP)), SNGL(AC2(1,1,IP)),& ! & IP = 1, MNP) CALL CHECKCONS(U,SUMAC2) IF (MINVAL(U) .LT. MINTEST) MINTEST = MINVAL(U) WRITE (*,*) 'VOLUMES AT T0, T1 and T2',SUMACt0, & & SUMAC1, SUMAC2, MINTEST WRITE (*,*) 'VOLUME ERROR: TOTAL and ACTUAL', & & 100.0_rkind-((SUMACt0-SUMAC2)/SUMACt0)*100.0_rkind, & & 100.0_rkind- ((SUMAC1-SUMAC2)/SUMAC1)*100.0_rkind END IF END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE CNEIMPS(IS,ID,DTMAX) USE DATAPOOL IMPLICIT NONE INTEGER, INTENT(IN) :: IS, ID REAL(rkind), INTENT(IN) :: DTMAX INTEGER :: I, J INTEGER :: IP, IE, POS INTEGER :: I1, I2, I3, IP_J, IP_K REAL(rkind) :: N(MNE), DT1, DT2, DT05, DT105, DT205 REAL(rkind) :: TMP1, TMP2, TMP3, TRIA03 REAL(rkind) :: LAMBDA(2), BTMP(3) REAL(rkind) :: K(3,MNE), KP(3,MNE), KM(3,MNE), KMAX REAL(rkind) :: DELTA(3,MNE), CRFS(3), KPN(3,MNE), KPNTMP REAL(rkind) :: U(MNP), X(MNP), ASPAR1(NNZ), B1(MNP), ASPAR2(NNZ), B2(MNP) REAL(rkind) :: FL11, FL12, FL21, FL22, FL31, FL32, C(2,MNP) INTEGER, PARAMETER :: IM = 30 INTEGER :: IPAR(16), MBLOC INTEGER :: IWKSP( 8*MNP ) ! Integer Workspace INTEGER :: JU(MNP), JAU(NNZ+1), IWK, LFIL INTEGER :: IERROR REAL(rkind) :: FPAR(16), DROPTOL, PERMTOL REAL(rkind) :: WKSP(8*MNP ) ! REAL WORKSPACES REAL(rkind) :: AU(NNZ+1) REAL(rkind) :: INIU(MNP) INTEGER :: POS_TRICK(3,2) external cgnr external bcg external bcgstab external tfqmr external fom external gmres external dqgmres external fgmres external dbcg POS_TRICK(1,1) = 2 POS_TRICK(1,2) = 3 POS_TRICK(2,1) = 3 POS_TRICK(2,2) = 1 POS_TRICK(3,1) = 1 POS_TRICK(3,2) = 2 IWKSP = 0 WKSP = ZERO INIU = ZERO U(:) = AC2(IS,ID,:) CALL CADVXY(IS,ID,C) DT1 = MIN(DT4A, DTMAX) DT2 = MAX(THR, DT4A - DT1) DT05 = 0.5 * DT4A DT105 = 0.5 * DT1 DT205 = 0.5 * DT2 DO IE = 1, MNE I1 = INE(1,IE) I2 = INE(2,IE) I3 = INE(3,IE) LAMBDA(1) = ONESIXTH * (C(1,I1)+C(1,I2)+C(1,I3)) LAMBDA(2) = ONESIXTH * (C(2,I1)+C(2,I2)+C(2,I3)) K(1,IE) = LAMBDA(1) * IEN(1,IE) + LAMBDA(2) * IEN(2,IE) K(2,IE) = LAMBDA(1) * IEN(3,IE) + LAMBDA(2) * IEN(4,IE) K(3,IE) = LAMBDA(1) * IEN(5,IE) + LAMBDA(2) * IEN(6,IE) KP(:,IE) = MAX(K(:,IE),ZERO) KM(:,IE) = MIN(K(:,IE),ZERO) N(IE) = - ONE/MIN(-THR,SUM(KM(:,IE))) KPN(:,IE) = KP(:,IE) * N(IE) FL11 = C(1,I2)*IEN(1,IE)+C(2,I2)*IEN(2,IE) FL12 = C(1,I3)*IEN(1,IE)+C(2,I3)*IEN(2,IE) FL21 = C(1,I3)*IEN(3,IE)+C(2,I3)*IEN(4,IE) FL22 = C(1,I1)*IEN(3,IE)+C(2,I1)*IEN(4,IE) FL31 = C(1,I1)*IEN(5,IE)+C(2,I1)*IEN(6,IE) FL32 = C(1,I2)*IEN(5,IE)+C(2,I2)*IEN(6,IE) CRFS(1) = - ONESIXTH * (2. *FL31 + FL32 + FL21 + 2. * FL22 ) CRFS(2) = - ONESIXTH * (2. *FL32 + 2. * FL11 + FL12 + FL31 ) CRFS(3) = - ONESIXTH * (2. *FL12 + 2. * FL21 + FL22 + FL11 ) DELTA(:,IE) = CRFS(:) - KP(:,IE) END DO B1(:) = 0. ASPAR1(:) = 0. B2(:) = 0. ASPAR2(:) = 0. J = 0 DO IP = 1, MNP KMAX = 0. DO I = 1, CCON(IP) J = J + 1 IE = IE_CELL(J) POS = POS_CELL(J) IP_J = INE(POS_TRICK(POS,1),IE) IP_K = INE(POS_TRICK(POS,2),IE) I1 = POSI(1,J) I2 = POSI(2,J) I3 = POSI(3,J) KPNTMP = KPN(POS,IE) TRIA03 = ONETHIRD * TRIA(IE) TMP1 = DT05 * KPNTMP TMP2 = DT105 * KPNTMP TMP3 = DT205 * KPNTMP ASPAR1(I1) = TRIA03 + DT05 * KP(POS,IE) - TMP1 * DELTA(POS,IE) + ASPAR1(I1) ASPAR1(I2) = - TMP1 * DELTA(POS_TRICK(POS,1),IE) + ASPAR1(I2) ASPAR1(I3) = - TMP1 * DELTA(POS_TRICK(POS,2),IE) + ASPAR1(I3) ASPAR2(I1) = TRIA03 + DT205 * KP(POS,IE) - TMP3 * DELTA(POS,IE) + ASPAR2(I1) ASPAR2(I2) = - TMP3 * DELTA(POS_TRICK(POS,1),IE) + ASPAR2(I2) ASPAR2(I3) = - TMP3 * DELTA(POS_TRICK(POS,2),IE) + ASPAR2(I3) BTMP(1) = TRIA03 - DT105 * KP(POS,IE) + TMP2 * DELTA(POS,IE) BTMP(2) = TMP2 * DELTA(POS_TRICK(POS,1),IE) BTMP(3) = TMP2 * DELTA(POS_TRICK(POS,2),IE) B1(IP) = ( BTMP(1) * U(IP) + BTMP(2) * U(IP_J) + BTMP(3) * U(IP_K) ) + B1(IP) END DO END DO IF (ICOMP .GT. 0 .AND. SMETHOD .GT. 0) THEN DO IP = 1, MNP ASPAR1(I_DIAG(IP)) = ASPAR1(I_DIAG(IP)) + IMATDAA(IS,ID,IP) * DT1 * SI(IP) ! Add sourcee term to the diagonal B1(IP) = B1(IP) + IMATRAA(IS,ID,IP) * DT1 * SI(IP) ! Add source term to the right hand side END DO END IF IPAR(1) = 0 ! always 0 to start an iterative solver IPAR(2) = 1 ! right preconditioning IPAR(3) = 1 ! use convergence test scheme 1 IPAR(4) = 100*MNP ! the 'w' has 10,000 elements IPAR(5) = 10 ! use *GMRES(10) (e.g. FGMRES(10)) IPAR(6) = 1000 ! use at most 100 matvec's FPAR(1) = 1.0E-6 ! relative tolerance 1.0E-6 FPAR(2) = 1.0E-8 ! absolute tolerance 1.0E-10 FPAR(11) = 0.0 ! clearing the FLOPS counter IWK = IPAR(4) DROPTOL = 10E-2 PERMTOL = 0.01 LFIL = MNP - 6 MBLOC = 4 INIU(:) = 0. CALL ILU0 (MNP, ASPAR1, JA, IA, AU, JAU, JU, IWKSP, IERROR) CALL RUNRC(MNP, NNZ, B1, X, IPAR, FPAR, WKSP, INIU, ASPAR1, JA, IA, AU, JAU, JU, BCGSTAB) U(:) = MAX(X(:),ZERO) J = 0 DO IP = 1, MNP DO I = 1, CCON(IP) J = J + 1 IE = IE_CELL(J) B2(IP) = ONETHIRD * TRIA(IE) * U(IP) + B2(IP) END DO END DO IF (ICOMP .GT. 0 .AND. SMETHOD .GT. 0) THEN DO IP = 1, MNP ASPAR2(I_DIAG(IP)) = ASPAR2(I_DIAG(IP)) + IMATDAA(IS,ID,IP) * DT2 * SI(IP) ! Add sourcee term to the diagonal B2(IP) = B2(IP) + IMATRAA(IS,ID,IP) * DT2 * SI(IP) ! Add source term to the right hand side END DO END IF IPAR(1) = 0 ! always 0 to start an iterative solver IPAR(2) = 1 ! right preconditioning IPAR(3) = 1 ! use convergence test scheme 1 IPAR(4) = 100*MNP ! the 'w' has 10,000 elements IPAR(5) = 10 ! use *GMRES(10) (e.g. FGMRES(10)) IPAR(6) = 1000 ! use at most 100 matvec's FPAR(1) = 1.0E-6 ! relative tolerance 1.0E-6 FPAR(2) = 1.0E-8 ! absolute tolerance 1.0E-10 FPAR(11) = 0.0 ! clearing the FLOPS counter IWK = IPAR(4) DROPTOL = 10E-2 PERMTOL = 0.1 LFIL = MNP - 3 MBLOC = MNP INIU(:) = U CALL ILU0 (MNP, ASPAR2, JA, IA, AU, JAU, JU, IWKSP, IERROR) CALL RUNRC(MNP, NNZ, B2, X, IPAR, FPAR, WKSP, INIU, ASPAR2, JA, IA, AU, JAU, JU, BCGSTAB) AC2(IS,ID,:) = MAX(ZERO,X) * IOBPD(ID,:) END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE CADVXY(IS,ID,C) USE DATAPOOL IMPLICIT NONE INTEGER, INTENT(IN) :: IS, ID REAL(rkind), INTENT(OUT) :: C(2,MNP) INTEGER :: IP REAL(rkind) :: DIFRU, USOC, WVC ! ! Loop over the resident nodes only ... exchange is done in the calling routine ! IF (LADVTEST) THEN C(1,:) = YP C(2,:) = - XP ELSE DO IP = 1, MNP IF (LSECU .OR. LSTCU) THEN C(1,IP) = CG(IS,IP)*COSTH(ID)+CURTXY(IP,1) C(2,IP) = CG(IS,IP)*SINTH(ID)+CURTXY(IP,2) ELSE C(1,IP) = CG(IS,IP)*COSTH(ID) C(2,IP) = CG(IS,IP)*SINTH(ID) END IF IF (LSPHE) THEN C(1,IP) = C(1,IP)*INVSPHTRANS(IP,1) C(2,IP) = C(2,IP)*INVSPHTRANS(IP,2) END IF IF (LDIFR) THEN C(1,IP) = C(1,IP)*DIFRM(IP) C(2,IP) = C(2,IP)*DIFRM(IP) IF (LSECU .OR. LSTCU) THEN IF (IDIFFR .GT. 1) THEN WVC = SPSIG(IS)/WK(IS,IP) USOC = (COSTH(ID)*CURTXY(IP,1) + SINTH(ID)*CURTXY(IP,2))/WVC DIFRU = ONE + USOC * (ONE - DIFRM(IP)) ELSE DIFRU = DIFRM(IP) END IF C(1,IP) = C(1,IP) + DIFRU*CURTXY(IP,1) C(2,IP) = C(2,IP) + DIFRU*CURTXY(IP,2) END IF END IF ! LDIFR END DO END IF END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE CADVXY_VECTOR(CX,CY) USE DATAPOOL IMPLICIT NONE REAL(rkind), INTENT(OUT) :: CX(MSC,MDC,MNP), CY(MSC,MDC,MNP) INTEGER :: IP, IS, ID REAL(rkind) :: DIFRU, USOC, WVC ! ! Loop over the resident nodes only ... exchange is done in the calling routine ! DO IS = 1, MSC DO ID = 1, MDC DO IP = 1, MNP IF (LSECU .OR. LSTCU) THEN CX(IS,ID,IP) = CG(IS,IP)*COSTH(ID)+CURTXY(IP,1) CY(IS,ID,IP) = CG(IS,IP)*SINTH(ID)+CURTXY(IP,2) ELSE CX(IS,ID,IP) = CG(IS,IP)*COSTH(ID) CY(IS,ID,IP) = CG(IS,IP)*SINTH(ID) END IF IF (LSPHE) THEN CX(IS,ID,IP) = CX(IS,ID,IP)*INVSPHTRANS(IP,1) CY(IS,ID,IP) = CY(IS,ID,IP)*INVSPHTRANS(IP,2) END IF IF (LDIFR) THEN CX(IS,ID,IP) = CX(IS,ID,IP)*DIFRM(IP) CY(IS,ID,IP) = CY(IS,ID,IP)*DIFRM(IP) IF (LSECU .OR. LSTCU) THEN IF (IDIFFR .GT. 1) THEN WVC = SPSIG(IS)/WK(IS,IP) USOC = (COSTH(ID)*CURTXY(IP,1) + SINTH(ID)*CURTXY(IP,2))/WVC DIFRU = ONE + USOC * (ONE - DIFRM(IP)) ELSE DIFRU = DIFRM(IP) END IF CX(IS,ID,IP) = CX(IS,ID,IP) + DIFRU*CURTXY(IP,1) CY(IS,ID,IP) = CY(IS,ID,IP) + DIFRU*CURTXY(IP,2) END IF END IF END DO !IP END DO !ID END DO !IS END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE INIT_FLUCT_ARRAYS USE DATAPOOL IMPLICIT NONE ALLOCATE(CCON(MNP), SI(MNP), ITER_EXP(MSC,MDC), ITER_EXPD(MSC), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_fluctsplit, allocate error 1') CCON = 0 SI = ZERO ITER_EXP = 0 ITER_EXPD = 0 ALLOCATE( I_DIAG(MNP), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_fluctsplit, allocate error 2') I_DIAG = 0 IF (LCFL) THEN ALLOCATE (CFLCXY(3,MNP), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_fluctsplit, allocate error 3') CFLCXY(1,:) = ZERO CFLCXY(2,:) = ZERO CFLCXY(3,:) = LARGE END IF END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE DEALLOC_FLUCT_ARRAYS USE DATAPOOL IMPLICIT NONE DEALLOCATE( CCON, SI, ITER_EXP, ITER_EXPD) IF ((ICOMP .GE. 1) .OR. LZETA_SETUP) THEN DEALLOCATE(I_DIAG) END IF IF (LCFL) THEN DEALLOCATE(CFLCXY) END IF END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE NEXT_IDX(len, i, iNext) IMPLICIT NONE integer, intent(in) :: len, i integer, intent(out) :: iNext IF (i .lt. len) THEN iNext=i+1 ELSE iNext=1 END IF END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE PREV_IDX(len, i, iPrev) IMPLICIT NONE integer, intent(in) :: len, i integer, intent(out) :: iPrev IF (i .gt. 1) THEN iPrev=i-1 ELSE iPrev=len END IF END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE INIT_FLUCT USE DATAPOOL IMPLICIT NONE INTEGER :: I, J, K INTEGER :: IP, IE, POS, POS_J, POS_K, IP_I, IP_J, IP_K INTEGER :: I1, I2, I3 INTEGER :: CHILF(MNP), COUNT_MAX INTEGER :: ITMP(MNP) INTEGER :: POS_TRICK(3,2) INTEGER :: IADJ, NB_ADJ INTEGER :: POS1, POS2, J1, J2 INTEGER :: FPOS, EPOS, JP INTEGER, ALLOCATABLE :: REV_BOOK(:) REAL(rkind) :: TRIA03, TMP1, TMP2 INTEGER, ALLOCATABLE :: CELLVERTEX(:,:,:) INTEGER, ALLOCATABLE :: PTABLE(:,:) INTEGER :: SUM_CCON, SizeAlloc INTEGER nbMatch, IE2, ICON, INEXT, IE_ADJ INTEGER IP_NEXT, IP_ADJ_NEXT, IP_ADJ_PREV INTEGER POS_NEXT, POS_PREV LOGICAL CHECK_COMBIN_ORIENT POS_TRICK(1,1) = 2 POS_TRICK(1,2) = 3 POS_TRICK(2,1) = 3 POS_TRICK(2,2) = 1 POS_TRICK(3,1) = 1 POS_TRICK(3,2) = 2 WRITE(STAT%FHNDL,'("+TRACE......",A)') 'CALCULATE CONNECTED AREA SI ' ! The situation is as follows with respect to MNP, NP_RES and friends. ! For sparse matrix, it makes sense to compute ASPAR (the sparse matrix ! elements) only for IP=1,NP_RES ! For the element NP_RES+1,MNP (those are ghost nodes) we only have part ! of the containing elements and any computation there will be partial, ! and so incorrect. ! However, we do need to be able to support matrix elements from 1 to MNP ! The reason is that in order to do ILU0 preconditioning, we need to be ! able to access such elements and so we need memory allocated for them ! (but computed from other nodes) !OPEN(5555, FILE = 'fluctgeo.dat', STATUS = 'UNKNOWN') ! ! Calculate the max. number of connected elements count_max ! WRITE(STAT%FHNDL,'("+TRACE......",A)') 'MEDIAN DUAL AREA and CCON' SI(:) = 0.0d0 ! Median Dual Patch Area of each Node CCON(:) = 0 ! Number of connected Elements DO IE = 1 , MNE I1 = INE(1,IE) I2 = INE(2,IE) I3 = INE(3,IE) CCON(I1) = CCON(I1) + 1 CCON(I2) = CCON(I2) + 1 CCON(I3) = CCON(I3) + 1 TRIA03 = ONETHIRD * TRIA(IE) SI(I1) = SI(I1) + TRIA03 SI(I2) = SI(I2) + TRIA03 SI(I3) = SI(I3) + TRIA03 !WRITE(STAT%FHNDL,*) IE, TRIA(IE) ENDDO !DO IP = 1, MNP !WRITE(STAT%FHNDL,*) IP, SI(IP) !ENDDO #ifdef MPI_PARALL_GRID CALL EXCHANGE_P2D(SI) #endif ! We don't need MAXMNECON from SCHISM/pdlib if we compute CCON itself ! #ifdef MPI_PARALL_GRID ! MAXMNECON = MNEI ! #else ! MAXMNECON = MAXVAL(CCON) ! #endif MAXMNECON = MAXVAL(CCON) ! check agains SCHISM to make sure that there is no problem #ifdef MPI_PARALL_GRID # ifndef PDLIB IF (MAXMNECON /= MNEI) THEN write(DBG%FHNDL,*) "WARNING", __FILE__ , "Line", __LINE__ write(DBG%FHNDL,*) "MAXMNECON from SCHISM does not match self calc value. This could be problems", MAXMNECON, MNEI END IF # endif #endif ! WRITE(STAT%FHNDL,'("+TRACE......",A)') 'CALCULATE FLUCTUATION POINTER' ALLOCATE(CELLVERTEX(MNP,MAXMNECON,2), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_fluctsplit, allocate error 4') ! CELLVERTEX(:,:,:) = 0 ! Stores for each node the Elementnumbers of the connected Elements ! and the Position of the position of the Node in the Element Index CHILF = 0 DO IE = 1, MNE DO J=1,3 I = INE(J,IE) CHILF(I) = CHILF(I)+1 CELLVERTEX(I,CHILF(I),1) = IE CELLVERTEX(I,CHILF(I),2) = J END DO ENDDO ! ! Emulates loop structure and counts max. entries in the different pointers that have to be designed ! J = 0 DO IP = 1, MNP DO I = 1, CCON(IP) J = J + 1 END DO END DO COUNT_MAX = J ! Max. Number of entries in the pointers used in the calculations IF (COUNT_MAX.ne.3*MNE) THEN WRITE(DBG%FHNDL,*) 'COUNT_MAX=', COUNT_MAX WRITE(DBG%FHNDL,*) 'MNE=', MNE CALL WWM_ABORT('Do Not Sleep Before solving the problem') ENDIF ALLOCATE(IE_CELL(COUNT_MAX), POS_CELL(COUNT_MAX), IE_CELL2(MNP,MAXMNECON), POS_CELL2(MNP,MAXMNECON), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_fluctsplit, allocate error 5') ! Just a remapping from CELLVERTEX ... Element number in the ! order of the occurence in the loop during runtime IE_CELL = 0 ! Just a remapping from CELLVERTEX ... Position of the node ! in the Element index -"- POS_CELL = 0 ! J = 0 DO IP = 1, MNP DO I = 1, CCON(IP) J = J + 1 IE_CELL(J) = CELLVERTEX(IP,I,1) POS_CELL(J) = CELLVERTEX(IP,I,2) IE_CELL2(IP,I) = CELLVERTEX(IP,I,1) POS_CELL2(IP,I) = CELLVERTEX(IP,I,2) END DO END DO DEALLOCATE(CELLVERTEX) CHECK_COMBIN_ORIENT=.TRUE. IF (CHECK_COMBIN_ORIENT) THEN DO IE=1,MNE DO I=1,3 INEXT=POS_TRICK(I,1) IP=INE(I, IE) IP_NEXT=INE(I, IE) nbMatch=0 IE_ADJ=-1 DO ICON=1,CCON(IP) IE2=IE_CELL2(IP,ICON) IF (IE .ne. IE2) THEN POS=POS_CELL2(IP, ICON) POS_NEXT=POS_TRICK(POS,1) IP_ADJ_NEXT=INE(POS_NEXT,IE2) IF (IP_ADJ_NEXT .eq. IP_NEXT) THEN WRITE(DBG%FHNDL,*) 'Combinatorial orientability problem' WRITE(DBG%FHNDL,*) 'IE=', IE, ' IE2=', IE2 WRITE(DBG%FHNDL,*) 'IP=', IP, ' IP_NEXT=', IP_NEXT FLUSH(DBG%FHNDL) CALL WWM_ABORT('Please correct the grid combinatorially 1') END IF POS_PREV=POS_TRICK(POS,2) IP_ADJ_PREV=INE(POS_PREV,IE2) IF (IP_ADJ_PREV .eq. IP_NEXT) THEN nbMatch=nbMatch+1 IE_ADJ=IE2 END IF END IF END DO IF (nbMatch .gt. 1) THEN WRITE(DBG%FHNDL,*) 'nbMatch is too large.' WRITE(DBG%FHNDL,*) 'Should be 0 for boundary edge' WRITE(DBG%FHNDL,*) 'Should be 1 for interior edges' WRITE(DBG%FHNDL,*) 'nbMatch=', nbMatch FLUSH(DBG%FHNDL) CALL WWM_ABORT('Please correct the grid combinatorially 2') END IF END DO END DO END IF IF (ICOMP .GT. 0 .OR. LEXPIMP .OR. LZETA_SETUP) THEN ALLOCATE(PTABLE(COUNT_MAX,7), JA_IE(3,3,MNE), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_fluctsplit, allocate error 6') J = 0 PTABLE(:,:) = 0. ! Table storing some other values needed to design the sparse matrix pointers. DO IP = 1, MNP DO I = 1, CCON(IP) J = J + 1 IE = IE_CELL(J) POS = POS_CELL(J) I1 = INE(1,IE) I2 = INE(2,IE) I3 = INE(3,IE) IF (POS == 1) THEN POS_J = 2 POS_K = 3 ELSE IF (POS == 2) THEN POS_J = 3 POS_K = 1 ELSE POS_J = 1 POS_K = 2 END IF IP_I = IP IP_J = INE(POS_J,IE) IP_K = INE(POS_K,IE) PTABLE(J,1) = IP_I ! Node numbers of the connected elements PTABLE(J,2) = IP_J PTABLE(J,3) = IP_K PTABLE(J,4) = POS ! Position of the nodes in the element index PTABLE(J,5) = POS_J PTABLE(J,6) = POS_K PTABLE(J,7) = IE ! Element numbers same as IE_CELL END DO END DO WRITE(STAT%FHNDL,'("+TRACE......",A)') 'SET UP SPARSE MATRIX POINTER ... COUNT NONZERO ENTRY' ! ! Count number of nonzero entries in the matrix ... ! Basically, each connected element may have two off-diagonal ! contribution and one diagonal related to the connected vertex itself ... ! J = 0 NNZ = 0 DO IP = 1, MNP !AR: bug below ... ITMP(:) = 0 DO I = 1, CCON(IP) J = J + 1 IP_J = PTABLE(J,2) IP_K = PTABLE(J,3) ITMP(IP) = 1 ITMP(IP_J) = 1 ITMP(IP_K) = 1 END DO NNZ = NNZ + SUM(ITMP) END DO WRITE(STAT%FHNDL,'("+TRACE......",A)') 'SET UP SPARSE MATRIX POINTER ... SETUP POINTER' ! ! Allocate sparse matrix pointers using the Compressed Sparse Row Format CSR ... this is now done only of MNP nodes ! The next step is to do it for the whole Matrix MNP * MSC * MDC ! see ...:x ! ! JA Pointer according to the convention in my thesis see p. 123 ! IA Pointer according to the convention in my thesis see p. 123 ALLOCATE (JA(NNZ), IA(MNP+1), POSI(3,COUNT_MAX), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_fluctsplit, allocate error 6') JA = 0 IA = 0 POSI = 0 ! Points to the position of the matrix entry in the mass matrix ! according to the CSR matrix format see p. 124 J = 0 K = 0 IA(1) = 1 MAX_DEG=0 DO IP = 1, MNP ! Run through all rows ITMP=0 DO I = 1, CCON(IP) ! Check how many entries there are ... J = J + 1 IP_J = PTABLE(J,2) IP_K = PTABLE(J,3) ITMP(IP) = 1 ITMP(IP_J) = 1 ITMP(IP_K) = 1 END DO IADJ=0 DO I = 1, MNP ! Run through all columns IF (ITMP(I) .GT. 0) THEN K = K + 1 IF (I .ne. IP) THEN IADJ=IADJ + 1 END IF JA(K) = I END IF END DO IF (IADJ .gt. MAX_DEG) THEN MAX_DEG=IADJ END IF IA(IP + 1) = K + 1 END DO J = 0 DO IP = 1, MNP DO I = 1, CCON(IP) J = J + 1 IP_J = PTABLE(J,2) IP_K = PTABLE(J,3) DO K = IA(IP), IA(IP+1) - 1 IF (IP == JA(K)) POSI(1,J) = K IF (IP == JA(K)) I_DIAG(IP) = K IF (IP_J == JA(K)) POSI(2,J) = K IF (IP_K == JA(K)) POSI(3,J) = K END DO END DO END DO J=0 DO IP=1,MNP DO I = 1, CCON(IP) J=J+1 IE = IE_CELL(J) POS = POS_CELL(J) I1 = POSI(1,J) I2 = POSI(2,J) I3 = POSI(3,J) JA_IE(POS,1,IE) = I1 JA_IE(POS,2,IE) = I2 JA_IE(POS,3,IE) = I3 END DO END DO ! ! Arrays for Jacobi-Gauss-Seidel solver ! IF (AMETHOD .eq. 7) THEN IF (ASPAR_LOCAL_LEVEL .eq. 2) THEN ALLOCATE(LIST_ADJ_VERT(MAX_DEG,MNP), VERT_DEG(MNP), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_fluctsplit, allocate error 6a') J = 0 K = 0 DO IP = 1, MNP ITMP=0 DO I = 1, CCON(IP) ! Check how many entries there are ... J = J + 1 IP_J = PTABLE(J,2) IP_K = PTABLE(J,3) ITMP(IP) = 1 ITMP(IP_J) = 1 ITMP(IP_K) = 1 END DO IADJ=0 DO I = 1, MNP IF (ITMP(I) .GT. 0) THEN K = K + 1 IF (I .ne. IP) THEN IADJ=IADJ + 1 LIST_ADJ_VERT(IADJ,IP)=I END IF JA(K) = I END IF END DO VERT_DEG(IP)=IADJ IA(IP + 1) = K + 1 END DO ALLOCATE(REV_BOOK(NNZ), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_fluctsplit, allocate error 6b') REV_BOOK=0 DO IP=1,MNP DO J=IA(IP),IA(IP+1)-1 JP=JA(J) IF (IP .ne. JP) THEN FPOS=-1 DO EPOS=1,MAX_DEG IF (LIST_ADJ_VERT(EPOS,IP) .eq. JP) THEN FPOS=EPOS END IF END DO IF (FPOS .eq. -1) THEN CALL WWM_ABORT('INDEXING FAILURE') END IF ELSE FPOS=-400 END IF REV_BOOK(J)=FPOS END DO END DO ALLOCATE(POS_IP_ADJ(2,3,MNE), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_fluctsplit, allocate error 6c') DO IE=1,MNE DO J=1,3 J1=JA_IE(J,2,IE) J2=JA_IE(J,3,IE) POS1=REV_BOOK(J1) POS2=REV_BOOK(J2) POS_IP_ADJ(1,J,IE)=POS1 POS_IP_ADJ(2,J,IE)=POS2 END DO END DO DEALLOCATE(REV_BOOK) END IF IF (ASPAR_LOCAL_LEVEL .le. 1) THEN ALLOCATE (ASPAR_JAC(MSC,MDC,NNZ), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_initio, allocate error 9') ASPAR_JAC = zero TMP1 = MyREAL(MSC)*MyREAL(MDC)*MyREAL(NNZ*8) TMP2 = 1024**2 WRITE(STAT%FHNDL,'("+TRACE......",A,F15.4,A)') 'MAX MEMORY SIZE OF ASPAR_JAC =', TMP1/TMP2, 'MB' TMP1 = TMP1 + MyREAL(MSC) * MyREAL(MDC) * MyREAL(MNP) * 4 * 8 WRITE(STAT%FHNDL,'("+TRACE......",A,F15.4,A)') 'TOTAL MEMORY SIZE =', TMP1/TMP2, 'MB' END IF IF ((ASPAR_LOCAL_LEVEL .ge. 5).and.(ASPAR_LOCAL_LEVEL .le. 7)) THEN SUM_CCON = 0 DO IP = 1, NP_RES SUM_CCON = SUM_CCON +CCON(IP) END DO ALLOCATE(K_CRFS_MSC(12, MSC,SUM_CCON), K_CRFS_U(6,SUM_CCON), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_initio, allocate error 9b') SizeAlloc=12*MSC*SUM_CCON WRITE(STAT%FHNDL,*) 'LEVEL 5, nb =', SizeAlloc SizeAlloc=MDC*MSC*NNZ WRITE(STAT%FHNDL,*) 'ASPAR_JAC, nb =', SizeAlloc SizeAlloc=MDC*MSC*MNP WRITE(STAT%FHNDL,*) 'U_JAC, nb =', SizeAlloc END IF ! IF ((.NOT. LNONL) .AND. SOURCE_IMPL .AND. (ASPAR_LOCAL_LEVEL.le.1)) THEN ALLOCATE (B_JAC(MSC,MDC,MNP), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_initio, allocate error 9') B_JAC = zero END IF END IF DEALLOCATE(PTABLE) ENDIF END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE DEALLOC_FLUCT USE DATAPOOL implicit none DEALLOCATE (IE_CELL, POS_CELL, IE_CELL2, POS_CELL2) IF (ICOMP .GT. 0 .OR. LEXPIMP) THEN DEALLOCATE (JA, IA, POSI) END IF END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE EXPLICIT_N_SCHEME_VECTOR USE DATAPOOL IMPLICIT NONE ! ! local integer ! INTEGER :: IP, IE, IT, IS, ID INTEGER :: I1, I2, I3 INTEGER :: NI(3), I, IPOS ! ! local double ! REAL(rkind) :: UTILDE REAL(rkind) :: DTMAX_GLOBAL_EXP, DTMAX_EXP #ifdef MPI_PARALL_GRID REAL(rkind) :: WILD(MNP), WILD2D(MDC,MNP) #endif REAL(rkind) :: REST, CFLXY REAL(rkind) :: LAMBDA(2,MSC,MDC), DT4AI REAL(rkind) :: FL11(MSC,MDC),FL12(MSC,MDC),FL21(MSC,MDC),FL22(MSC,MDC),FL31(MSC,MDC),FL32(MSC,MDC) REAL(rkind) :: KTMP(3,MSC,MDC) REAL(rkind) :: U3(3) REAL(rkind) :: KKSUM(MNP,MSC,MDC), ST(MNP,MSC,MDC), N(MNE,MSC,MDC) REAL(rkind) :: CX(MSC,MDC,MNP), CY(MSC,MDC,MNP) REAL(rkind) :: FLALL(3,MNE,MSC,MDC) REAL(rkind) :: KELEM(3,MNE,MSC,MDC) REAL(rkind) :: FL111(MSC,MDC), FL112(MSC,MDC), FL211(MSC,MDC), FL212(MSC,MDC), FL311(MSC,MDC), FL312(MSC,MDC) REAL(rkind) :: UTILDE3(MNE) REAL(rkind) :: USOC, WVC, DIFRU REAL(rkind) :: TIME1, TIME2 ! ! local parameter ! REAL(rkind) :: TMP(MSC,MDC) ! ! Calculate phase speeds for the certain spectral component ... ! FLALL = ZERO KELEM = ZERO KKSUM = ZERO ST = ZERO N = ZERO DO IS = 1, MSC DO ID = 1, MDC DO IP = 1, MNP IF (LSECU .OR. LSTCU) THEN CX(IS,ID,IP) = CG(IS,IP)*COSTH(ID)+CURTXY(IP,1) CY(IS,ID,IP) = CG(IS,IP)*SINTH(ID)+CURTXY(IP,2) ELSE CX(IS,ID,IP) = CG(IS,IP)*COSTH(ID) CY(IS,ID,IP) = CG(IS,IP)*SINTH(ID) END IF IF (LSPHE) THEN CX(IS,ID,IP) = CX(IS,ID,IP)*INVSPHTRANS(IP,1) CY(IS,ID,IP) = CY(IS,ID,IP)*INVSPHTRANS(IP,2) END IF IF (LDIFR) THEN CX(IS,ID,IP) = CX(IS,ID,IP)*DIFRM(IP) CY(IS,ID,IP) = CY(IS,ID,IP)*DIFRM(IP) IF (LSECU .OR. LSTCU) THEN IF (IDIFFR .GT. 1) THEN WVC = SPSIG(IS)/WK(IS,IP) USOC = (COSTH(ID)*CURTXY(IP,1) + SINTH(ID)*CURTXY(IP,2))/WVC DIFRU = ONE + USOC * (ONE - DIFRM(IP)) ELSE DIFRU = DIFRM(IP) END IF CX(IS,ID,IP) = CX(IS,ID,IP) + DIFRU*CURTXY(IP,1) CY(IS,ID,IP) = CY(IS,ID,IP) + DIFRU*CURTXY(IP,2) END IF END IF END DO END DO END DO ! ! Calculate K-Values and contour based quantities ... ! !!$OMP DO PRIVATE(IE,I1,I2,I3,LAMBDA,KTMP,TMP,FL11,FL12,FL21,FL22,FL31,FL32,FL111,FL112,FL211,FL212,FL311,FL312) DO IE = 1, MNE I1 = INE(1,IE) I2 = INE(2,IE) I3 = INE(3,IE) LAMBDA(1,:,:) = ONESIXTH *(CX(:,:,I1)+CX(:,:,I2)+CX(:,:,I3)) LAMBDA(2,:,:) = ONESIXTH *(CY(:,:,I1)+CY(:,:,I2)+CY(:,:,I3)) KELEM(1,IE,:,:) = LAMBDA(1,:,:) * IEN(1,IE) + LAMBDA(2,:,:) * IEN(2,IE) KELEM(2,IE,:,:) = LAMBDA(1,:,:) * IEN(3,IE) + LAMBDA(2,:,:) * IEN(4,IE) KELEM(3,IE,:,:) = LAMBDA(1,:,:) * IEN(5,IE) + LAMBDA(2,:,:) * IEN(6,IE) KTMP(1,:,:) = KELEM(1,IE,:,:) KTMP(2,:,:) = KELEM(2,IE,:,:) KTMP(3,:,:) = KELEM(3,IE,:,:) TMP(:,:) = SUM(MIN(ZERO,KTMP(:,:,:)),DIM=1) N(IE,:,:) = -ONE/MIN(-THR,TMP(:,:)) KELEM(1,IE,:,:) = MAX(ZERO,KTMP(1,:,:)) KELEM(2,IE,:,:) = MAX(ZERO,KTMP(2,:,:)) KELEM(3,IE,:,:) = MAX(ZERO,KTMP(3,:,:)) ! WRITE(DBG%FHNDL,'(3I10,3F15.4)') IS, ID, IE, KELEM(:,IE) FL11 = CX(:,:,I2) * IEN(1,IE) + CY(:,:,I2) * IEN(2,IE) FL12 = CX(:,:,I3) * IEN(1,IE) + CY(:,:,I3) * IEN(2,IE) FL21 = CX(:,:,I3) * IEN(3,IE) + CY(:,:,I3) * IEN(4,IE) FL22 = CX(:,:,I1) * IEN(3,IE) + CY(:,:,I1) * IEN(4,IE) FL31 = CX(:,:,I1) * IEN(5,IE) + CY(:,:,I1) * IEN(6,IE) FL32 = CX(:,:,I2) * IEN(5,IE) + CY(:,:,I2) * IEN(6,IE) FL111 = TWO*FL11+FL12 FL112 = TWO*FL12+FL11 FL211 = TWO*FL21+FL22 FL212 = TWO*FL22+FL21 FL311 = TWO*FL31+FL32 FL312 = TWO*FL32+FL31 FLALL(1,IE,:,:) = (FL311 + FL212) * ONESIXTH + KELEM(1,IE,:,:) FLALL(2,IE,:,:) = (FL111 + FL312) * ONESIXTH + KELEM(2,IE,:,:) FLALL(3,IE,:,:) = (FL211 + FL112) * ONESIXTH + KELEM(3,IE,:,:) END DO IF (LCALC) THEN KKSUM = ZERO DO IE = 1, MNE NI = INE(:,IE) KKSUM(NI(1),:,:) = KKSUM(NI(1),:,:) + KELEM(1,IE,:,:) KKSUM(NI(2),:,:) = KKSUM(NI(2),:,:) + KELEM(2,IE,:,:) KKSUM(NI(3),:,:) = KKSUM(NI(3),:,:) + KELEM(3,IE,:,:) END DO IF (IVECTOR == 1) THEN DO ID = 1, MDC DO IS = 1, MSC DTMAX_GLOBAL_EXP = VERYLARGE #ifdef MPI_PARALL_GRID DO IP = 1, MNP DTMAX_EXP = SI(IP)/MAX(VERYSMALL,KKSUM(IP,IS,ID)) DTMAX_GLOBAL_EXP = MIN ( DTMAX_GLOBAL_EXP, DTMAX_EXP) END DO DTMAX_EXP=DTMAX_GLOBAL_EXP call mpi_allreduce(DTMAX_EXP,DTMAX_GLOBAL_EXP,1,rtype,MPI_MIN,comm,ierr) #else DO IP = 1, MNP DTMAX_EXP = SI(IP)/MAX(VERYSMALL,KKSUM(IP,IS,ID)) DTMAX_GLOBAL_EXP = MIN ( DTMAX_GLOBAL_EXP, DTMAX_EXP) END DO #endif CFLXY = DT4A/DTMAX_GLOBAL_EXP REST = ABS(MOD(CFLXY,ONE)) IF (REST .LT. THR) THEN ITER_EXP(IS,ID) = ABS(NINT(CFLXY)) ELSE IF (REST .GT. THR .AND. REST .LT. ONEHALF) THEN ITER_EXP(IS,ID) = ABS(NINT(CFLXY)) + 1 ELSE ITER_EXP(IS,ID) = ABS(NINT(CFLXY)) END IF END DO END DO ELSE IF (IVECTOR .EQ. 2) THEN DTMAX_GLOBAL_EXP = VERYLARGE #ifdef MPI_PARALL_GRID DO IP = 1, NP_RES DTMAX_EXP = SI(IP)/MAX(THR,MAXVAL(KKSUM(IP,:,:))) DTMAX_GLOBAL_EXP = MIN ( DTMAX_GLOBAL_EXP, DTMAX_EXP) END DO DTMAX_EXP=DTMAX_GLOBAL_EXP call mpi_allreduce(DTMAX_EXP,DTMAX_GLOBAL_EXP,1,rtype,MPI_MIN,comm,ierr) #else DO IP = 1, MNP DTMAX_EXP = SI(IP)/MAX(THR,MAXVAL(KKSUM(IP,:,:))) DTMAX_GLOBAL_EXP = MIN ( DTMAX_GLOBAL_EXP, DTMAX_EXP) END DO #endif CFLXY = DT4A/DTMAX_GLOBAL_EXP REST = ABS(MOD(CFLXY,ONE)) IF (REST .LT. THR) THEN ITER_MAX = ABS(NINT(CFLXY)) ELSE IF (REST .GT. THR .AND. REST .LT. ONEHALF) THEN ITER_MAX = ABS(NINT(CFLXY)) + 1 ELSE ITER_MAX = ABS(NINT(CFLXY)) END IF ELSE IF (IVECTOR .EQ. 3) THEN DO IS = 1, MSC DTMAX_GLOBAL_EXP = VERYLARGE #ifdef MPI_PARALL_GRID DO IP = 1, NP_RES DTMAX_EXP = SI(IP)/MAX(VERYSMALL,MAXVAL(KKSUM(IP,IS,:))) DTMAX_GLOBAL_EXP = MIN ( DTMAX_GLOBAL_EXP, DTMAX_EXP) END DO DTMAX_EXP=DTMAX_GLOBAL_EXP call mpi_allreduce(DTMAX_EXP,DTMAX_GLOBAL_EXP,1,rtype,MPI_MIN,comm,ierr) #else DO IP = 1, MNP DTMAX_EXP = SI(IP)/MAX(VERYSMALL,MAXVAL(KKSUM(IP,IS,:))) DTMAX_GLOBAL_EXP = MIN ( DTMAX_GLOBAL_EXP, DTMAX_EXP) END DO #endif CFLXY = DT4A/DTMAX_GLOBAL_EXP REST = ABS(MOD(CFLXY,1._rkind)) IF (REST .LT. THR) THEN ITER_EXPD(IS) = ABS(NINT(CFLXY)) ELSE IF (REST .GT. THR .AND. REST .LT. ONEHALF) THEN ITER_EXPD(IS) = ABS(NINT(CFLXY)) + 1 ELSE ITER_EXPD(IS) = ABS(NINT(CFLXY)) END IF END DO ELSE IF (IVECTOR .EQ. 4) THEN DO IS = 1, MSC DTMAX_GLOBAL_EXP = VERYLARGE #ifdef MPI_PARALL_GRID DO IP = 1, NP_RES DTMAX_EXP = SI(IP)/MAX(VERYSMALL,MAXVAL(KKSUM(IP,IS,:))) DTMAX_GLOBAL_EXP = MIN ( DTMAX_GLOBAL_EXP, DTMAX_EXP) END DO DTMAX_EXP=DTMAX_GLOBAL_EXP call mpi_allreduce(DTMAX_EXP,DTMAX_GLOBAL_EXP,1,rtype,MPI_MIN,comm,ierr) #else DO IP = 1, MNP DTMAX_EXP = SI(IP)/MAX(VERYSMALL,MAXVAL(KKSUM(IP,IS,:))) DTMAX_GLOBAL_EXP = MIN ( DTMAX_GLOBAL_EXP, DTMAX_EXP) END DO #endif CFLXY = DT4A/DTMAX_GLOBAL_EXP REST = ABS(MOD(CFLXY,1._rkind)) IF (REST .LT. THR) THEN ITER_EXPD(IS) = ABS(NINT(CFLXY)) ELSE IF (REST .GT. THR .AND. REST .LT. ONEHALF) THEN ITER_EXPD(IS) = ABS(NINT(CFLXY)) + 1 ELSE ITER_EXPD(IS) = ABS(NINT(CFLXY)) END IF END DO ELSE IF (IVECTOR .EQ. 5) THEN DTMAX_GLOBAL_EXP = VERYLARGE #ifdef MPI_PARALL_GRID DO IP = 1, NP_RES DTMAX_EXP = SI(IP)/MAX(THR,MAXVAL(KKSUM(IP,:,:))) DTMAX_GLOBAL_EXP = MIN ( DTMAX_GLOBAL_EXP, DTMAX_EXP) END DO DTMAX_EXP=DTMAX_GLOBAL_EXP call mpi_allreduce(DTMAX_EXP,DTMAX_GLOBAL_EXP,1,rtype,MPI_MIN,comm,ierr) #else DO IP = 1, MNP DTMAX_EXP = SI(IP)/MAX(THR,MAXVAL(KKSUM(IP,:,:))) DTMAX_GLOBAL_EXP = MIN ( DTMAX_GLOBAL_EXP, DTMAX_EXP) END DO #endif CFLXY = DT4A/DTMAX_GLOBAL_EXP REST = ABS(MOD(CFLXY,ONE)) IF (REST .LT. THR) THEN ITER_MAX = ABS(NINT(CFLXY)) ELSE IF (REST .GT. THR .AND. REST .LT. ONEHALF) THEN ITER_MAX = ABS(NINT(CFLXY)) + 1 ELSE ITER_MAX = ABS(NINT(CFLXY)) END IF END IF !IVECTOR WRITE(STAT%FHNDL,*) 'MAX. ITERATIONS USED IN ADV. SCHEME', ITER_MAX, MAXVAL(ITER_EXP) FLUSH(STAT%FHNDL) END IF !LCALC #ifdef TIMINGS CALL WAV_MY_WTIME(TIME1) #endif IF (IVECTOR == 1) THEN DO ID = 1, MDC DO IS = 1, MSC DT4AI = DT4A/ITER_EXP(IS,ID) DO IT = 1, ITER_EXP(IS,ID) ST(:,IS,ID) = ZERO ! Init. ... only used over the residual nodes see IP loop DO IE = 1, MNE NI = INE(:,IE) U3(:) = AC2(IS,ID,NI) UTILDE = N(IE,IS,ID) * ( FLALL(1,IE,IS,ID) * U3(1) + FLALL(2,IE,IS,ID) * U3(2) + FLALL(3,IE,IS,ID) * U3(3) ) ST(NI(1),IS,ID) = ST(NI(1),IS,ID) + KELEM(1,IE,IS,ID) * (U3(1) - UTILDE) ST(NI(2),IS,ID) = ST(NI(2),IS,ID) + KELEM(2,IE,IS,ID) * (U3(2) - UTILDE) ST(NI(3),IS,ID) = ST(NI(3),IS,ID) + KELEM(3,IE,IS,ID) * (U3(3) - UTILDE) END DO !IE AC2(IS,ID,:) = MAX(ZERO,AC2(IS,ID,:)-DT4AI/SI*ST(:,IS,ID)*IOBWB)*IOBPD(ID,:) #ifdef MPI_PARALL_GRID WILD = AC2(IS,ID,:) CALL EXCHANGE_P2D(WILD) AC2(IS,ID,:) = WILD #endif END DO ! IT----> End Iteration END DO !IS END DO !ID ELSE IF (IVECTOR == 2) THEN DT4AI = DT4A/ITER_MAX DO IT = 1, ITER_MAX DO ID = 1, MDC DO IS = 1, MSC ST(:,IS,ID) = ZERO ! Init. ... only used over the residual nodes see IP loop DO IE = 1, MNE NI = INE(:,IE) U3(:) = AC2(IS,ID,NI) UTILDE = N(IE,IS,ID) * ( FLALL(1,IE,IS,ID) * U3(1) + FLALL(2,IE,IS,ID) * U3(2) + FLALL(3,IE,IS,ID) * U3(3) ) ST(NI(1),IS,ID) = ST(NI(1),IS,ID) + KELEM(1,IE,IS,ID) * (U3(1) - UTILDE) ST(NI(2),IS,ID) = ST(NI(2),IS,ID) + KELEM(2,IE,IS,ID) * (U3(2) - UTILDE) ST(NI(3),IS,ID) = ST(NI(3),IS,ID) + KELEM(3,IE,IS,ID) * (U3(3) - UTILDE) END DO !IE AC2(IS,ID,:) = MAX(ZERO,AC2(IS,ID,:)-DT4AI/SI*ST(:,IS,ID)*IOBWB)*IOBPD(ID,:) END DO !IS END DO !ID #ifdef MPI_PARALL_GRID CALL EXCHANGE_P4D_WWM(AC2) #endif END DO !IT ELSE IF (IVECTOR == 3) THEN DO IS = 1, MSC ITER_MAX = ITER_EXPD(IS) DT4AI = DT4A/ITER_MAX DO IT = 1, ITER_MAX DO ID = 1, MDC ST(:,IS,ID) = ZERO ! Init. ... only used over the residual nodes see IP loop DO IE = 1, MNE NI = INE(:,IE) U3(:) = AC2(IS,ID,NI) UTILDE = N(IE,IS,ID) * ( FLALL(1,IE,IS,ID) * U3(1) + FLALL(2,IE,IS,ID) * U3(2) + FLALL(3,IE,IS,ID) * U3(3) ) ST(NI(1),IS,ID) = ST(NI(1),IS,ID) + KELEM(1,IE,IS,ID) * (U3(1) - UTILDE) ST(NI(2),IS,ID) = ST(NI(2),IS,ID) + KELEM(2,IE,IS,ID) * (U3(2) - UTILDE) ST(NI(3),IS,ID) = ST(NI(3),IS,ID) + KELEM(3,IE,IS,ID) * (U3(3) - UTILDE) END DO !IE AC2(IS,ID,:) = MAX(ZERO,AC2(IS,ID,:)-DT4AI/SI*ST(:,IS,ID)*IOBWB)*IOBPD(ID,:) END DO! ID #ifdef MPI_PARALL_GRID WILD2D = AC2(IS,:,:) CALL EXCHANGE_P3D_WWM(WILD2D) AC2(IS,:,:) = WILD2D #endif END DO !IT END DO !IS ELSE IF (IVECTOR == 4) THEN DO IS = 1, MSC ITER_MAX = ITER_EXPD(IS) DT4AI = DT4A/ITER_MAX DO IT = 1, ITER_MAX DO ID = 1, MDC DO IE = 1, MNE NI = INE(:,IE) U3(:) = AC2(IS,ID,NI) UTILDE3(IE) = N(IE,IS,ID) * ( FLALL(1,IE,IS,ID) * U3(1) + FLALL(2,IE,IS,ID) * U3(2) + FLALL(3,IE,IS,ID) * U3(3) ) END DO ST(:,IS,ID) = ZERO DO IP = 1, MNP DO I = 1, CCON(IP) IE = IE_CELL2(IP,I) IPOS = POS_CELL2(IP,I) ST(IP,IS,ID) = ST(IP,IS,ID) + KELEM(IPOS,IE,IS,ID) * (AC2(IS,ID,IP) - UTILDE3(IE)) END DO AC2(IS,ID,IP) = MAX(ZERO,AC2(IS,ID,IP)-DT4AI/SI(IP)*ST(IP,IS,ID)*IOBWB(IP))*IOBPD(ID,IP) END DO END DO #ifdef MPI_PARALL_GRID WILD2D = AC2(IS,:,:) CALL EXCHANGE_P3D_WWM(WILD2D) AC2(IS,:,:) = WILD2D #endif END DO !IT END DO !IS ELSE IF (IVECTOR == 5) THEN DT4AI = DT4A/ITER_MAX DO IT = 1, ITER_MAX DO IS = 1, MSC DO ID = 1, MDC DO IE = 1, MNE NI = INE(:,IE) U3(:) = AC2(IS,ID,NI) UTILDE3(IE) = N(IE,IS,ID) * ( FLALL(1,IE,IS,ID) * U3(1) + FLALL(2,IE,IS,ID) * U3(2) + FLALL(3,IE,IS,ID) * U3(3) ) END DO !IE ST(:,IS,ID) = ZERO DO IP = 1, MNP DO I = 1, CCON(IP) IE = IE_CELL2(IP,I) IPOS = POS_CELL2(IP,I) ST(IP,IS,ID) = ST(IP,IS,ID) + KELEM(IPOS,IE,IS,ID) * (AC2(IS,ID,IP) - UTILDE3(IE)) END DO AC2(IS,ID,IP) = MAX(ZERO,AC2(IS,ID,IP)-DT4AI/SI(IP)*ST(IP,IS,ID)*IOBWB(IP))*IOBPD(ID,IP) END DO !IP END DO !ID END DO !IS #ifdef MPI_PARALL_GRID CALL EXCHANGE_P4D_WWM(AC2) #endif END DO !IT END IF END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE EXPLICIT_N_SCHEME_VECTOR_HPCF USE DATAPOOL IMPLICIT NONE ! ! local integer ! INTEGER :: IP, IE, IT, IS, ID INTEGER :: I1, I2, I3 INTEGER :: NI(3), I, IPOS ! ! local double ! REAL(rkind) :: DTMAX_GLOBAL_EXP, DTMAX_EXP #ifdef MPI_PARALL_GRID REAL(rkind) :: DTMAX_GLOBAL_EXP_LOC #endif REAL(rkind) :: REST, CFLXY REAL(rkind) :: LAMBDA(2), DT4AI REAL(rkind) :: FL11,FL12,FL21,FL22,FL31,FL32 REAL(rkind) :: U3(3), UIP(MNP) REAL(rkind) :: KKSUM, ST, N REAL(rkind) :: CX(3), CY(3) REAL(rkind) :: FLALL(3) REAL(rkind) :: KELEM(3) REAL(rkind) :: FL111, FL112, FL211, FL212, FL311, FL312 REAL(rkind) :: UTILDE3 REAL(rkind) :: KSUM(MNP), KMAX(MNP) REAL(rkind) :: TIME1, TIME2 ! ! local parameter ! ! Calculate phase speeds for the certain spectral component ... ! FLALL = ZERO KELEM = ZERO KKSUM = ZERO ST = ZERO N = ZERO IF (LCALC) THEN DO ID = 1, MDC DO IS = 1, MSC KMAX = ZERO KSUM = ZERO DO IP = 1, MNP DO I = 1, CCON(IP) IE = IE_CELL2(IP,I) IPOS = POS_CELL2(IP,I) ! get node indices from the element table ... NI = INE(:,IE) ! estimate speed in WAE CX = (CG(NI,IS)*COSTH(ID)+CURTXY(NI,1))* INVSPHTRANS(IP,1) CY =( CG(NI,IS)*SINTH(ID)+CURTXY(NI,2))* INVSPHTRANS(IP,2) ! upwind indicators LAMBDA(1) = ONESIXTH * SUM(CX) LAMBDA(2) = ONESIXTH * SUM(CY) ! flux jacobians KELEM(1) = MAX(ZERO, LAMBDA(1) * IEN(1,IE) + LAMBDA(2) * IEN(2,IE) )! K KELEM(2) = MAX(ZERO, LAMBDA(1) * IEN(3,IE) + LAMBDA(2) * IEN(4,IE) ) KELEM(3) = MAX(ZERO, LAMBDA(1) * IEN(5,IE) + LAMBDA(2) * IEN(6,IE) ) KSUM(IP) = KSUM(IP) + KELEM(IPOS) !2do check if also stable when abs removed IF ( KELEM(IPOS) > KMAX(IP) ) KMAX(IP) = KELEM(IPOS) END DO END DO #ifdef MPI_PARALL_GRID DTMAX_GLOBAL_EXP = VERYLARGE DTMAX_GLOBAL_EXP_LOC = VERYLARGE DO IP = 1, NP_RES DTMAX_EXP = SI(IP)/MAX(THR,KSUM(IP)) DTMAX_GLOBAL_EXP_LOC = MIN ( DTMAX_GLOBAL_EXP_LOC, DTMAX_EXP) END DO CALL MPI_ALLREDUCE(DTMAX_GLOBAL_EXP_LOC,DTMAX_GLOBAL_EXP,1,rtype,MPI_MIN,COMM,IERR) #else !2do pack it in the loop above ... DTMAX_GLOBAL_EXP = VERYLARGE DO IP = 1, MNP DTMAX_EXP = SI(IP)/MAX(THR,KSUM(IP)) !DTMAX_EXP = SI(IP)/MAX(THR,KMAX(IP)) DTMAX_GLOBAL_EXP = MIN ( DTMAX_GLOBAL_EXP, DTMAX_EXP) END DO #endif CFLXY = DT4A/DTMAX_GLOBAL_EXP REST = ABS(MOD(CFLXY,1._rkind)) IF (REST .LT. THR) THEN ITER_EXP(IS,ID) = ABS(NINT(CFLXY)) ELSE IF (REST .GT. THR .AND. REST .LT. ONEHALF) THEN ITER_EXP(IS,ID) = ABS(NINT(CFLXY)) + 1 ELSE ITER_EXP(IS,ID) = ABS(NINT(CFLXY)) END IF END DO ! IS END DO ! ID WRITE(STAT%FHNDL,*) 'MAX. ITERATIONS USED IN ADV. SCHEME', ITER_MAX, MAXVAL(ITER_EXP) FLUSH(STAT%FHNDL) END IF ! LCALC #ifdef MPI_PARALL_GRID CALL EXCHANGE_P4D_WWM(AC2) #endif #ifdef TIMINGS CALL WAV_MY_WTIME(TIME1) #endif ITER_MAX = MAXVAL(ITER_EXP) DT4AI = DT4A/ITER_MAX DO IT = 1, ITER_MAX DO ID = 1, MDC DO IS = 1, MSC UIP = AC2(IS,ID,:) !!$OMP DO PRIVATE(IP,I,IE,IPOS) DO IP = 1, MNP ST = ZERO DO I = 1, CCON(IP) ! get element and the position of IP in the element index IE = IE_CELL2(IP,I) IPOS = POS_CELL2(IP,I) ! get node indices from the element table ... NI = INE(:,IE) I1 = NI(1) I2 = NI(2) I3 = NI(3) ! estimate speed in WAE CX = (CG(NI,IS)*COSTH(ID)+CURTXY(NI,1))*INVSPHTRANS(IP,1) CY = (CG(NI,IS)*SINTH(ID)+CURTXY(NI,2))*INVSPHTRANS(IP,2) ! flux integration using simpson rule ... FL11 = CX(2) * IEN(1,IE) + CY(2) * IEN(2,IE) FL12 = CX(3) * IEN(1,IE) + CY(3) * IEN(2,IE) FL21 = CX(3) * IEN(3,IE) + CY(3) * IEN(4,IE) FL22 = CX(1) * IEN(3,IE) + CY(1) * IEN(4,IE) FL31 = CX(1) * IEN(5,IE) + CY(1) * IEN(6,IE) FL32 = CX(2) * IEN(5,IE) + CY(2) * IEN(6,IE) FL111 = TWO*FL11+FL12 FL112 = TWO*FL12+FL11 FL211 = TWO*FL21+FL22 FL212 = TWO*FL22+FL21 FL311 = TWO*FL31+FL32 FL312 = TWO*FL32+FL31 ! upwind indicators LAMBDA(1) = ONESIXTH * SUM(CX) LAMBDA(2) = ONESIXTH * SUM(CY) ! flux jacobians KELEM(1) = LAMBDA(1) * IEN(1,IE) + LAMBDA(2) * IEN(2,IE) ! K KELEM(2) = LAMBDA(1) * IEN(3,IE) + LAMBDA(2) * IEN(4,IE) KELEM(3) = LAMBDA(1) * IEN(5,IE) + LAMBDA(2) * IEN(6,IE) ! inverse of the positive sum ... N = -ONE/MIN(-THR,SUM(MIN(ZERO,KELEM))) ! N ! positive flux jacobians KELEM(1) = MAX(ZERO,KELEM(1)) ! K+ KELEM(2) = MAX(ZERO,KELEM(2)) KELEM(3) = MAX(ZERO,KELEM(3)) ! simposon integration last step ... FLALL(1) = (FL311 + FL212) * ONESIXTH + KELEM(1) FLALL(2) = (FL111 + FL312) * ONESIXTH + KELEM(2) FLALL(3) = (FL211 + FL112) * ONESIXTH + KELEM(3) ! flux conserving upwind contribution U3 = UIP(NI) UTILDE3 = N * ( FLALL(1) * U3(1) + FLALL(2) * U3(2) + FLALL(3) * U3(3) ) ! coefficient for the integration in time ST = ST + KELEM(IPOS) * (UIP(IP) - UTILDE3) END DO ! time stepping ... UIP(IP) = MAX(ZERO,UIP(IP)-DT4AI/SI(IP)*ST*IOBWB(IP))*IOBPD(ID,IP) END DO !IP AC2(IS,ID,:) = UIP END DO !ID END DO !IS #ifdef MPI_PARALL_GRID CALL EXCHANGE_P4D_WWM(AC2) #endif END DO !IT END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE EXPLICIT_N_SCHEME_VECTOR_HPCF_VECTOPER USE DATAPOOL IMPLICIT NONE ! ! local integer ! INTEGER :: IP, IE, IT, IS, ID INTEGER :: I1, I2, I3 INTEGER :: NI(3), I, IPOS ! ! local double ! REAL(rkind) :: DTMAX_GLOBAL_EXP(MSC,MDC), DTMAX_EXP REAL(rkind) :: DTMAX_GLOBAL_EXP_LOC(MSC,MDC) REAL(rkind) :: REST, CFLXY REAL(rkind) :: LAMBDA(MSC,MDC,2), DT4AI REAL(rkind) :: FL11(MSC,MDC), FL12(MSC,MDC) REAL(rkind) :: FL21(MSC,MDC), FL22(MSC,MDC) REAL(rkind) :: FL31(MSC,MDC), FL32(MSC,MDC) REAL(rkind) :: ST(MSC,MDC), N(MSC,MDC) REAL(rkind) :: CX(MSC,MDC,3), CY(MSC,MDC,3) REAL(rkind) :: FLALL(MSC,MDC,3) REAL(rkind) :: KELEM(MSC,MDC,3) REAL(rkind) :: KM(MSC,MDC,3), KP(MSC,MDC,3) REAL(rkind) :: FL111(MSC,MDC), FL112(MSC,MDC) REAL(rkind) :: FL211(MSC,MDC), FL212(MSC,MDC) REAL(rkind) :: FL311(MSC,MDC), FL312(MSC,MDC) REAL(rkind) :: UTILDE3(MSC,MDC) REAL(rkind) :: KSUM(MSC,MDC), KMAX(MSC,MDC) REAL(rkind) :: TIME1, TIME2 ! ! local parameter ! ! Calculate phase speeds for the certain spectral component ... ! IF (LCALC) THEN DTMAX_GLOBAL_EXP_LOC=VERYLARGE DO IP = 1, MNP KMAX = ZERO KSUM = ZERO DO I = 1, CCON(IP) IE = IE_CELL2(IP,I) IPOS = POS_CELL2(IP,I) ! get node indices from the element table ... NI = INE(:,IE) ! estimate speed in WAE DO ID = 1, MDC DO IS = 1, MSC CX(IS,ID,:) = (CG(IS,NI)*COSTH(ID)+CURTXY(NI,1))* INVSPHTRANS(IP,1) CY(IS,ID,:) = (CG(IS,NI)*SINTH(ID)+CURTXY(NI,2))* INVSPHTRANS(IP,2) END DO END DO ! upwind indicators LAMBDA(:,:,1) = ONESIXTH * (CX(:,:,1) + CX(:,:,1) + CX(:,:,1)) LAMBDA(:,:,2) = ONESIXTH * (CY(:,:,1) + CY(:,:,1) + CY(:,:,1)) ! flux jacobians KELEM(:,:,1) = MAX(ZERO, LAMBDA(:,:,1) * IEN(1,IE) + LAMBDA(:,:,2) * IEN(2,IE) )! K KELEM(:,:,2) = MAX(ZERO, LAMBDA(:,:,1) * IEN(3,IE) + LAMBDA(:,:,2) * IEN(4,IE) ) KELEM(:,:,3) = MAX(ZERO, LAMBDA(:,:,1) * IEN(5,IE) + LAMBDA(:,:,2) * IEN(6,IE) ) KSUM(:,:) = KSUM(:,:) + KELEM(:,:,IPOS) !2do check if also stable when abs removed IF (IP .le. NP_RES) THEN DO ID=1,MDC DO IS=1,MSC IF ( KELEM(IS,ID,IPOS) > KMAX(IS,ID) ) KMAX(IS,ID) = KELEM(IS,ID,IPOS) DTMAX_EXP=SI(IP)/MAX(THR,KSUM(IS,ID)) !DTMAX_EXP=SI(IP)/MAX(THR,KMAX(IS,ID)) DTMAX_GLOBAL_EXP_LOC(IS,ID)=MIN(DTMAX_GLOBAL_EXP_LOC(IS,ID),DTMAX_EXP) END DO END DO END IF #ifdef MPI_PARALL_GRID CALL MPI_ALLREDUCE(DTMAX_GLOBAL_EXP_LOC,DTMAX_GLOBAL_EXP,MSC*MDC,rtype,MPI_MIN,COMM,IERR) #else DTMAX_GLOBAL_EXP=DTMAX_GLOBAL_EXP_LOC #endif DO ID=1,MDC DO IS=1,MSC CFLXY = DT4A/DTMAX_GLOBAL_EXP(IS,ID) REST = ABS(MOD(CFLXY,ONE)) IF (REST .LT. THR) THEN ITER_EXP(IS,ID) = ABS(NINT(CFLXY)) ELSE IF (REST .GT. THR .AND. REST .LT. ONEHALF) THEN ITER_EXP(IS,ID) = ABS(NINT(CFLXY)) + 1 ELSE ITER_EXP(IS,ID) = ABS(NINT(CFLXY)) END IF END DO ! IS END DO ! ID WRITE(STAT%FHNDL,*) 'MAX. ITERATIONS USED IN ADV. SCHEME', ITER_MAX, MAXVAL(ITER_EXP) FLUSH(STAT%FHNDL) END DO END DO END IF #ifdef TIMINGS CALL WAV_MY_WTIME(TIME1) #endif ITER_MAX = MAXVAL(ITER_EXP) DT4AI = DT4A/ITER_MAX DO IT = 1, ITER_MAX AC1 = AC2 DO IP = 1, MNP ST = ZERO DO I = 1, CCON(IP) ! get element and the position of IP in the element index IE = IE_CELL2(IP,I) IPOS = POS_CELL2(IP,I) ! get node indices from the element table ... NI = INE(:,IE) I1 = NI(1) I2 = NI(2) I3 = NI(3) ! estimate speed in WAE DO ID=1,MDC DO IS=1,MSC CX(IS,ID,:) = (CG(IS,NI)*COSTH(ID)+CURTXY(NI,1))*INVSPHTRANS(IP,1) CY(IS,ID,:) = (CG(IS,NI)*SINTH(ID)+CURTXY(NI,2))*INVSPHTRANS(IP,2) END DO END DO ! flux integration using simpson rule ... FL11(:,:) = CX(:,:,2) * IEN(1,IE) + CY(:,:,2) * IEN(2,IE) FL12(:,:) = CX(:,:,3) * IEN(1,IE) + CY(:,:,3) * IEN(2,IE) FL21(:,:) = CX(:,:,3) * IEN(3,IE) + CY(:,:,3) * IEN(4,IE) FL22(:,:) = CX(:,:,1) * IEN(3,IE) + CY(:,:,1) * IEN(4,IE) FL31(:,:) = CX(:,:,1) * IEN(5,IE) + CY(:,:,1) * IEN(6,IE) FL32(:,:) = CX(:,:,2) * IEN(5,IE) + CY(:,:,2) * IEN(6,IE) FL111(:,:) = TWO*FL11(:,:)+FL12(:,:) FL112(:,:) = TWO*FL12(:,:)+FL11(:,:) FL211(:,:) = TWO*FL21(:,:)+FL22(:,:) FL212(:,:) = TWO*FL22(:,:)+FL21(:,:) FL311(:,:) = TWO*FL31(:,:)+FL32(:,:) FL312(:,:) = TWO*FL32(:,:)+FL31(:,:) ! upwind indicators LAMBDA(:,:,1) = ONESIXTH * (CX(:,:,1) + CX(:,:,2) + CX(:,:,3)) LAMBDA(:,:,2) = ONESIXTH * (CY(:,:,1) + CY(:,:,2) + CY(:,:,3)) ! flux jacobians KELEM(:,:,1) = LAMBDA(:,:,1) * IEN(1,IE) + LAMBDA(:,:,2) * IEN(2,IE) ! K KELEM(:,:,2) = LAMBDA(:,:,1) * IEN(3,IE) + LAMBDA(:,:,2) * IEN(4,IE) KELEM(:,:,3) = LAMBDA(:,:,1) * IEN(5,IE) + LAMBDA(:,:,2) * IEN(6,IE) ! inverse of the positive sum ... KM=MIN(ZERO,KELEM) KP=MAX(ZERO,KELEM) N(:,:) = -ONE/MIN(-THR,KM(:,:,1) + KM(:,:,2) + KM(:,:,3)) ! simposon integration last step ... FLALL(:,:,1) = (FL311(:,:) + FL212(:,:)) * ONESIXTH + KP(:,:,1) FLALL(:,:,2) = (FL111(:,:) + FL312(:,:)) * ONESIXTH + KP(:,:,2) FLALL(:,:,3) = (FL211(:,:) + FL112(:,:)) * ONESIXTH + KP(:,:,3) ! flux conserving upwind contribution UTILDE3(:,:) = N(:,:) * ( FLALL(:,:,1) * AC2(:,:,I1) + FLALL(:,:,2) * AC2(:,:,I2) + FLALL(:,:,3) * AC2(:,:,3) ) ! coefficient for the integration in time ST(:,:) = ST(:,:) + KP(:,:,IPOS) * (AC2(:,:,IP) - UTILDE3(:,:)) ! time stepping ... DO ID=1,MDC AC1(:,ID,IP) = MAX(ZERO,AC1(:,ID,IP)-DT4AI/SI(IP)*ST(:,ID)*IOBWB(IP))*IOBPD(ID,IP) END DO END DO END DO AC2 = AC1 #ifdef MPI_PARALL_GRID CALL EXCHANGE_P4D_WWM(AC2) #endif END DO END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE EXPLICIT_N_SCHEME_HPCF2 USE DATAPOOL IMPLICIT NONE ! ! local integer ! INTEGER :: IP, IE, IT, IS, ID INTEGER :: I1, I2, I3 INTEGER :: NI(3), I, IPOS ! ! local double ! REAL(rkind) :: DTMAX_GLOBAL_EXP, DTMAX_EXP REAL(rkind) :: DTMAX_GLOBAL_EXP_LOC REAL(rkind) :: REST, CFLXY REAL(rkind) :: LLAMBDA(2), DT4AI REAL(rkind) :: FL11,FL12,FL21,FL22,FL31,FL32 REAL(rkind) :: KTMP(3) REAL(rkind) :: ST, N, KSUM, KMAX REAL(rkind) :: CX(MNP), CY(MNP)! 20000 * 8 / 1024**2 REAL(rkind) :: FLALL(3), TMP REAL(rkind) :: KELEM(3,MNE), KKELEM(3) ! 3 * 20000 / 1024**2 REAL(rkind) :: FL111, FL112, FL211, FL212, FL311, FL312 REAL(rkind) :: UTILDE3(MNE) ! 20000 * 8 / 1024**2. REAL(rkind) :: TIME1, TIME2 ! ! Calculate K-Values and contour based quantities ... ! IF (LCALC) THEN DO ID = 1, MDC DO IS = 1, MSC CX = (CG(IS,:)*COSTH(ID)+CURTXY(:,1))*INVSPHTRANS(:,1) CY = (CG(IS,:)*SINTH(ID)+CURTXY(:,2))*INVSPHTRANS(:,2) DTMAX_GLOBAL_EXP = VERYLARGE DTMAX_GLOBAL_EXP_LOC = VERYLARGE DO IP = 1, MNP KSUM = ZERO KMAX = ZERO DO I = 1, CCON(IP) IE = IE_CELL2(IP,I) IPOS = POS_CELL2(IP,I) ! get node indices from the element table ... NI = INE(:,IE) ! upwind indicators LLAMBDA(1) = ONESIXTH * SUM(CX(NI)) LLAMBDA(2) = ONESIXTH * SUM(CY(NI)) ! flux jacobians KKELEM(1) = MAX(ZERO, LLAMBDA(1) * IEN(1,IE) + LLAMBDA(2) * IEN(2,IE) )! K KKELEM(2) = MAX(ZERO, LLAMBDA(1) * IEN(3,IE) + LLAMBDA(2) * IEN(4,IE) ) KKELEM(3) = MAX(ZERO, LLAMBDA(1) * IEN(5,IE) + LLAMBDA(2) * IEN(6,IE) ) ! sum over connected nodes KSUM = KSUM + KKELEM(IPOS) !2do check if also stable when abs removed IF ( KKELEM(IPOS) > KMAX ) KMAX = KKELEM(IPOS) END DO DTMAX_EXP = SI(IP)/MAX(THR,KSUM) !DTMAX_EXP = SI(IP)/MAX(THR,KMAX) #ifdef MPI_PARALL_GRID DTMAX_GLOBAL_EXP_LOC = MIN ( DTMAX_GLOBAL_EXP_LOC, DTMAX_EXP) #else DTMAX_GLOBAL_EXP = MIN ( DTMAX_GLOBAL_EXP, DTMAX_EXP) #endif END DO #ifdef MPI_PARALL_GRID CALL MPI_ALLREDUCE(DTMAX_GLOBAL_EXP_LOC,DTMAX_GLOBAL_EXP,1,rtype,MPI_MIN,COMM,IERR) #endif CFLXY = DT4A/DTMAX_GLOBAL_EXP REST = ABS(MOD(CFLXY,ONE)) IF (REST .LT. THR) THEN ITER_EXP(IS,ID) = ABS(NINT(CFLXY)) ELSE IF (REST .GT. THR .AND. REST .LT. ONEHALF) THEN ITER_EXP(IS,ID) = ABS(NINT(CFLXY)) + 1 ELSE ITER_EXP(IS,ID) = ABS(NINT(CFLXY)) END IF END DO ! IS END DO ! ID ITER_MAX = MAXVAL(ITER_EXP) WRITE(STAT%FHNDL,*) 'MAX. ITERATIONS USED IN ADV. SCHEME', ITER_MAX FLUSH(STAT%FHNDL) END IF #ifdef TIMINGS CALL WAV_MY_WTIME(TIME1) #endif DT4AI = DT4A/ITER_MAX DO IT = 1, ITER_MAX DO ID = 1, MSC DO IS = 1, MDC !!$OMP WORKSHARE CX = (CG(IS,:)*COSTH(ID)+CURTXY(:,1))*INVSPHTRANS(:,1) CY = (CG(IS,:)*SINTH(ID)+CURTXY(:,2))*INVSPHTRANS(:,2) !!$OMP END WORKSHARE !!$OMP DO PRIVATE(IE,NI,I1,I2,I3,LLAMBDA,FL11,FL12,FL21,FL22,FL31,FL32,FL111,FL112,FL211,FL212,FL311,FL312,FLALL) DO IE = 1, MNE ! must go over the augmented domain ... NI = INE(:,IE) I1 = INE(1,IE) I2 = INE(2,IE) I3 = INE(3,IE) LLAMBDA(1) = ONESIXTH *(CX(I1)+CX(I2)+CX(I3)) LLAMBDA(2) = ONESIXTH *(CY(I1)+CY(I2)+CY(I3)) KELEM(1,IE) = LLAMBDA(1) * IEN(1,IE) + LLAMBDA(2) * IEN(2,IE) KELEM(2,IE) = LLAMBDA(1) * IEN(3,IE) + LLAMBDA(2) * IEN(4,IE) KELEM(3,IE) = LLAMBDA(1) * IEN(5,IE) + LLAMBDA(2) * IEN(6,IE) KTMP(1) = KELEM(1,IE) KTMP(2) = KELEM(2,IE) KTMP(3) = KELEM(3,IE) TMP = SUM(MIN(ZERO,KTMP)) N = -1._rkind/MIN(-THR,TMP) KELEM(1,IE) = MAX(ZERO,KTMP(1)) KELEM(2,IE) = MAX(ZERO,KTMP(2)) KELEM(3,IE) = MAX(ZERO,KTMP(3)) ! WRITE(DBG%FHNDL,'(3I10,3F15.4)') IS, ID, IE, KELEM(:,IE) FL11 = CX(I2) * IEN(1,IE) + CY(I2) * IEN(2,IE) FL12 = CX(I3) * IEN(1,IE) + CY(I3) * IEN(2,IE) FL21 = CX(I3) * IEN(3,IE) + CY(I3) * IEN(4,IE) FL22 = CX(I1) * IEN(3,IE) + CY(I1) * IEN(4,IE) FL31 = CX(I1) * IEN(5,IE) + CY(I1) * IEN(6,IE) FL32 = CX(I2) * IEN(5,IE) + CY(I2) * IEN(6,IE) FL111 = TWO*FL11+FL12 FL112 = TWO*FL12+FL11 FL211 = TWO*FL21+FL22 FL212 = TWO*FL22+FL21 FL311 = TWO*FL31+FL32 FL312 = TWO*FL32+FL31 FLALL(1) = (FL311 + FL212) * ONESIXTH + KELEM(1,IE) FLALL(2) = (FL111 + FL312) * ONESIXTH + KELEM(2,IE) FLALL(3) = (FL211 + FL112) * ONESIXTH + KELEM(3,IE) UTILDE3(IE) = N * (DOT_PRODUCT(FLALL,AC2(IS,ID,NI(:)))) !UTILDE3(IE) = N * ( FLALL(1) * AC2(IS,ID,NI(1)) + FLALL(2) * AC2(IS,ID,NI(2) + FLALL(3) * AC2(IS,ID,NI(3))) END DO !IE !!$OMP DO PRIVATE(IP,ST,IE,IPOS) DO IP = 1, NP_RES ST = 0 DO I = 1, CCON(IP) IE = IE_CELL2(IP,I) IPOS = POS_CELL2(IP,I) ST = ST + KELEM(IPOS,IE) * (AC2(IS,ID,IP) - UTILDE3(IE)) END DO AC2(IS,ID,IP) = MAX(ZERO,AC2(IS,ID,IP)-DT4AI/SI(IP)*ST*IOBWB(IP))*IOBPD(ID,IP) END DO !IP END DO !IS END DO !ID #ifdef MPI_PARALL_GRID CALL EXCHANGE_P4D_WWM(AC2) #endif END DO END SUBROUTINE !********************************************************************** !* * !**********************************************************************