#include "wwm_functions.h" !********************************************************************** !* * !********************************************************************** SUBROUTINE GRADDEP() USE DATAPOOL IMPLICIT NONE INTEGER :: IP REAL(rkind) :: GDL, GDD DDEP(:,:) = 0.0 SELECT CASE (DIMMODE) CASE (1) #ifdef MPI_PARALL_GRID call parallel_abort('WWM - 1d mode cannot work with SCHISM ') #endif CALL DIFFERENTIATE_XDIR(DEP,DDEP(:,1)) IF (LSPHE) THEN DO IP = 1, MNP DDEP(IP,1) = DDEP(IP,1)/( DEGRAD*REARTH*COS(YP(IP)*DEGRAD) ) END DO END IF IF (LTEST .AND. ITEST > 100) THEN WRITE(STAT%FHNDL,*) 'Gradients of depth ' WRITE(STAT%FHNDL,*) '@D/@X ' DO IP = 1, MNP WRITE(STAT%FHNDL,'(1X,I5,3F10.5)') IP, DDEP(IP,1) END DO FLUSH(STAT%FHNDL) END IF CASE (2) CALL DIFFERENTIATE_XYDIR(DEP,DDEP(:,1),DDEP(:,2)) IF (LSPHE) THEN DO IP = 1, MNP DDEP(IP,1) = DDEP(IP,1)/( DEGRAD*REARTH*COS(YP(IP)*DEGRAD) ) DDEP(IP,2) = DDEP(IP,2)/( DEGRAD*REARTH ) END DO END IF IF (LSLOP) THEN DO IP = 1, MNP GDL = SQRT(DDEP(IP,1)**2 + DDEP(IP,2)**2) ! Achtung Gradient mit SQRT berechnet ... GDD = MyATAN2(DDEP(IP,2), DDEP(IP,1)) IF (GDL < 1.0E-8) CYCLE GDL = SQRT(GDL) IF (GDL > SLMAX) THEN DDEP(IP,1) = SLMAX*COS(GDD) DDEP(IP,2) = SLMAX*SIN(GDD) WRITE(STAT%FHNDL,*) IP, SLMAX, GDL, GDD , 'MAXSLOPE' END IF END DO FLUSH(STAT%FHNDL) END IF CASE DEFAULT END SELECT END SUBROUTINE !********************************************************************** !* * !********************************************************************** #ifdef TIMINGS SUBROUTINE WAV_MY_WTIME(wtime) USE DATAPOOL, only : rkind implicit none real(rkind), intent(inout) :: wtime # ifdef MPI_PARALL_GRID real(8) mpi_wtime wtime=mpi_wtime() # else CALL cpu_time(wtime) # endif END SUBROUTINE #endif !********************************************************************** !* * !********************************************************************** SUBROUTINE GRADCURT() USE DATAPOOL IMPLICIT NONE INTEGER :: IP DCUX(:,:) = 0.0 DCUY(:,:) = 0.0 SELECT CASE (DIMMODE) CASE (1) CALL DIFFERENTIATE_XDIR(CURTXY(:,1),DCUX(:,1)) CALL DIFFERENTIATE_XDIR(CURTXY(:,2),DCUY(:,1)) IF (LSPHE) THEN DO IP = 1, MNP DCUX(IP,1) = DCUX(IP,1)/( DEGRAD*REARTH*COS(YP(IP)*DEGRAD) ) DCUY(IP,1) = DCUY(IP,1)/( DEGRAD*REARTH*COS(YP(IP)*DEGRAD) ) END DO END IF IF (LTEST .AND. ITEST > 100) THEN WRITE(STAT%FHNDL,*) 'Gradients of depth and current' WRITE(STAT%FHNDL,*) '@U/@X @V/@X' DO IP = 1, MNP WRITE(STAT%FHNDL,'(1X,I5,3F10.5)') IP, DCUX, DCUY END DO FLUSH(STAT%FHNDL) END IF CASE (2) CALL DIFFERENTIATE_XYDIR(CURTXY(:,1),DCUX(:,1),DCUX(:,2)) CALL DIFFERENTIATE_XYDIR(CURTXY(:,2),DCUY(:,1),DCUY(:,2)) IF (LSPHE) THEN DO IP = 1, MNP DCUX(IP,1) = DCUX(IP,1)/( DEGRAD*REARTH*COS(YP(IP)*DEGRAD) ) DCUY(IP,1) = DCUY(IP,1)/( DEGRAD*REARTH*COS(YP(IP)*DEGRAD) ) DCUX(IP,2) = DCUX(IP,2)/( DEGRAD*REARTH ) DCUY(IP,2) = DCUY(IP,2)/( DEGRAD*REARTH ) END DO END IF IF (LTEST .AND. ITEST > 100) THEN WRITE(STAT%FHNDL,*) 'The Gradient of Depth and Current' WRITE(STAT%FHNDL,*) ' @U/@X @U/@Y @V/@X @V/@Y' DO IP = 1, MNP WRITE(STAT%FHNDL,'(1X,I5,4F15.7)') IP, DCUX(IP,1), DCUX(IP,2), DCUY(IP,1), DCUY(IP,2) END DO FLUSH(STAT%FHNDL) END IF CASE DEFAULT END SELECT END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE GRAD_CG_K() USE DATAPOOL IMPLICIT NONE INTEGER :: IS DCGDX(:,:) = 0.0 DCGDY(:,:) = 0.0 DWKDX(:,:) = 0.0 DWKDY(:,:) = 0.0 SELECT CASE (DIMMODE) CASE (1) DO IS = 1, MSC CALL DIFFERENTIATE_XDIR(CG(IS,:),DCGDX(:,IS)) CALL DIFFERENTIATE_XDIR(WK(IS,:),DWKDX(:,IS)) ENDDO CASE (2) DO IS = 1, MSC CALL DIFFERENTIATE_XYDIR(CG(IS,:),DCGDX(:,IS),DCGDY(:,IS)) CALL DIFFERENTIATE_XYDIR(WK(IS,:),DWKDX(:,IS),DWKDY(:,IS)) ENDDO CASE DEFAULT END SELECT END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE SMOOTH( BETA, MNP, XP, YP, VAR ) USE DATAPOOL, ONLY : RKIND, ZERO, ONE, TWO IMPLICIT NONE INTEGER :: MNP REAL(rkind), INTENT(IN) :: XP(MNP), YP(MNP) REAL(rkind), INTENT(INOUT) :: VAR(MNP) REAL(rkind), INTENT(IN) :: BETA REAL(rkind) :: VART(MNP) REAL(rkind) :: SW, SWQ, DISX, DISY, DIST, DIS INTEGER :: I, J DO I = 1, MNP SW = ZERO SWQ = ZERO DO J = 1, MNP DISX = (XP(I) - XP(J))**2 DISY = (YP(I) - YP(J))**2 DIST = DISX + DISY IF (DIST > TINY(1.)) THEN DIS = SQRT(DIST)**BETA ELSE DIS = ONE END IF SW = SW + DIS SWQ = SWQ + DIS*VAR(J) END DO VART(I) = SWQ / SW END DO VAR(:) = VART(:) END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE SMOOTH_V2(VAR_IN, VAR_OUT) USE DATAPOOL IMPLICIT NONE REAL(rkind), INTENT(IN) :: VAR_IN(MNP) REAL(rkind), INTENT(OUT) :: VAR_OUT(MNP) INTEGER IP, IADJ, IP_ADJ REAL(rkind) :: SumVAR, SumSI, eVal DO IP = 1, NP_RES SumVAR=SI(IP) * VAR_IN(IP) SumSI =SI(IP) DO IADJ=1,VERT_DEG(IP) IP_ADJ=LIST_ADJ_VERT(IADJ,IP) SumVAR=SumVAR + SI(IP_ADJ)*VAR_IN(IP_ADJ) SumSI =SumSI + SI(IP_ADJ) END DO eVal=SumVAR/SumSI VAR_OUT(IP)=eVal END DO #ifdef MPI_PARALL_GRID CALL exchange_p2d(VAR_OUT) #endif END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE SMOOTH_ON_TRIANGLE(VAR_IN, VAR_OUT, nbIter) USE DATAPOOL IMPLICIT NONE integer, intent(in) :: nbIter REAL(rkind), INTENT(IN) :: VAR_IN(MNE) REAL(rkind), INTENT(OUT) :: VAR_OUT(MNE) INTEGER IE, IEadj, I, nb, iIter REAL(rkind) :: eVal, eValB REAL(rkind) :: VARextent(MNEextent) VAR_OUT=VAR_IN DO iIter=1,nbIter VARextent(1:MNE)=VAR_OUT #ifdef MPI_PARALL_GRID CALL TRIG_SYNCHRONIZATION(VARextent) #endif DO IE=1,MNE eVal=ZERO nb=0 DO I=1,3 IEadj=IEneighbor(I,IE) IF (IEadj .gt. 0) THEN eVal=eVal + VARextent(IEadj) nb=nb+1 END IF END DO eValB=(MyREAL(6-nb)*VARextent(IE) + eVal)/6.0_rkind VAR_OUT(IE)=eValB END DO END DO END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE DIFFERENTIATE_XDIR(VAR, DVDX) USE DATAPOOL IMPLICIT NONE REAL(rkind), INTENT(IN) :: VAR(MNP) REAL(rkind), INTENT(OUT) :: DVDX(MNP) INTEGER :: IP REAL(rkind) :: TMP1, TMP2 DVDX(1) = (VAR(2)-VAR(1))/(XP(2)-XP(1)) DVDX(MNP) = (VAR(MNP)-VAR(MNP-1))/(XP(MNP)-XP(MNP-1)) DO IP = 2, MNP-1 ! DVDX(IP) = ( VAR(IP)-VAR(IP-1) ) / ( XP(IP)-XP(IP-1) ) ! DVDX(IP) = (VAR(IP+1)-VAR(IP))/(XP(IP+1)-XP(IP)) ! TMP1 = (VAR(IP)-VAR(IP-1))/(XP(IP)-XP(IP-1)) ! TMP2 = (VAR(IP+1)-VAR(IP))/(XP(IP+1)-XP(IP)) ! DVDX(IP) = 0.5 * (TMP1 + TMP2) TMP1 = VAR(IP+1)-VAR(IP-1) TMP2 = XP(IP+1)-XP(IP-1) DVDX(IP) = TMP1/TMP2 END DO END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE DIFFERENTIATE_XYDIR(VAR, DVDX, DVDY) USE DATAPOOL IMPLICIT NONE REAL(rkind), INTENT(IN) :: VAR(MNP) REAL(rkind), INTENT(OUT) :: DVDX(MNP), DVDY(MNP) INTEGER :: NI(3) INTEGER :: IE, I1, I2, I3, IP REAL(rkind) :: DEDY(3),DEDX(3) REAL(rkind) :: DVDXIE, DVDYIE REAL(rkind) :: WEI(MNP) WEI(:) = 0.0_rkind DVDX(:) = 0.0_rkind DVDY(:) = 0.0_rkind DO IE = 1, MNE NI = INE(:,IE) I1 = INE(1,IE) I2 = INE(2,IE) I3 = INE(3,IE) WEI(NI) = WEI(NI) + 2.*TRIA(IE) DEDX(1) = IEN(1,IE) DEDX(2) = IEN(3,IE) DEDX(3) = IEN(5,IE) DEDY(1) = IEN(2,IE) DEDY(2) = IEN(4,IE) DEDY(3) = IEN(6,IE) DVDXIE = DOT_PRODUCT( VAR(NI),DEDX) DVDYIE = DOT_PRODUCT( VAR(NI),DEDY) DVDX(NI) = DVDX(NI) + DVDXIE DVDY(NI) = DVDY(NI) + DVDYIE END DO DVDX(:) = DVDX(:)/WEI(:) DVDY(:) = DVDY(:)/WEI(:) DO IP = 1, MNP IF (DEP(IP) .LT. DMIN) THEN DVDX(IP) = 0. DVDY(IP) = 0. END IF END DO #ifdef MPI_PARALL_GRID CALL exchange_p2d(DVDX) CALL exchange_p2d(DVDY) #endif IF (.FALSE.) THEN OPEN(2305, FILE = 'erggrad.bin' , FORM = 'UNFORMATTED') WRITE(2305) 1. WRITE(2305) (DVDX(IP), DVDY(IP), SQRT(DVDY(IP)**2+DVDY(IP)**2), IP = 1, MNP) ENDIF END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE DIFFERENTIATE_XYDIR_NEIGHADJ(VAR, DVDX, DVDY) USE DATAPOOL IMPLICIT NONE REAL(rkind), INTENT(IN) :: VAR(MNP) REAL(rkind), INTENT(OUT) :: DVDX(MNP), DVDY(MNP) INTEGER :: IP, IADJ, IP_ADJ REAL(rkind) :: xpdiff, ypdiff, dist, eSum REAL(rkind) :: XPcomp, YPcomp, eXP, eYP, eCoeff DO IP=1,NP_RES eSum=ZERO XPcomp=ZERO YPcomp=ZERO DO IADJ=1,VERT_DEG(IP) IP_ADJ=LIST_ADJ_VERT(IADJ,IP) xpdiff=XP(IP_ADJ) - XP(IP) ypdiff=YP(IP_ADJ) - YP(IP) dist=sqrt(xpdiff*xpdiff+ypdiff*ypdiff) eSum=eSum + ONE/dist eCoeff=(VAR(IP_ADJ) - VAR(IP))/(dist**3) XPcomp = XPcomp + xpdiff*eCoeff YPcomp = YPcomp + ypdiff*eCoeff END DO eXP=XPcomp / eSum eYP=YPcomp / eSum DVDX(IP)=eXP DVDY(IP)=eYP END DO #ifdef MPI_PARALL_GRID CALL exchange_p2d(DVDX) CALL exchange_p2d(DVDY) #endif END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE DIFFERENTIATE_XYDIR_AND_SMOOTH(VAR, DVDX, DVDY, nbIter) USE DATAPOOL IMPLICIT NONE integer, intent(in) :: nbIter REAL(rkind), INTENT(IN) :: VAR(MNP) REAL(rkind), INTENT(OUT) :: DVDX(MNP), DVDY(MNP) INTEGER :: NI(3) INTEGER :: IE, I1, I2, I3, IP REAL(rkind) :: DEDY(3),DEDX(3) REAL(rkind) :: DVDXIE, DVDYIE REAL(rkind) :: WEI(MNP) REAL(rkind) :: IE_DY(MNE),IE_DX(MNE) REAL(rkind) :: IE_DY_SMO(MNE),IE_DX_SMO(MNE) WEI(:) = 0.0_rkind DO IE = 1, MNE NI = INE(:,IE) I1 = INE(1,IE) I2 = INE(2,IE) I3 = INE(3,IE) WEI(NI) = WEI(NI) + 2.*TRIA(IE) DEDX(1) = IEN(1,IE) DEDX(2) = IEN(3,IE) DEDX(3) = IEN(5,IE) DEDY(1) = IEN(2,IE) DEDY(2) = IEN(4,IE) DEDY(3) = IEN(6,IE) DVDXIE = DOT_PRODUCT( VAR(NI),DEDX) DVDYIE = DOT_PRODUCT( VAR(NI),DEDY) IE_DX(IE)=DVDXIE IE_DY(IE)=DVDYIE END DO ! ! Now smoothing ! CALL SMOOTH_ON_TRIANGLE(IE_DX, IE_DX_SMO, nbIter) CALL SMOOTH_ON_TRIANGLE(IE_DY, IE_DY_SMO, nbIter) ! ! Mapping to nodes ! DVDX(:) = 0.0_rkind DVDY(:) = 0.0_rkind DO IE = 1, MNE NI = INE(:,IE) I1 = INE(1,IE) I2 = INE(2,IE) I3 = INE(3,IE) DVDX(NI) = DVDX(NI) + IE_DX_SMO(IE) DVDY(NI) = DVDY(NI) + IE_DY_SMO(IE) END DO DVDX(:) = DVDX(:)/WEI(:) DVDY(:) = DVDY(:)/WEI(:) DO IP = 1, MNP IF (DEP(IP) .LT. DMIN) THEN DVDX(IP) = 0. DVDY(IP) = 0. END IF END DO #ifdef MPI_PARALL_GRID CALL exchange_p2d(DVDX) CALL exchange_p2d(DVDY) #endif END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE WAVEKCG(DEPLOC, SIGIN, WN, WVC, WVK, WVCG2) USE DATAPOOL IMPLICIT NONE REAL(rkind), INTENT(IN) :: DEPLOC, SIGIN REAL(rkind), INTENT(OUT) :: WVC, WVK, WVCG2, WN REAL(rkind) :: SGDLS , AUX1, AUX2 REAL(rkind) :: WKDEP ! WVC = 0. WVK = 0. WVCG2 = 0. WN = 0. IF (SIGIN .LT. VERYSMALL) THEN WN = 0. WVK=10. WVCG2=0. RETURN END IF IF (DEPLOC .GT. DMIN) THEN SGDLS = SIGIN*SIGIN*DEPLOC/G9 AUX1 = 1.0+0.6522*SGDLS+0.4622*(SGDLS**2.0)+0.0864*(SGDLS**4.0)+0.0675*(SGDLS**5.0) AUX2 = 1.0/(SGDLS+1.0/AUX1) WVC = SQRT(AUX2*G9*DEPLOC) WVK = SIGIN/WVC WKDEP = WVK*DEPLOC IF (WKDEP > 13.0_rkind) THEN WN = 0.5_rkind ELSE WN = 0.5_rkind*(ONE+TWO*WKDEP/SINH(MIN(TWO*KDMAX,TWO*WKDEP))) END IF WVCG2 = WN*WVC ELSE WVC = 0. WVK = 10. WVCG2 = 0. END IF END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE MAKE_WAVE_TABLE() USE DATAPOOL IMPLICIT NONE INTEGER :: IS REAL(rkind) :: SIGIN, WVN, WVC, WVK, WVCG REAL(rkind) :: DEPTH, SIGMAMAX NMAX = IDISPTAB - 1 DEPTH = 1. SIGMAMAX = SQRT (G9 * DEPFAC) DSIGTAB = SIGMAMAX / MyREAL(NMAX) TABK(0) = 0. TABCG(0) = SQRT(G9) DO IS = 1, NMAX SIGIN = MyREAL(IS)*DSIGTAB CALL WAVEKCG(DEPTH, SIGIN, WVN, WVC, WVK, WVCG) TABK(IS) = WVK TABCG(IS) = WVCG END DO IS = NMAX + 1 SIGIN = MyREAL(IS)*DSIGTAB CALL WAVEKCG(DEPTH, SIGIN, WVN, WVC, WVK, WVCG) TABK(IS) = WVK TABCG(IS) = WVCG END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE ALL_FROM_TABLE(SIGIN,DEPTH,WVK,WVCG,WVKDEP,WVN,WVC) USE DATAPOOL IMPLICIT NONE REAL(rkind), INTENT(IN) :: SIGIN, DEPTH REAL(rkind), INTENT(OUT) :: WVK, WVCG, WVKDEP, WVN, WVC REAL(rkind) :: SQRTDEP, SIGSQDEP, R1, R2, DEPLOC INTEGER :: ITAB, ITAB2 DEPLOC = MAX(DMIN,DEPTH) SQRTDEP = MySQRT(DEPLOC) SIGSQDEP = SIGIN * SQRTDEP ITAB = INT(SIGSQDEP/DSIGTAB) ! IF (ITAB.LE.NMAX.AND.ITAB.GE.0) THEN ITAB2 = ITAB + 1 R1 = SIGSQDEP/DSIGTAB - MyREAL(ITAB) R2 = ONE - R1 WVK = ( R2*TABK(ITAB) + R1*TABK(ITAB2) ) / DEPLOC WVCG = ( R2*TABCG(ITAB) + R1*TABCG(ITAB2) ) * SQRTDEP WVKDEP = WVK * DEPLOC WVN = ZEROFIVE*(ONE+TWO*WVKDEP/MySINH(MIN(KDMAX,TWO*WVKDEP))) WVC = WVCG/WVN ELSE WVK = SIGIN*SIGIN/G9 WVCG = ZEROFIVE*G9/SIGIN WVKDEP = WVK * DEPLOC WVN = ZEROFIVE WVC = TWO*WVCG END IF ! END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE CG_FROM_TABLE(SIGIN,DEPTH,WVCG) USE DATAPOOL IMPLICIT NONE REAL(rkind), INTENT(IN) :: SIGIN, DEPTH REAL(rkind), INTENT(OUT) :: WVCG REAL(rkind) :: SQRTDEP, SIGSQDEP, R1, DEPLOC INTEGER :: ITAB DEPLOC = MAX(DMIN,DEPTH) SQRTDEP = SQRT(DEPLOC) SIGSQDEP = SIGIN*SQRTDEP ITAB = INT(SIGSQDEP/DSIGTAB) IF (ITAB.LE.NMAX.AND.ITAB.GE.0) THEN R1 = SIGSQDEP/DSIGTAB - MyREAL(ITAB) WVCG = ((1.-R1)*TABCG(ITAB)+R1*TABCG(ITAB+1))*SQRTDEP ELSE WVCG = .5*G9/SIGIN END IF END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE K_FROM_TABLE(SIGIN,DEPTH,WVK) USE DATAPOOL IMPLICIT NONE REAL(rkind), INTENT(IN) :: SIGIN, DEPTH REAL(rkind), INTENT(OUT) :: WVK REAL(rkind) :: SQRTDEP, SIGSQDEP, R1, DEPLOC INTEGER :: ITAB DEPLOC = MAX(DMIN,DEPTH) SQRTDEP = SQRT(DEPLOC) SIGSQDEP = SIGIN * SQRTDEP ITAB = INT(SIGSQDEP/DSIGTAB) IF (ITAB.LE.NMAX.AND.ITAB.GE.0) THEN R1 = SIGSQDEP/DSIGTAB - MyREAL(ITAB) WVK = ((1.-R1)*TABK(ITAB)+R1*TABK(ITAB + 1))/DEPLOC ELSE WVK = SIGIN*SIGIN/G9 END IF END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE WAVE_K_C_CG USE DATAPOOL IMPLICIT NONE INTEGER :: IP, IS REAL(rkind) :: DEPLOC REAL(rkind) :: WVK,WVCG,WVKDEP,WVN,WVC,SPSIGLOC DO IP = 1, MNP DEPLOC = MAX(DMIN,DEP(IP)) DO IS = 1, MSC SPSIGLOC = SPSIG(IS) CALL ALL_FROM_TABLE(SPSIGLOC,DEPLOC,WVK,WVCG,WVKDEP,WVN,WVC) ! CALL WAVEKCG(DEPLOC,SPSIGLOC,WVN,WVC,WVK,WVCG) WK(IS,IP) = WVK CG(IS,IP) = WVCG WC(IP,IS) = WVC END DO END DO END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE CHECK_STEADY(TIME, CONV1, CONV2, CONV3, CONV4, CONV5) USE DATAPOOL IMPLICIT NONE ! !AR: Joseph please check this code ... ! REAL(rkind), INTENT(IN) :: TIME REAL(rkind), INTENT(OUT) :: CONV1, CONV2, CONV3, CONV4, CONV5 INTEGER :: I, IP, IE, IS, ID, NI(3) INTEGER :: IPCONV1, IPCONV2, IPCONV3, IPCONV4, IPCONV5, ISCONV(MNP) REAL(rkind) :: SUMAC, ACLOC(MSC,MDC) REAL(rkind) :: ETOT, EAD, DS, HS2, KD REAL(rkind) :: ETOTF3, ETOTF4, TP, KHS2, EFTOT, TM02 REAL(rkind) :: FP, CP, KPP, CGP, WNP, UXD, OMEG, OMEG2 REAL(rkind) :: CONVK1, CONVK2, CONVK3, CONVK4, CONVK5 INTEGER :: ITMP IPCONV1 = 0 IPCONV2 = 0 IPCONV3 = 0 IPCONV4 = 0 IPCONV5 = 0 #ifndef MPI_PARALL_GRID DO IP = 1, MNP #else DO IP = 1, NP_RES IF(ASSOCIATED(IPGL(IPLG(IP))%NEXT)) THEN !interface node IF(IPGL(IPLG(ip))%NEXT%RANK < MYRANK) CYCLE !already in the sum so skip ENDIF #endif ACLOC(:,:) = AC2(:,:,IP) SUMAC = SUM(ACLOC) ETOT = 0.0 DO ID = 1, MDC DO IS = 2, MSC DS = SPSIG(IS) - SPSIG(IS-1) EAD = 0.5*(SPSIG(IS)*ACLOC(IS,ID)+SPSIG(IS-1)*ACLOC(IS-1,ID))*DS*DDIR ETOT = ETOT + EAD END DO END DO HS2 = 4.0_rkind*SQRT(ETOT) ETOTF3 = 0. ETOTF4 = 0. DO IS = 1, MSC DO ID = 1, MDC ETOTF3 = ETOTF3 + SPSIG(IS) * ACLOC(IS,ID)**4 * DDIR * DS_BAND(IS) ETOTF4 = ETOTF4 + ACLOC(IS,ID)**4 * DDIR * DS_BAND(IS) END DO END DO IF(ETOTF4 .GT. THR8 .AND. ETOTF3 .GT. THR8) THEN FP = ETOTF3/ETOTF4 ! CALL WAVEKCG(DEP(IP), FP, WNP, CP, KPP, CGP) CALL ALL_FROM_TABLE(FP,DEP(IP),KPP,CGP,KD,WNP,CP) TP = 1.0_rkind/FP/PI2 KHS2 = HS2 * KPP ELSE KHS2 = 0.0_rkind END IF ETOT = 0. EFTOT = 0. DO ID=1, MDC IF (LSECU .OR. LSTCU) THEN UXD = CURTXY(IP,1)*COSTH(ID) + CURTXY(IP,2)*SINTH(ID) ENDIF DO IS=1,MSC EAD = SIGPOW(IS,2) * ACLOC(IS,ID) * FRINTF IF (LSECU .OR. LSTCU) THEN OMEG = SPSIG(IS) + WK(IS,IP) * UXD OMEG2 = OMEG**2 ELSE OMEG2 = SIGPOW(IS,2) ENDIF ETOT = ETOT + EAD EFTOT = EFTOT + EAD * OMEG2 ENDDO ENDDO IF (ETOT/EFTOT .GT. THR) THEN TM02 = PI2 * SQRT(ETOT/EFTOT) ELSE TM02 = 0. END IF IF (DEP(IP) .LT. DMIN .OR. SUMAC .LT. THR .OR. IOBP(IP) .EQ. 2 .OR. HS2 .LT. THR) THEN IPCONV1 = IPCONV1 + 1 ! Summation of the converged grid points ... IPCONV2 = IPCONV2 + 1 IPCONV3 = IPCONV3 + 1 IPCONV4 = IPCONV4 + 1 IPCONV5 = IPCONV5 + 1 ISCONV(IP) = 1 !write(dbg%fhndl,*) 'boundary -------1------', TIME, IP, IP_IS_STEADY(IP) CYCLE ELSE CONVK1 = ABS(HSOLD(IP)-HS2)/HS2 CONVK2 = ABS(HS2-HSOLD(IP)) CONVK3 = ABS(SUMACOLD(IP)-SUMAC)/SUMAC CONVK4 = ABS(KHS2-KHSOLD(IP))/KHSOLD(IP) CONVK5 = ABS(TM02-TM02OLD(IP))/TM02OLD(IP) IF (CONVK1 .LT. EPSH1) IPCONV1 = IPCONV1 + 1 IF (CONVK2 .LT. EPSH2) IPCONV2 = IPCONV2 + 1 IF (CONVK3 .LT. EPSH3) IPCONV3 = IPCONV3 + 1 IF (CONVK4 .LT. EPSH4) IPCONV4 = IPCONV4 + 1 IF (CONVK5 .LT. EPSH5) IPCONV5 = IPCONV5 + 1 IF (CONVK1 .LT. EPSH1 .AND. CONVK2 .LT. EPSH2 .AND. CONVK3 .LT. EPSH3 .AND. CONVK4 .LT. EPSH4 .AND. CONVK5 .LT. EPSH5) THEN ISCONV(IP) = 1 !write(dbg%fhndl,*) 'converged -------2------', TIME, IP, IP_IS_STEADY(IP) ELSE ISCONV(IP) = 0 !write(dbg%fhndl,*) 'not converged -------3------', TIME, IP, IP_IS_STEADY(IP) ENDIF END IF HSOLD(IP) = HS2 SUMACOLD(IP) = SUMAC KHSOLD(IP) = KHS2 TM02OLD(IP) = TM02 END DO ! IP !write(*,*) time, maxval(IP_IS_STEADY), minval(IP_IS_STEADY) DO IE = 1, MNE NI = INE(:,IE) IF (SUM(ISCONV(NI)) .EQ. 3) THEN IE_IS_STEADY(IE) = IE_IS_STEADY(IE) + 1 ELSE IE_IS_STEADY(IE) = 0 ENDIF !WRITE(*,*) IE, IE_IS_STEADY(IE) ENDDO DO IP = 1, MNP ITMP = 0 DO I = 1, CCON(IP) IF (IE_IS_STEADY(IE_CELL2(IP,I)) .GE. 1) ITMP = ITMP + 1 ENDDO IF (ITMP .EQ. CCON(IP)) THEN IP_IS_STEADY(IP) = IP_IS_STEADY(IP) + 1 !IF (IP_IS_STEADY(IP) .GT. 2) WRITE(*,*) TIME, IP, IP_IS_STEADY(IP) ELSE IP_IS_STEADY(IP) = 0 ENDIF ENDDO #ifdef MPI_PARALL_GRID CALL MPI_ALLREDUCE(IPCONV1, itmp, 1, itype, MPI_SUM, COMM, ierr) IPCONV1 = itmp CALL MPI_ALLREDUCE(IPCONV2, itmp, 1, itype, MPI_SUM, COMM, ierr) IPCONV2 = itmp CALL MPI_ALLREDUCE(IPCONV3, itmp, 1, itype, MPI_SUM, COMM, ierr) IPCONV3 = itmp CALL MPI_ALLREDUCE(IPCONV4, itmp, 1, itype, MPI_SUM, COMM, ierr) IPCONV4 = itmp CALL MPI_ALLREDUCE(IPCONV5, itmp, 1, itype, MPI_SUM, COMM, ierr) IPCONV5 = itmp CONV1 = MyREAL(IPCONV1)/MyREAL(NP_GLOBAL)*100.0 CONV2 = MyREAL(IPCONV2)/MyREAL(NP_GLOBAL)*100.0 CONV3 = MyREAL(IPCONV3)/MyREAL(NP_GLOBAL)*100.0 CONV4 = MyREAL(IPCONV4)/MyREAL(NP_GLOBAL)*100.0 CONV5 = MyREAL(IPCONV5)/MyREAL(NP_GLOBAL)*100.0 #else CONV1 = MyREAL(IPCONV1)/MyREAL(MNP)*100.0 CONV2 = MyREAL(IPCONV2)/MyREAL(MNP)*100.0 CONV3 = MyREAL(IPCONV3)/MyREAL(MNP)*100.0 CONV4 = MyREAL(IPCONV4)/MyREAL(MNP)*100.0 CONV5 = MyREAL(IPCONV5)/MyREAL(MNP)*100.0 #endif #ifdef MPI_PARALL_GRID IF (myrank == 0) THEN WRITE(STAT%FHNDL,*) 'CONVERGENCE CRIT. 1 REACHED IN', CONV1, '% GRIDPOINTS' WRITE(STAT%FHNDL,*) 'CONVERGENCE CRIT. 2 REACHED IN', CONV2, '% GRIDPOINTS' WRITE(STAT%FHNDL,*) 'CONVERGENCE CRIT. 3 REACHED IN', CONV3, '% GRIDPOINTS' WRITE(STAT%FHNDL,*) 'CONVERGENCE CRIT. 4 REACHED IN', CONV4, '% GRIDPOINTS' WRITE(STAT%FHNDL,*) 'CONVERGENCE CRIT. 5 REACHED IN', CONV5, '% GRIDPOINTS' END IF #else WRITE(STAT%FHNDL,*) 'CONVERGENCE CRIT. 1 REACHED IN', CONV1, '% GRIDPOINTS' WRITE(STAT%FHNDL,*) 'CONVERGENCE CRIT. 2 REACHED IN', CONV2, '% GRIDPOINTS' WRITE(STAT%FHNDL,*) 'CONVERGENCE CRIT. 3 REACHED IN', CONV3, '% GRIDPOINTS' WRITE(STAT%FHNDL,*) 'CONVERGENCE CRIT. 4 REACHED IN', CONV4, '% GRIDPOINTS' WRITE(STAT%FHNDL,*) 'CONVERGENCE CRIT. 5 REACHED IN', CONV5, '% GRIDPOINTS' #endif FLUSH(STAT%FHNDL) END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE TWOD2ONED(ACLOC,AC1D) USE DATAPOOL IMPLICIT NONE REAL(rkind), INTENT(IN) :: ACLOC(MSC,MDC) REAL(rkind), INTENT(OUT) :: AC1D(MSC*MDC) INTEGER :: IS, ID DO IS = 1, MSC DO ID = 1, MDC AC1D(ID + (IS-1) * MDC) = ACLOC(IS,ID) END DO END DO END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE ONED2TWOD(AC1D,ACLOC) USE DATAPOOL IMPLICIT NONE REAL(rkind), INTENT(OUT) :: ACLOC(MSC,MDC) REAL(rkind), INTENT(IN) :: AC1D(MSC*MDC) INTEGER :: IS, ID DO IS = 1, MSC DO ID = 1, MDC ACLOC(IS,ID) = AC1D(ID + (IS-1) * MDC) END DO END DO END SUBROUTINE !********************************************************************** !* * !********************************************************************** ! REAL(rkind) FUNCTION GAMMA_FUNC(XX) FUNCTION GAMMA_FUNC(XX) USE DATAPOOL, ONLY: rkind IMPLICIT NONE REAL(rkind) :: GAMMA_FUNC ! ! Purpose: ! Compute the transcendental function Gamma ! ! Subroutines used ! GAMMLN (Numerical Recipes) ! real(rkind) GAMMLN REAL(rkind) XX, YY, ABIG SAVE ABIG DATA ABIG /30./ YY = GAMMLN(XX) IF (YY > ABIG) YY = ABIG IF (YY < -ABIG) YY = -ABIG GAMMA_FUNC = EXP(YY) END FUNCTION !********************************************************************** !* * !********************************************************************** FUNCTION GAMMLN(XX) USE DATAPOOL, ONLY: rkind, ONEHALF, ONE IMPLICIT NONE real(rkind) :: GAMMLN ! ! Method: ! function is copied from: Press et al., "Numerical Recipes" ! real(rkind) XX INTEGER J real(rkind) COF(6),STP,FPF,X,TMP,SER DATA COF,STP/76.18009173_rkind,-86.50532033_rkind, & & 24.01409822_rkind,-1.231739516_rkind, & & .120858003E-2_rkind,-.536382E-5_rkind, & & 2.50662827465_rkind/ DATA FPF/5.5_rkind/ X = XX-ONE TMP = X+FPF TMP = (X+ONEHALF)*LOG(TMP)-TMP SER = ONE DO J = 1, 6 X = X+ONE SER = SER+COF(J)/X END DO GAMMLN = TMP+LOG(STP*SER) END FUNCTION !********************************************************************** !* * !********************************************************************** FUNCTION VEC2RAD(U,V) USE DATAPOOL, ONLY: rkind, pi IMPLICIT NONE REAL(RKIND) :: VEC2RAD REAL(rkind) :: U,V VEC2RAD = MyATAN2(V,U) * 180/PI IF (VEC2RAD < 0.0) VEC2RAD = VEC2RAD + 360.0_rkind VEC2RAD = VEC2RAD * PI/180. END FUNCTION !********************************************************************** !* * !********************************************************************** FUNCTION VEC2DEG(U,V) USE DATAPOOL, ONLY: rkind, pi IMPLICIT NONE REAL(RKIND) :: VEC2DEG REAL(rkind), intent(in) :: U,V VEC2DEG = MyATAN2(V,U) * 180./PI IF (VEC2DEG < 0.0) VEC2DEG = VEC2DEG + 360.0_rkind END FUNCTION !********************************************************************** !* * !********************************************************************** FUNCTION DVEC2RAD(U,V) USE DATAPOOL, ONLY: rkind, pi IMPLICIT NONE REAL(rkind) :: U,V REAL(rkind) :: DVEC2RAD DVEC2RAD = MyATAN2(V,U) * 180.0_rkind/PI IF (DVEC2RAD < 0.0_rkind) DVEC2RAD = DVEC2RAD + 360.0_rkind DVEC2RAD = DVEC2RAD * PI/180.0_rkind END FUNCTION !********************************************************************** !* * !********************************************************************** FUNCTION DVEC2DEG(U,V) USE DATAPOOL, ONLY : PI, rkind IMPLICIT NONE REAL(rkind) :: DVEC2DEG REAL(rkind) :: U,V DVEC2DEG = MyATAN2(V,U) * 180.0_rkind/PI IF (DVEC2DEG < 0.0_rkind) DVEC2DEG = DVEC2DEG + 360.0_rkind END FUNCTION !********************************************************************** !* * !********************************************************************** SUBROUTINE DEG2NAUT (DEGREE, DEG, LTRANS) USE DATAPOOL, ONLY: rkind IMPLICIT NONE ! ! Nautical convention Cartesian convention ! (Where the wind/waves come from) (Where the wind/waves go to) ! 0 90 ! | | ! | | ! | | ! | | ! 270 --------+-------- 90 180 --------+-------- 0 ! | | ! | | ! | | ! | | ! 180 270 ! LOGICAL :: LTRANS REAL(rkind), INTENT(IN) :: DEGREE REAL(rkind), INTENT(OUT) :: DEG REAL(rkind) :: DNORTH IF ( LTRANS ) THEN DNORTH = 90.0 DEG = 180. + DNORTH - DEGREE ELSE DEG = DEGREE END IF ! IF (DEG .GE. 360.0_rkind) THEN DEG = MOD (DEG, 360.0_rkind) ELSE IF (DEG .LT. 0.) THEN DEG = MOD (DEG, 360.0_rkind) + 360.0_rkind ELSE ! DEG between 0 and 360; do nothing endif ! END SUBROUTINE !********************************************************************** !* * !********************************************************************** logical function my_isnan(a) use datapool, only : rkind real(rkind) :: a if (a.ne.a) then my_isnan = .true. else my_isnan = .false. end if return end !********************************************************************** !* * !********************************************************************** logical function my_isinf(a) use datapool, only : rkind real(rkind) :: a if (int(a*0) .ne. 0) then my_isinf = .true. else my_isinf = .false. end if return end !********************************************************************** !* * !********************************************************************** logical function iseq0(a) use datapool, only : rkind real(rkind) :: a if (abs(a) .lt. tiny(1.0)) then iseq0 = .true. else iseq0 = .false. end if return end !********************************************************************** !* * !********************************************************************** SUBROUTINE CHKVOL(VAL,PO,NE,POSNEG) USE DATAPOOL IMPLICIT NONE REAL(rkind), INTENT(IN) :: VAL(MNP) REAL(rkind), INTENT(OUT) :: PO,NE,POSNEG REAL(rkind) :: TMP INTEGER :: IE,I1,I2,I3 PO = ZERO NE = ZERO DO IE = 1, MNE I1 = INE(1,IE) I2 = INE(2,IE) I3 = INE(3,IE) TMP = ONETHIRD * (VAL(I1)+VAL(I2)+VAL(I3)) * TRIA(IE) IF (TMP .GT. ZERO) THEN PO = PO + TMP ELSE NE = NE + TMP END IF POSNEG = PO + NE END DO END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE CHECKCONS(VAL,SUMAC) USE DATAPOOL IMPLICIT NONE REAL(rkind), INTENT(IN) :: VAL(MNP) REAL(rkind), INTENT(OUT) :: SUMAC REAL(rkind) :: TMP INTEGER :: IE,I1,I2,I3 SUMAC = 0. DO IE = 1, MNE I1 = INE(1,IE) I2 = INE(2,IE) I3 = INE(3,IE) TMP = 1./3. * (VAL(I1)+VAL(I2)+VAL(I3)) * TRIA(IE) SUMAC = SUMAC + TMP END DO END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE ADVTEST(INIT) USE DATAPOOL IMPLICIT NONE INTEGER :: IP REAL(rkind) :: R REAL(rkind), INTENT(OUT) :: INIT(MNP) INIT = 0.0_rkind DO IP = 1, MNP R = SQRT( (XP(IP)+ONEHALF)**2 + YP(IP)**2 ) IF ( R <= 0.25_rkind ) THEN IF (LZYLINDER) THEN INIT(IP) = ONE ELSE INIT(IP) = COS(PI2*R) END IF ELSE INIT(IP) = ZERO END IF END DO WRITE(4001) SNGL(RTIME) WRITE(4001) (1., 1., SNGL(INIT(IP)), IP = 1, MNP) END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE ERG2WWM(STEPS) USE DATAPOOL IMPLICIT NONE INTEGER :: T, IP, STEPS REAL(rkind) :: HEADSP REAL(rkind), PARAMETER :: FRFAK = 0.9 REAL(rkind) :: VEC2RAD, ANG, VEL, WAVEL, DEPTH REAL(rkind), ALLOCATABLE :: HP(:), QU(:), QV(:), UP(:), VP(:) HEADSP = 0.0 #ifdef MPI_PARALL_GRID CALL WWM_ABORT('ERG2WWM CANNOT BE CALLED FROM SCHISM') #endif ALLOCATE (QU(MNP), QV(MNP), UP(MNP), VP(MNP), HP(MNP), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_aux, allocate error 1') QU = 0. QV = 0. UP = 0. VP = 0. HP = 0. OPEN(1230, FILE = 'quqvh.bin', FORM = 'UNFORMATTED') OPEN(1240, FILE = 'current.dat') OPEN(1250, FILE = 'wlevel.dat') DO T = 1, STEPS READ(1230) HEADSP READ(1230) (QU(IP), QV(IP), HP(IP) , IP = 1, MNP) WRITE (1240, *) HEADSP WRITE (1250, *) HEADSP DO IP = 1, MNP DEPTH = DEP(IP)+HP(IP) IF ( DEPTH .GT. DMIN ) THEN UP(IP) = QU(IP)/DEPTH VP(IP) = QV(IP)/DEPTH VEL = SQRT(UP(IP)**2 + VP(IP)**2) ANG = VEC2RAD(UP(IP),VP(IP)) WAVEL = SQRT(G9*DEPTH) IF (VEL .GT. FRFAK * WAVEL) THEN UP(IP) = COS(ANG)*FRFAK*WAVEL VP(IP) = SIN(ANG)*FRFAK*WAVEL END IF ELSE UP(IP) = 0.0 VP(IP) = 0.0 END IF END DO WRITE (1240, 1200) UP(:) WRITE (1240, 1200) VP(:) WRITE (1250, 1200) HP(:) END DO DEALLOCATE (QU,QV,HP,UP,VP) 1200 FORMAT(8F9.4) CLOSE (1240) CLOSE (1250) END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE INTELEMENT_COEF(X,Y,XP,YP,WI) USE DATAPOOL, ONLY : RKIND IMPLICIT NONE REAL(rkind), INTENT(IN) :: X(3), Y(3) REAL(rkind), INTENT(IN) :: XP, YP REAL(rkind), INTENT(OUT) :: WI(3) REAL(rkind) :: y1,y2,y3,x1,x2,x3 REAL(rkind) :: n1, d1, n2, d2, n3, d3 x1 = X(1) x2 = X(2) x3 = X(3) y1 = Y(1) y2 = Y(2) y3 = Y(3) n1=(XP-x2)*(y3-y2) - (YP-y2)*(x3-x2) d1=(x1-x2)*(y3-y2) - (y1-y2)*(x3-x2) n2=(XP-x1)*(y3-y1) - (YP-y1)*(x3-x1) d2=(x2-x1)*(y3-y1) - (y2-y1)*(x3-x1) n3=(XP-x1)*(y2-y1) - (YP-y1)*(x2-x1) d3=(x3-x1)*(y2-y1) - (y3-y1)*(x2-x1) Wi(1)=n1/d1 Wi(2)=n2/d2 Wi(3)=n3/d3 END SUBROUTINE INTELEMENT_COEF !********************************************************************** !* * !********************************************************************** SUBROUTINE INTELEMENT(X,Y,Z,XP,YP,Wi,Zi,LSAME) USE DATAPOOL, ONLY : RKIND, THR8 IMPLICIT NONE LOGICAL, INTENT(IN) :: LSAME REAL(rkind), INTENT(IN) :: X(3), Y(3), Z(3) REAL(rkind), INTENT(IN) :: XP, YP REAL(rkind), INTENT(OUT) :: Zi REAL(rkind), INTENT(OUT) :: WI(3) REAL(rkind) :: y1,y2,y3,x1,x2,x3,z1,z2,z3 REAL(rkind), SAVE :: A,B,C,D IF (.NOT. LSAME) THEN x1 = X(1) x2 = X(2) x3 = X(3) y1 = Y(1) y2 = Y(2) y3 = Y(3) z1 = Z(1) z2 = Z(2) z3 = Z(3) A = y1*(z2 - z3) + y2*(z3 - z1) + y3*(z1 - z2) B = z1*(x2 - x3) + z2*(x3 - x1) + z3*(x1 - x2) C = x1*(y2 - y3) + x2*(y3 - y1) + x3*(y1 - y2) D = -A*x1 - B*y1 - C*z1 IF (ABS(C) .GT. THR8 ) THEN WI(1) = -A/C WI(2) = -B/C WI(3) = -D/C ELSE WI = 0.0_rkind endif END IF Zi = WI(1) * XP + WI(2) * YP + WI(3) END SUBROUTINE INTELEMENT !********************************************************************** !* * !********************************************************************** SUBROUTINE INTELEMENT_AC_LOC(I,ACLOC,CURTXYLOC,DEPLOC,WATLEVLOC,WKLOC) USE DATAPOOL IMPLICIT NONE INTEGER, INTENT(IN) :: I REAL(rkind), INTENT(OUT) :: ACLOC(MSC,MDC) REAL(rkind), INTENT(OUT) :: CURTXYLOC(2), DEPLOC, WATLEVLOC, WKLOC(MSC) REAL(rkind), SAVE :: WI(3) INTEGER :: IS, NI(3), IE REAL(rkind) :: WVN, WVC, WVK, WVCG, WVKDEP integer IP, J IE = STATION(I)%ELEMENT WI = STATION(I)%WI NI = INE(:,IE) DEPLOC = ZERO WATLEVLOC = ZERO CURTXYLOC = ZERO ACLOC = ZERO DO J=1,3 IP = INE(J,IE) DEPLOC = DEPLOC + WI(J) * DEP(IP) CURTXYLOC(:) = CURTXYLOC(:) + WI(J) * CURTXY(IP,:) WATLEVLOC = WATLEVLOC + WI(J) * WATLEV(IP) ACLOC(:,:) = ACLOC(:,:) + WI(J) * AC2(:,:,IP) !write(*,'(3I10,5F15.8)') j, ie, ip, wi(j), DEP(IP), WATLEV(IP), CURTXY(IP,:) END DO DO IS = 1, MSC !CALL WAVEKCG(DEPLOC, SPSIG(IS), WVN, WVC, WVK, WVCG) CALL ALL_FROM_TABLE(SPSIG(IS),DEPLOC,WVK,WVCG,WVKDEP,WVN,WVC) WKLOC(IS) = WVK END DO END SUBROUTINE INTELEMENT_AC_LOC !********************************************************************** !* * !********************************************************************** SUBROUTINE INTELEMENT_WW3GLOBAL_LOC(IE,XPC,YPC,WW3LOCAL) USE DATAPOOL IMPLICIT NONE INTEGER, INTENT(IN) :: IE REAL(rkind), INTENT(IN) :: XPC, YPC REAL(rkind), INTENT(OUT) :: WW3LOCAL(8) REAL(rkind) :: XOUTELE(3), YOUTELE(3) REAL(rkind), SAVE :: WI(3) INTEGER :: I, NI(3) LOGICAL :: LSAME NI = INE(:,IE) XOUTELE = XP(NI) YOUTELE = YP(NI) LSAME = .FALSE. DO I = 1, 8 CALL INTELEMENT(XOUTELE,YOUTELE,WW3GLOBAL(I,NI),XPC,YPC,WI,WW3LOCAL(I),LSAME) LSAME = .TRUE. END DO END SUBROUTINE INTELEMENT_WW3GLOBAL_LOC !********************************************************************** !* * !********************************************************************** SUBROUTINE CSEVAL ( NFU, FILEN, LSE, DIMS, SEVAL, MULTIPLE_IN) USE DATAPOOL IMPLICIT NONE INTEGER, INTENT(IN) :: NFU LOGICAL, INTENT(IN) :: MULTIPLE_IN CHARACTER(LEN=*), INTENT(IN) :: FILEN LOGICAL, INTENT(IN) :: LSE INTEGER, INTENT(IN) :: DIMS REAL(rkind), INTENT(INOUT) :: SEVAL(MNP, DIMS) REAL(rkind) :: SEVAL2(NP_TOTAL, DIMS) #ifdef MPI_PARALL_GRID INTEGER :: IP REAL(rkind) :: Vtotal(np_total), Vlocal(MNP) #endif INTEGER :: IC, IFSTAT, IDIM CHARACTER(LEN=128) :: HEADLN #ifdef MPI_PARALL_GRID IF (MULTIPLE_IN .or. (myrank .eq. 0)) THEN #endif IF (LSE) THEN READ(NFU,*) HEADLN WRITE(STAT%FHNDL,'("+TRACE...",2A)') 'Reading the header of the serial file ... HEADER ', TRIM(HEADLN) DO IC = 1, DIMS READ( NFU, *, IOSTAT = IFSTAT ) SEVAL2(:, IC) IF ( IFSTAT /= 0 ) CALL WWM_ABORT('unexpected error reading the serial file in CSEVAL 1') END DO ELSE OPEN( NFU, FILE = TRIM(FILEN), STATUS = 'OLD') WRITE(STAT%FHNDL,'("+TRACE...",2A)') 'Reading the file of the request ... ', TRIM(FILEN) READ(NFU,*) HEADLN DO IC = 1, DIMS READ( NFU, *, IOSTAT = IFSTAT ) SEVAL2(:, IC) IF ( IFSTAT /= 0 ) CALL WWM_ABORT(' unexpected error reading the serial file in CSEVAL 2') END DO CLOSE( NFU ) END IF #ifdef MPI_PARALL_GRID END IF #endif #ifdef MPI_PARALL_GRID IF (MULTIPLE_IN) THEN DO IP=1,MNP SEVAL(IP,:) = SEVAL2(IPLG(IP),:) END DO ELSE DO IDIM=1,DIMS IF (myrank .eq. 0) THEN Vtotal=SEVAL2(:,IDIM) END IF CALL SCATTER_ONED_ARRAY(Vtotal, Vlocal) SEVAL(:,IDIM)=Vlocal END DO END IF #else SEVAL(:,:) = SEVAL2(:,:) #endif END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE ERROR(X,ERR) USE DATAPOOL, ONLY : rkind,pi IMPLICIT NONE REAL(rkind) :: EPS,X,X2,ERR,ER,R,C0 INTEGER :: K ! ========================================= ! Purpose: Compute error function erf(x) ! Input: x --- Argument of erf(x) ! Output: ERR --- erf(x) ! ========================================= EPS=1.E-15_rkind X2=X*X IF (MyABS(X).LT.3.5D0) THEN ER=1.0D0 R=1.0D0 DO 10 K=1,50 R=R*X2/(K+0.5D0) ER=ER+R IF (MyABS(R).LE.MyABS(ER)*EPS) GO TO 15 10 CONTINUE 15 C0=2.0D0/MySQRT(PI)*X*MyEXP(-X2) ERR=C0*ER ELSE ER=1.0D0 R=1.0D0 DO K=1,12 R=-R*(K-0.5D0)/X2 ER=ER+R END DO C0=MyEXP(-X2)/(MyABS(X)*MySQRT(PI)) ERR=1.0D0-C0*ER IF (X.LT.0.0) ERR=-ERR endif END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE WRINPGRD_XFN() #ifdef MPI_PARALL_GRID USE DATAPOOL, ONLY : myrank, comm, ierr, MNP, XP, YP, DEP, MNE, INE #else USE DATAPOOL, ONLY : MNP, XP, YP, DEP, MNE, INE #endif IMPLICIT NONE INTEGER :: I #ifdef MPI_PARALL_GRID if (myrank == 0) then #endif OPEN(2222, FILE='systest.dat', STATUS='UNKNOWN') CALL XFNHEADER_1(2222,0,MNP) DO I = 1, MNP WRITE(2222,'(I10,2F20.8,F15.4)') I-1, XP(I), YP(I), DEP(I) END DO CALL XFNHEADER_2(2222,MNE) DO I = 1, MNE WRITE(2222,'(5I10)') INE(1,I)-1, INE(2,I)-1, INE(3,I)-1, 0, I-1 END DO #ifdef MPI_PARALL_GRID endif #endif END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE WRINPGRD_SHP() USE DATAPOOL IMPLICIT NONE INTEGER :: I #ifdef MPI_PARALL_GRID if (myrank == 0) then #endif OPEN(2222, FILE='systest.dat', STATUS='UNKNOWN') CALL XFNHEADER_1(2222,0,MNP) DO I = 1, MNP WRITE(2222,'(I10,2F20.8,F15.4)') I-1, XP(I), YP(I), DEP(I) END DO CALL XFNHEADER_2(2222,MNE) DO I = 1, MNE WRITE(2222,'(5I10)') INE(1,I)-1, INE(2,I)-1, INE(3,I)-1, 0, I-1 END DO #ifdef MPI_PARALL_GRID endif #endif END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE XFNHEADER_1(IFILE,NKR,NKG) IMPLICIT NONE INTEGER, INTENT(IN) :: IFILE, NKR, NKG WRITE(IFILE,'(A)') 'C system.dat, made by tri2sys' WRITE(IFILE,'(A)') 'C Number of Boundary Nodes:' WRITE(IFILE,'(I10)') NKR WRITE(IFILE,'(A)') 'C Number of Domain Nodes:' WRITE(IFILE,'(I10)') NKG WRITE(IFILE,'(A)') 'C Koordinaten und Skalarwerte der Knoten' WRITE(IFILE,'(A)') 'C --------------------------------------' WRITE(IFILE,'(A)') 'C Zuerst die Randknoten (Anzahl s.o.),' WRITE(IFILE,'(A)') 'C dann die Gebietsknoten (Anzahl s.o.).' WRITE(IFILE,'(A)') 'C ------------+-------------+-------------+---------------' WRITE(IFILE,'(A)') 'C Nr. | x-Koord. | y-Koord. | Skalarwert' WRITE(IFILE,'(A)') 'C ------------+-------------+-------------+---------------' END SUBROUTINE XFNHEADER_1 !********************************************************************** !* * !********************************************************************** SUBROUTINE XFNHEADER_2(IFILE,NELEM) IMPLICIT NONE INTEGER, INTENT(IN) :: IFILE, NELEM WRITE(IFILE,'(A)') "C ------------------------------------------------------------" WRITE(IFILE,'(A)') "C Anzahl der Elemente:" WRITE(IFILE,'(I11)') NELEM WRITE(IFILE,'(A)') "C Elementverzeichnis" WRITE(IFILE,'(A)') "C ------------------------------------------------------------" WRITE(IFILE,'(A)') "C Knoten i Knoten j Knoten k Kennung Nr." END SUBROUTINE XFNHEADER_2 !********************************************************************** !* * !********************************************************************** SUBROUTINE RHEADER_NODE(IFILE,NKR,NKG) IMPLICIT NONE INTEGER, INTENT(IN) :: IFILE INTEGER, INTENT(OUT):: NKR, NKG READ(IFILE,*) READ(IFILE,*) READ(IFILE,*) NKR READ(IFILE,*) READ(IFILE,*) NKG READ(IFILE,*) READ(IFILE,*) READ(IFILE,*) READ(IFILE,*) READ(IFILE,*) READ(IFILE,*) READ(IFILE,*) END SUBROUTINE RHEADER_NODE !********************************************************************** !* * !********************************************************************** SUBROUTINE RHEADER_ELEMENT(IFILE,NELEM) IMPLICIT NONE INTEGER, INTENT(IN) :: IFILE INTEGER, INTENT(OUT):: NELEM READ(IFILE,*) READ(IFILE,*) READ(IFILE,*) NELEM READ(IFILE,*) READ(IFILE,*) READ(IFILE,*) END SUBROUTINE RHEADER_ELEMENT !********************************************************************** !* * !********************************************************************** SUBROUTINE FIND_ELE ( M,N,Lelem,Xkno,Xp,Yp,Ele ) USE DATAPOOL, ONLY : THR8, RKIND IMPLICIT NONE ! ! Dummy arguments ! INTEGER, INTENT(IN) :: M, N REAL(rkind), INTENT(IN) :: Xkno(2,N), Xp , Yp INTEGER, INTENT(IN) :: Lelem(3,M) INTEGER, INTENT(INOUT) :: Ele ! ! Local variables ! REAL(rkind), SAVE :: xi, xj, xk, yi, yj, yk, dx, dy, f, Xp8, Yp8 REAL(rkind), SAVE :: xmax , xmin , ymax , ymin INTEGER, SAVE :: i , i0 , idx , ielem , if0 , if1 , ijk , k , ki , kj , kk , l DATA idx/0/ , i0/0/ , ielem/1/ , if0/0/ , if1/0/ ! ! Laengste Kannte (DX,DY) bestimmen ! Dieser Programmabschnitt wird nur beim ersten Aufruf von PLO160 ! durchlaufen! ! Xp8 = Xp Yp8 = Yp IF ( idx/=1 ) THEN idx = 1 DO i = 1 , M ki = Lelem(1,i)! + 1 kj = Lelem(2,i)! + 1 kk = Lelem(3,i)! + 1 xi = Xkno(1,ki) yi = Xkno(2,ki) xj = Xkno(1,kj) yj = Xkno(2,kj) xk = Xkno(1,kk) yk = Xkno(2,kk) IF ( i==1 ) THEN dx = MAX(ABS(xi-xj),ABS(xi-xk),ABS(xj-xk)) dy = MAX(ABS(yi-yj),ABS(yi-yk),ABS(yj-yk)) ELSE dx = MAX(dx,ABS(xi-xj),ABS(xi-xk),ABS(xj-xk)) dy = MAX(dy,ABS(yi-yj),ABS(yi-yk),ABS(yj-yk)) endif ENDDO endif ! ------------------------------------------------------------------ ! TEST, OB DER PUNKT IM ZULETZT ANGESPROCHENEN ELEMENT LIEGT ! ------------------------------------------------------------------ IF ( i0==1 .AND. Ele/=0 ) THEN IF ( Yp8-ymin > THR8 ) THEN IF ( Yp8-ymax < -THR8 ) THEN IF ( Xp8-xmin > THR8 ) THEN IF ( Xp8-xmax < -THR8 ) THEN f = xi*(yj-Yp8) + xj*(Yp8-yi) + Xp8*(yi-yj) IF ( f > THR8 ) THEN f = xj*(yk-Yp8) + xk*(Yp8-yj) + Xp8*(yj-yk) IF ( f > THR8 ) THEN f = xk*(yi-Yp8) + xi*(Yp8-yk) + Xp8*(yk-yi) IF ( f > THR8 ) THEN Ele = ielem ! Element gefunden -->RETURN RETURN endif endif endif endif endif endif endif endif ! ------------------------------------------------------------------ ! Element suchen ! ------------------------------------------------------------------ i0 = 1 i = ielem IF ( i<1 ) i = 1 k = i l = i ijk = 0 100 DO ijk = ijk + 1 !..... ABFRAGE AUF X-Richtung ki = Lelem(1,i)! + 1 xi = Xkno(1,ki) IF ( MyABS(xi-Xp8)<=dx ) THEN kj = Lelem(2,i)! + 1 kk = Lelem(3,i)! + 1 xj = Xkno(1,kj) xk = Xkno(1,kk) !..... Punkt ausserhalb Element: xmin = MIN(xi,xj,xk) IF ( Xp8>=xmin ) THEN xmax = MAX(xi,xj,xk) IF ( Xp8<=xmax ) THEN !..... ABFRAGE AUF Y-Richtung yi = Xkno(2,ki) IF ( MyABS(yi-Yp8)<=dy ) THEN yj = Xkno(2,kj) yk = Xkno(2,kk) !..... Punkt ausserhalb Element: ymin = MIN(yi,yj,yk) IF ( Yp8>=ymin ) THEN ymax = MAX(yi,yj,yk) IF ( Yp8<=ymax ) THEN !..... Bis jetzt liegt Punkt innerhalb des das Element ! umschlieszenden Rechtecks XMIN/XMAX, YMIN/YMAX ! Pruefen, ob Punkt wirklich innerhalb DREIECK-Element ! liegt: BERECHNUNG DER TEILFLAECHEN (ohne 0.5) f = xi*(yj-Yp8) + xj*(Yp8-yi) + Xp8*(yi-yj) IF ( f>=0.0_rkind) THEN f = xj*(yk-Yp8) + xk*(Yp8-yj) + Xp8*(yj-yk) IF ( f>=0.0_rkind ) THEN f = xk*(yi-Yp8) + xi*(Yp8-yk) + Xp8*(yk-yi) IF ( f>=0.0_rkind ) THEN Ele = i ielem = Ele RETURN endif endif endif endif endif endif endif endif endif ! SCHLEIFE UEBER ALLE ELEMENTE wird hier folgendermassen hochgezaehlt: ! beginnend bei IEALT, im Wechsel nach vorn und rueckwaerts suchend IF ( k1 ) if1 = 1 k = k + 1 i = k IF ( ijk<=M ) CYCLE endif CONTINUE EXIT ENDDO IF ( l>1 .AND. if0==0 ) THEN if1 = 0 IF ( k THR8 ) THEN IF ( Yp8-ymax < -THR8 ) THEN IF ( Xp8-xmin > THR8 ) THEN IF ( Xp8-xmax < -THR8 ) THEN f = xi*(yj-Yp8) + xj*(Yp8-yi) + Xp8*(yi-yj) IF ( f > THR8 ) THEN f = xj*(yk-Yp8) + xk*(Yp8-yj) + Xp8*(yj-yk) IF ( f > THR8 ) THEN f = xk*(yi-Yp8) + xi*(Yp8-yk) + Xp8*(yk-yi) IF ( f > THR8 ) THEN Ele = ielem ! Element gefunden -->RETURN RETURN endif endif endif endif endif endif endif endif ! ------------------------------------------------------------------ ! Element suchen ! ------------------------------------------------------------------ i0 = 1 i = ielem IF ( i<1 ) i = 1 k = i l = i ijk = 0 100 DO ijk = ijk + 1 !..... ABFRAGE AUF X-Richtung ki = Lelem(1,i)! + 1 xi = Xkno(1,ki) IF ( MyABS(xi-Xp8)<=dx ) THEN kj = Lelem(2,i)! + 1 kk = Lelem(3,i)! + 1 xj = Xkno(1,kj) xk = Xkno(1,kk) !..... Punkt ausserhalb Element: xmin = MIN(xi,xj,xk) IF ( Xp8>=xmin ) THEN xmax = MAX(xi,xj,xk) IF ( Xp8<=xmax ) THEN !..... ABFRAGE AUF Y-Richtung yi = Xkno(2,ki) IF ( MyABS(yi-Yp8)<=dy ) THEN yj = Xkno(2,kj) yk = Xkno(2,kk) !..... Punkt ausserhalb Element: ymin = MIN(yi,yj,yk) IF ( Yp8>=ymin ) THEN ymax = MAX(yi,yj,yk) IF ( Yp8<=ymax ) THEN !..... Bis jetzt liegt Punkt innerhalb des das Element ! umschlieszenden Rechtecks XMIN/XMAX, YMIN/YMAX ! Pruefen, ob Punkt wirklich innerhalb DREIECK-Element ! liegt: BERECHNUNG DER TEILFLAECHEN (ohne 0.5) f = xi*(yj-Yp8) + xj*(Yp8-yi) + Xp8*(yi-yj) IF ( f>=0.0_rkind) THEN f = xj*(yk-Yp8) + xk*(Yp8-yj) + Xp8*(yj-yk) IF ( f>=0.0_rkind ) THEN f = xk*(yi-Yp8) + xi*(Yp8-yk) + Xp8*(yk-yi) IF ( f>=0.0_rkind ) THEN Ele = i ielem = Ele RETURN endif endif endif endif endif endif endif endif endif ! SCHLEIFE UEBER ALLE ELEMENTE wird hier folgendermassen hochgezaehlt: ! beginnend bei IEALT, im Wechsel nach vorn und rueckwaerts suchend IF ( k1 ) if1 = 1 k = k + 1 i = k IF ( ijk<=M ) CYCLE endif CONTINUE EXIT ENDDO IF ( l>1 .AND. if0==0 ) THEN if1 = 0 IF ( k ZERO .AND. Y1 < 90._rkind) .AND. (Y2 < 360._rkind .AND. Y2 > 270._rkind)) THEN YINTER=(Y1+360._rkind)+(Y2-(Y1+360._rkind))/DX*DIFFDX IF (YINTER > 360._rkind) YINTER = YINTER - 360._rkind ELSE IF ((Y2 > ZERO .AND. Y2 < 90._rkind) .AND. (Y1 < 360._rkind .AND. Y1 > 270._rkind)) THEN YINTER = Y1+((Y2+360._rkind)-Y1)/DX*DIFFDX IF (YINTER > 360._rkind) YINTER = YINTER - 360._rkind ELSE YINTER = Y1+(Y2-Y1)/DX*DIFFDX END IF IF (ABS(YINTER) .LT. TINY(1.)) YINTER = Y1 END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE INTERLIN (NX1, NX2, X1, X2, Y1, Y2) USE DATAPOOL, ONLY : THR, RKIND IMPLICIT NONE INTEGER, INTENT(IN) :: NX1, NX2 REAL(rkind), INTENT(IN) :: X1(NX1), Y1(NX1) REAL(rkind), INTENT(IN) :: X2(NX2) REAL(rkind), INTENT(OUT) :: Y2(NX2) INTEGER :: I, J REAL(rkind) :: DX1(NX1-1) DO I = 1, NX1 - 1 DX1(I) = X1(I+1) - X1(I) END DO DO I = 1, NX2 DO J = 1, NX1 - 1 IF (ABS(X2(I) - X1(J)) .LT. THR) THEN Y2(I) = Y1(J) ELSE IF (X2(I) .GT. X1(J) .AND. X2(I) .LT. X1(J+1)) THEN Y2(I) = Y1(J) + (Y1(J+1)-Y1(J))/DX1(J)*(X2(I)-X1(J)) END IF END DO END DO END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE WWM_ABORT(string) #ifdef MPI_PARALL_GRID USE DATAPOOL, only : parallel_abort, DBG #else USE DATAPOOL, only : DBG #endif IMPLICIT NONE character(*), intent(in) :: string WRITE(DBG%FHNDL, *) TRIM(string) FLUSH(DBG%FHNDL) #ifdef MPI_PARALL_GRID CALL PARALLEL_ABORT(TRIM(string)) #else Print *, 'We have to abort. Reason:' Print *, TRIM(string) STOP 'WWM_ABORT' #endif END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE TEST_FILE_EXIST_DIE(string1, string2) CHARACTER(LEN=*) :: string1 CHARACTER(LEN=*) :: string2 CHARACTER(LEN=512) :: ErrMsg LOGICAL :: LFLIVE ! Print *, 'string1=', TRIM(string1) ! Print *, 'string2=', TRIM(string2) INQUIRE( FILE = TRIM(string2), EXIST = LFLIVE ) IF ( .NOT. LFLIVE ) THEN WRITE(ErrMsg,10) TRIM(string1), TRIM(string2) 10 FORMAT(a, ' ', a) CALL WWM_ABORT(TRIM(ErrMsg)) END IF END SUBROUTINE !********************************************************************** !* Simple Gauss method for solving systems * !* We simply solve the equation A X = Y * !********************************************************************** SUBROUTINE GAUSS_SOLVER(N, EMAT, X, Y) USE DATAPOOL, ONLY: RKIND, ONE IMPLICIT NONE integer, intent(in) :: N real(rkind), intent(inout) :: EMAT(N,N) real(rkind), intent(out) :: X(N) real(rkind), intent(in) :: Y(N) ! real(rkind) :: EMATW(N,N), YW(N) integer U(N), Jsel, I, J, I2 real(rkind) :: alpha, ThePivot, eVal, MaxVal U=0 EMATW=EMAT YW=Y DO I=1,N MaxVal=0 Jsel=-1 DO J=1,N IF (U(J) .eq. 0) THEN eVal=EMATW(J,I) IF (abs(eVal) .gt. MaxVal) THEN MaxVal=abs(eVal) Jsel=J END IF END IF END DO IF (Jsel == -1) THEN CALL WWM_ABORT('Error in GAUSS_SOLVER, singular matrix') END IF ThePivot=ONE/EMATW(Jsel,I) U(Jsel)=I DO I2=I,N EMATW(Jsel,I2)=EMATW(Jsel,I2)*ThePivot END DO YW(Jsel)=YW(Jsel)*ThePivot DO J=1,N IF (J .ne. Jsel) THEN alpha=EMATW(J,I) DO I2=I,N EMATW(J,I2)=EMATW(J,I2) - alpha*EMATW(Jsel,I2) END DO YW(J)=YW(J) - alpha*YW(Jsel) END IF END DO END DO DO Jsel=1,N I=U(Jsel) X(I)=YW(Jsel) END DO END SUBROUTINE !********************************************************************** !* from http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm * !********************************************************************** subroutine SOLVE_TRIDIAG(a,b,c,d,x,n) use datapool, only : rkind implicit none ! a - sub-diagonal (means it is the diagonal below the main diagonal) ! b - the main diagonal ! c - sup-diagonal (means it is the diagonal above the main diagonal) ! d - right part ! x - the answer ! n - number of equations integer,intent(in) :: n real(rkind), intent(in) :: a(n), b(n), c(n), d(n) real(rkind), intent(out) :: x(n) real(rkind) :: cp(n),dp(n) real(rkind) :: m integer i ! initialize c-prime and d-prime cp(1) = c(1)/b(1) dp(1) = d(1)/b(1) ! solve for vectors c-prime and d-prime do i = 2,n m = b(i)-cp(i-1)*a(i) cp(i) = c(i)/m dp(i) = (d(i)-dp(i-1)*a(i))/m enddo ! initialize x x(n) = dp(n) ! solve for x from the vectors c-prime and d-prime do i = n-1, 1, -1 x(i) = dp(i)-cp(i)*x(i+1) end do end subroutine SOLVE_TRIDIAG !********************************************************************** !* * !********************************************************************** SUBROUTINE GAUS1D( M, AMAT, R, AM1, A1M ) USE DATAPOOL, ONLY: RKIND IMPLICIT NONE INTEGER :: M, I, K REAL(rkind) :: AMAT(3,M), R(M), A1M, AM1, FAK REAL(rkind) :: SPALTE(M), ZEILE(M) ! AMAT: TRIDIAGONAL: 1. LEFT DIAGONALE; 2. MAIN DIAGONAL; 3. RIGHT DIAGONAL ! DIMENSION OF MATRIX ! AM1 LOWER LEFT ELEMENT ! A1M UPPER RIGHT ELEMENT ! M ~ EQUATIONS ! R ~ RIGHT HAND SIDE DO K = 1, M ZEILE(K) = 0.0 SPALTE(K)= 0.0 END DO SPALTE(1) = A1M SPALTE(M-1) = AMAT(3,M-1) SPALTE(M) = AMAT(2,M) ZEILE(1) = AM1 ! R DO I = 2, M-1 FAK = AMAT(1,I)/AMAT(2,I-1) R(I)= R(I) - R(I-1)*FAK AMAT(2,I) = AMAT(2,I) - AMAT(3,I-1)*FAK SPALTE(I) = SPALTE(I) - SPALTE(I-1)*FAK END DO ! I = M ZEILE(M-1) = AMAT(1,M) DO K = 2, M-1 FAK = ZEILE(K-1) / AMAT(2,K-1) ZEILE(K) = ZEILE(K) - AMAT(3,K-1)*FAK R(M) = R(M) - R(K-1)*FAK SPALTE(M) = SPALTE(M) - SPALTE(K-1)*FAK END DO K = M FAK = ZEILE(K-1) / AMAT(2,K-1) R(M) = R(M) - R(K-1)*FAK SPALTE(M) = SPALTE(M) - SPALTE(K-1)*FAK ! R I = M - 1 FAK = SPALTE(I) / SPALTE(M) R(I) = R(I) - R(M)*FAK AMAT(3,I) = 0.0 DO I = M-2, 1, -1 FAK = SPALTE(I) / SPALTE(M) R(I) = R(I) - R(M)*FAK AMAT(3,I) = AMAT(3,I) - R(M)*FAK FAK = AMAT(3,I) / AMAT(2,I+1) R(I) = R(I) - R(I+1)*FAK END DO DO I = 1, M-1 R(I) = R(I) / AMAT(2,I) END DO R(M) = R(M) / SPALTE(M) END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE SSORT2 (X, Y, Z, N, KFLAG) USE DATAPOOL, only : rkind !TS: Modified second array z(*) to carry along !***BEGIN PROLOGUE SSORT !***PURPOSE Sort an array and optionally make the same interchanges in ! an auxiliary array. The array may be sorted in increasing ! or decreasing order. A slightly modified QUICKSORT ! algorithm is used. !***LIBRARY SLATEC !***CATEGORY N6A2B !***TYPE SINGLE PRECISION (SSORT-S, DSORT-D, ISORT-I) !***KEYWORDS SINGLETON QUICKSORT, SORT, SORTING !***AUTHOR Jones, R. E., (SNLA) ! Wisniewski, J. A., (SNLA) !***DESCRIPTION ! ! SSORT sorts array X and optionally makes the same interchanges in ! array Y. The array X may be sorted in increasing order or ! decreasing order. A slightly modified quicksort algorithm is used. ! ! Description of Parameters ! X - array of values to be sorted (usually abscissas) ! Y - array to be (optionally) carried along ! N - number of values in array X to be sorted ! KFLAG - control parameter ! = 2 means sort X in increasing order and carry Y along. ! = 1 means sort X in increasing order (ignoring Y) ! = -1 means sort X in decreasing order (ignoring Y) ! = -2 means sort X in decreasing order and carry Y along. ! !***REFERENCES R. C. Singleton, Algorithm 347, An efficient algorithm ! for sorting with minimal storage, Communications of ! the ACM, 12, 3 (1969), pp. 185-187. !***REVISION HISTORY (YYMMDD) ! 761101 DATE WRITTEN ! 761118 Modified to use the Singleton quicksort algorithm. (JAW) ! 890531 Changed all specific intrinsics to generic. (WRB) ! 890831 Modified array declarations. (WRB) ! 891009 Removed unreferenced statement labels. (WRB) ! 891024 Changed category. (WRB) ! 891024 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) ! 901012 Declared all variables; changed X,Y to SX,SY. (M. McClain) ! 920501 Reformatted the REFERENCES section. (DWL, WRB) ! 920519 Clarified error messages. (DWL) ! 920801 Declarations section rebuilt and code restructured to use ! IF-THEN-ELSE-ENDIF. (RWC, WRB) !***END PROLOGUE SSORT ! .. Scalar Arguments .. INTEGER KFLAG, N ! .. Array Arguments .. REAL(rkind) X(*), Y(*), Z(*) ! .. Local Scalars .. REAL(rkind) R, T, TT, TTY, TY, TZ, TTZ INTEGER I, IJ, J, K, KK, L, M, NN ! .. Local Arrays .. INTEGER IL(21), IU(21) ! .. External Subroutines .. ! None ! .. Intrinsic Functions .. INTRINSIC ABS, INT !***FIRST EXECUTABLE STATEMENT SSORT NN = N IF (NN .LT. 1) THEN WRITE (*,*) 'The number of values to be sorted is not positive.' RETURN ENDIF ! KK = ABS(KFLAG) IF (KK.NE.1 .AND. KK.NE.2) THEN WRITE (*,*) 'The sort control parameter, K, is not 2, 1, -1, or -2.' RETURN ENDIF ! ! Alter array X to get decreasing order if needed ! IF (KFLAG .LE. -1) THEN DO 10 I=1,NN X(I) = -X(I) 10 CONTINUE ENDIF ! IF (KK .EQ. 2) GO TO 100 ! ! Sort X only ! M = 1 I = 1 J = NN R = 0.375E0 ! 20 IF (I .EQ. J) GO TO 60 IF (R .LE. 0.5898437E0) THEN R = R+3.90625E-2 ELSE R = R-0.21875E0 ENDIF ! 30 K = I ! ! Select a central element of the array and save it in location T ! IJ = I + INT((J-I)*R) T = X(IJ) ! ! If first element of array is greater than T, interchange with T ! IF (X(I) .GT. T) THEN X(IJ) = X(I) X(I) = T T = X(IJ) ENDIF L = J ! ! If last element of array is less than than T, interchange with T ! IF (X(J) .LT. T) THEN X(IJ) = X(J) X(J) = T T = X(IJ) ! ! If first element of array is greater than T, interchange with T ! IF (X(I) .GT. T) THEN X(IJ) = X(I) X(I) = T T = X(IJ) ENDIF ENDIF ! ! Find an element in the second half of the array which is smaller ! than T ! 40 L = L-1 IF (X(L) .GT. T) GO TO 40 ! ! Find an element in the first half of the array which is greater ! than T ! 50 K = K+1 IF (X(K) .LT. T) GO TO 50 ! ! Interchange these elements ! IF (K .LE. L) THEN TT = X(L) X(L) = X(K) X(K) = TT GO TO 40 ENDIF ! ! Save upper and lower subscripts of the array yet to be sorted ! IF (L-I .GT. J-K) THEN IL(M) = I IU(M) = L I = K M = M+1 ELSE IL(M) = K IU(M) = J J = L M = M+1 ENDIF GO TO 70 ! ! Begin again on another portion of the unsorted array ! 60 M = M-1 IF (M .EQ. 0) GO TO 190 I = IL(M) J = IU(M) ! 70 IF (J-I .GE. 1) GO TO 30 IF (I .EQ. 1) GO TO 20 I = I-1 ! 80 I = I+1 IF (I .EQ. J) GO TO 60 T = X(I+1) IF (X(I) .LE. T) GO TO 80 K = I ! 90 X(K+1) = X(K) K = K-1 IF (T .LT. X(K)) GO TO 90 X(K+1) = T GO TO 80 ! ! Sort X and carry Y along ! 100 M = 1 I = 1 J = NN R = 0.375E0 ! 110 IF (I .EQ. J) GO TO 150 IF (R .LE. 0.5898437E0) THEN R = R+3.90625E-2 ELSE R = R-0.21875E0 ENDIF ! 120 K = I ! ! Select a central element of the array and save it in location T ! IJ = I + INT((J-I)*R) T = X(IJ) TY = Y(IJ) TZ = Z(IJ) ! ! If first element of array is greater than T, interchange with T ! IF (X(I) .GT. T) THEN X(IJ) = X(I) X(I) = T T = X(IJ) Y(IJ) = Y(I) Y(I) = TY TY = Y(IJ) Z(IJ) = Z(I) Z(I) = TZ TZ = Z(IJ) ENDIF L = J ! ! If last element of array is less than T, interchange with T ! IF (X(J) .LT. T) THEN X(IJ) = X(J) X(J) = T T = X(IJ) Y(IJ) = Y(J) Y(J) = TY TY = Y(IJ) Z(IJ) = Z(J) Z(J) = TZ TZ = Z(IJ) ! ! If first element of array is greater than T, interchange with T ! IF (X(I) .GT. T) THEN X(IJ) = X(I) X(I) = T T = X(IJ) Y(IJ) = Y(I) Y(I) = TY TY = Y(IJ) Z(IJ) = Z(I) Z(I) = TZ TZ = Z(IJ) ENDIF ENDIF ! ! Find an element in the second half of the array which is smaller ! than T ! 130 L = L-1 IF (X(L) .GT. T) GO TO 130 ! ! Find an element in the first half of the array which is greater ! than T ! 140 K = K+1 IF (X(K) .LT. T) GO TO 140 ! ! Interchange these elements ! IF (K .LE. L) THEN TT = X(L) X(L) = X(K) X(K) = TT TTY = Y(L) Y(L) = Y(K) Y(K) = TTY TTZ = Z(L) Z(L) = Z(K) Z(K) = TTZ GO TO 130 ENDIF ! ! Save upper and lower subscripts of the array yet to be sorted ! IF (L-I .GT. J-K) THEN IL(M) = I IU(M) = L I = K M = M+1 ELSE IL(M) = K IU(M) = J J = L M = M+1 ENDIF GO TO 160 ! ! Begin again on another portion of the unsorted array ! 150 M = M-1 IF (M .EQ. 0) GO TO 190 I = IL(M) J = IU(M) ! 160 IF (J-I .GE. 1) GO TO 120 IF (I .EQ. 1) GO TO 110 I = I-1 ! 170 I = I+1 IF (I .EQ. J) GO TO 150 T = X(I+1) TY = Y(I+1) TZ = Z(I+1) IF (X(I) .LE. T) GO TO 170 K = I ! 180 X(K+1) = X(K) Y(K+1) = Y(K) Z(K+1) = Z(K) K = K-1 IF (T .LT. X(K)) GO TO 180 X(K+1) = T Y(K+1) = TY Z(K+1) = TZ GO TO 170 ! ! Clean up ! 190 IF (KFLAG .LE. -1) THEN DO 200 I=1,NN X(I) = -X(I) 200 CONTINUE ENDIF END SUBROUTINE !********************************************************************** !* * !********************************************************************** FUNCTION POSITION_BEFORE_POINT(FILEIN) IMPLICIT NONE character(len=*), intent(in) :: FILEIN character(len=1), parameter :: ePoint = '.' integer POSITION_BEFORE_POINT integer len, LPOS, I len=LEN_TRIM(FILEIN) LPOS=-1 DO I=1,len IF (FILEIN(I:I).eq.ePoint) THEN LPOS=I-1 END IF ENDDO IF (LPOS.eq.-1) THEN LPOS=len ENDIF POSITION_BEFORE_POINT=LPOS END FUNCTION !********************************************************************** !* * !********************************************************************** SUBROUTINE GETSTRING (ThePow, TheNb, eStr) IMPLICIT NONE integer, intent(in) :: TheNb, ThePow character(len=ThePow), intent(out) :: eStr character(len=ThePow) :: eStrZero character(len=40) :: eStrNb INTEGER :: ePower, iPower, WeFind, I, nb, IFIRST, h IF (TheNb.lt.0) THEN CALL WWM_ABORT('GETSTRING - call with negative value') END IF ePower=1 WeFind=0 DO I=1,ThePow ePower=ePower*10 IF ((WeFind.eq.0).and.(TheNb.lt.ePower)) THEN WeFind=1 iPower=I END IF END DO IF (WeFind.eq.0) THEN CALL WWM_ABORT('GETSTRING - call with non-representable number') END IF nb=ThePow - iPower DO I=1,ThePow eStrZero(I:I)=' ' END DO DO I=1,nb eStrZero(I:I)='0' END DO WRITE(eStrNb,*) TheNb IFIRST=-1 DO I=1,40 IF ((eStrNb(I:I).ne.' ').and.(IFIRST.eq.-1)) THEN IFIRST=I ENDIF END DO nb=40+1-IFIRST DO I=1,nb h=I+IFIRST-1 eStrNb(I:I)=eStrNb(h:h) END DO DO I=nb+1,40 eStrNb(I:I)=' ' END DO WRITE(eStr,40) TRIM(eStrZero),TRIM(eStrNb) 40 FORMAT (a,a) END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE SPECSTAT(TIME) USE DATAPOOL IMPLICIT NONE REAL, INTENT(IN) :: TIME INTEGER :: IS, ID, IP REAL(RKIND) :: ETOT ! Specstat 2D definiton 1 based on the total energy ... may be to less specific for certain spatial points e.g. coastal points in shadowing areas .... STAT2D = 0. DO IP = 1, MNP ETOT = SUM(AC2(:,:,IP)) DO IS = 1, MSC DO ID = 1, MDC IF (ETOT .GT. THR) THEN IF (AC2(IS,ID,IP)/ETOT .LT. 1.E-5) THEN STAT2D(IS,ID) = STAT2D(IS,ID) + 1./REAL(MNP)*100. ! WRITE(*,'(3I10,4E15.4)') IP, IS, ID, STAT2D(IS,ID), AC2(IS,ID,IP)/ETOT, AC2(IS,ID,IP), ETOT END IF ELSE STAT2D(IS,ID) = STAT2D(IS,ID) + 1./REAL(MNP)*100. END IF END DO END DO END DO DO IS = 1, MSC DO ID = 1, MDC IF (STAT2D(IS,ID) .GT. 99.9) STAT2D(IS,ID) = 0. END DO END DO ! CALL SSORTORG(STAT1D,SIGTMP,MSC,1) ! DO IS = 1, MSC ! DO ID = 1, MDC ! WRITE(*,'(2I10,F15.4)') IS, ID, STAT2D(IS,ID) ! END DO ! WRITE(*,*) IS, STAT1D(IS) ! END DO ! CALL DISPLAY_GRAPH(SPSIG,SUM1D,MSC,MINVAL(SPSIG),MAXVAL(SPSIG),MINVAL(SUM1D),MAXVAL(SUM1D),'','','') END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE INTER_STRUCT_DATA(NDX,NDY,DX,DY,OFFSET_X,OFFSET_Y,MAT,VAL) USE DATAPOOL IMPLICIT NONE INTEGER, INTENT(IN) :: NDX, NDY REAL(rkind), INTENT(IN) :: MAT(NDX,NDY) REAL(rkind), INTENT(IN) :: DX, DY, OFFSET_X, OFFSET_Y REAL(rkind), INTENT(OUT) :: VAL(MNP_WIND) INTEGER :: IP, J_INT, I_INT REAL(rkind) :: WX1, WX2, WX3, WX4, HX1, HX2 REAL(rkind) :: DELTA_X, DELTA_Y, LEN_X, LEN_Y DO IP = 1, MNP_WIND LEN_X = XP_WIND(IP)-OFFSET_X LEN_Y = YP_WIND(IP)-OFFSET_Y I_INT = INT( LEN_X/DX ) + 1 J_INT = INT( LEN_Y/DY ) + 1 DELTA_X = LEN_X - (I_INT - 1) * DX ! Abstand X u. Y DELTA_Y = LEN_Y - (J_INT - 1) * DY ! WX1 = MAT( I_INT , J_INT ) ! Unten Links WX2 = MAT( I_INT , J_INT+1) ! Oben Links WX3 = MAT( I_INT+1, J_INT+1) ! Oben Rechts WX4 = MAT( I_INT+1, J_INT ) ! Unten Rechts IF (IWINDFORMAT == 2) THEN HX1 = WX1 + (WX2-WX1)/DX * DELTA_X HX2 = WX4 + (WX3-WX4)/DX * DELTA_X ELSE IF (IWINDFORMAT == 3) THEN HX1 = WX1 + (WX4-WX1)/DX * DELTA_X HX2 = WX2 + (WX3-WX2)/DX * DELTA_X ELSE CALL WWM_ABORT('Write your HX1/HX2 code here') END IF VAL(IP) = HX1 + (HX2-HX1)/DY * DELTA_Y END DO END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE CreateAngleMatrix(eta_rho, xi_rho, ANG_rho, LON_rho, LAT_rho) implicit none integer, intent(in) :: eta_rho, xi_rho REAL*8, DIMENSION(eta_rho, xi_rho), intent(in) :: LON_rho, LAT_rho REAL*8, DIMENSION(eta_rho, xi_rho), intent(out) :: ANG_rho ! integer eta_u, xi_u, iEta, iXi real*8, allocatable :: LONrad_u(:,:) real*8, allocatable :: LATrad_u(:,:) real*8, allocatable :: azim(:,:) real*8 :: eAzim, fAzim, dlam, eFact1, eFact2 real*8 :: signAzim, signDlam, ThePi, DegTwoRad real*8 :: eLon, eLat, phi1, phi2, xlam1, xlam2 real*8 :: TPSI2, cta12 integer istat eta_u=eta_rho xi_u=xi_rho-1 allocate(LONrad_u(eta_u,xi_u), LATrad_u(eta_u,xi_u), azim(eta_u,xi_u-1), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 33') ThePi=3.141592653589792 DegTwoRad=ThePi/180 DO iEta=1,eta_u DO iXi=1,xi_u eLon=(LON_rho(iEta,iXi)+LON_rho(iEta,iXi+1))*0.5 eLat=(LAT_rho(iEta,iXi)+LAT_rho(iEta,iXi+1))*0.5 LONrad_u(iEta,iXi)=eLon*DegTwoRad LATrad_u(iEta,iXi)=eLat*DegTwoRad END DO END DO DO iEta=1,eta_u DO iXi=1,xi_u-1 phi1=LATrad_u(iEta,iXi) xlam1=LONrad_u(iEta,iXi) phi2=LATrad_u(iEta,iXi+1) xlam2=LONrad_u(iEta,iXi+1) TPSI2=TAN(phi2) dlam=xlam2-xlam1 CALL TwoPiNormalization(dlam) cta12=(cos(phi1)*TPSI2 - sin(phi1)*cos(dlam))/sin(dlam) eAzim=ATAN(1./cta12) CALL MySign(eAzim, signAzim) CALL MySign(dlam, signDlam) IF (signDlam.ne.signAzim) THEN eFact2=1 ELSE eFact2=0 END IF eFact1=-signAzim fAzim=eAzim+ThePi*eFact1*eFact2 azim(iEta,iXi)=fAzim END DO END DO DO iEta=1,eta_u DO iXi=2,xi_u ANG_rho(iEta,iXi)=ThePi*0.5 - azim(iEta,iXi-1) END DO END DO DO iEta=1,eta_u ANG_rho(iEta,1)=ANG_rho(iEta,2) ANG_rho(iEta,xi_rho)=ANG_rho(iEta,xi_u) END DO deallocate(LONrad_u, LATrad_u, azim) END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE CreateAngleMatrix_v(eta_rho,xi_rho,ANG_rho,LON_rho,LAT_rho) USE DATAPOOL, only : rkind, istat implicit none integer, intent(in) :: eta_rho, xi_rho REAL(rkind), DIMENSION(eta_rho, xi_rho), intent(in) :: LON_rho, LAT_rho REAL(rkind), DIMENSION(eta_rho, xi_rho), intent(out) :: ANG_rho ! integer eta_v, xi_v, iEta, iXi real(rkind), allocatable :: LONrad_v(:,:) real(rkind), allocatable :: LATrad_v(:,:) real(rkind), allocatable :: azim(:,:) real(rkind) :: eAzim, fAzim, dlam, eFact1, eFact2 real(rkind) :: signAzim, signDlam, ThePi, DegTwoRad real(rkind) :: eLon, eLat, phi1, phi2, xlam1, xlam2 real(rkind) :: TPSI2, cta12 eta_v=eta_rho-1 xi_v=xi_rho allocate(LONrad_v(eta_v,xi_v), LATrad_v(eta_v,xi_v), azim(eta_v-1,xi_v), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 34') ThePi=3.141592653589792 DegTwoRad=ThePi/180 DO iEta=1,eta_v DO iXi=1,xi_v eLon=(LON_rho(iEta,iXi)+LON_rho(iEta+1,iXi))*0.5 eLat=(LAT_rho(iEta,iXi)+LAT_rho(iEta+1,iXi))*0.5 LONrad_v(iEta,iXi)=eLon*DegTwoRad LATrad_v(iEta,iXi)=eLat*DegTwoRad END DO END DO DO iEta=1,eta_v-1 DO iXi=1,xi_v phi1=LATrad_v(iEta,iXi) xlam1=LONrad_v(iEta,iXi) phi2=LATrad_v(iEta+1,iXi) xlam2=LONrad_v(iEta+1,iXi) TPSI2=TAN(phi2) dlam=xlam2-xlam1 CALL TwoPiNormalization(dlam) cta12=(cos(phi1)*TPSI2 - sin(phi1)*cos(dlam))/sin(dlam) eAzim=ATAN(1./cta12) CALL MySign(eAzim, signAzim) CALL MySign(dlam, signDlam) IF (signDlam.ne.signAzim) THEN eFact2=1 ELSE eFact2=0 END IF eFact1=-signAzim fAzim=eAzim+ThePi*eFact1*eFact2 azim(iEta,iXi)=fAzim END DO END DO DO iEta=2,eta_v DO iXi=1,xi_v ANG_rho(iEta,iXi)=ThePi*0.5 - azim(iEta-1,iXi) END DO END DO DO iXi=1,xi_v ANG_rho(1,iXi)=ANG_rho(2,iXi) ANG_rho(eta_rho,iXi)=ANG_rho(eta_v,iXi) END DO deallocate(LONrad_v) deallocate(LATrad_v) deallocate(azim) END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE TwoPiNormalization(TheAng) USE DATAPOOL, only : rkind implicit none REAL(rkind), intent(inout) :: TheAng ! REAL(rkind) :: ThePi ThePi=3.141592653589792 IF (TheAng < -ThePi) THEN TheAng=TheAng + 2*ThePi ENDIF IF (TheAng > ThePi) THEN TheAng=TheAng - 2*ThePi ENDIF END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE MySign(TheVal, TheSign) USE DATAPOOL, only : rkind implicit none REAL(rkind), intent(in) :: TheVal REAL(rkind), intent(out) :: TheSign IF (TheVal > 0) THEN TheSign=1 ELSEIF (TheVal < 0) THEN TheSign=-1 ELSE TheSign=0 END IF END SUBROUTINE !********************************************************************** !* * !**********************************************************************