#include "wwm_functions.h" #define DEBUG #undef DEBUG !********************************************************************** !* * !********************************************************************** SUBROUTINE SPECTRAL_SHAPE(SPPAR,ACLOC,LDEBUG,CALLFROM, OPTI) USE DATAPOOL IMPLICIT NONE REAL(rkind), INTENT(OUT) :: ACLOC(MSC,MDC) REAL(rkind), INTENT(INOUT) :: SPPAR(8) CHARACTER(LEN=*), INTENT(IN) :: CALLFROM LOGICAL, INTENT(IN) :: LDEBUG, OPTI IF (OPTI) THEN CALL OPTI_SPECTRAL_SHAPE(SPPAR,ACLOC,LDEBUG,CALLFROM) ELSE CALL KERNEL_SPECTRAL_SHAPE(SPPAR,ACLOC,LDEBUG,CALLFROM) END IF END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE COMPUTE_ESTIMATE_PER_DIR_SHAPE(SPPAR, ACLOC, HS, TM, DM) USE DATAPOOL IMPLICIT NONE REAL(rkind), INTENT(IN) :: SPPAR(8) REAL(rkind), INTENT(IN) :: ACLOC(MSC,MDC) REAL(rkind), INTENT(OUT) :: HS, TM, DM REAL(rkind) :: DEPLOC, CURTXYLOC(2) REAL(rkind) :: WKLOC(MSC) REAL(rkind) :: FPP, TPP, CPP, WNPP, CGPP, KPP, LPP REAL(rkind) :: PEAKDSPR, PEAKDM, DPEAK, TPPD, KPPD, CGPD, CPPD REAL(rkind) :: TM01, TM02, TM10, KLM, WLM REAL(rkind) :: ETOTS, ETOTC, DSPR REAL(rkind) :: SPSIGLOC, WVN, WVC, WVK, WVCG integer ISMAX, IS DEPLOC=10000 ISMAX=MSC CURTXYLOC=ZERO DO IS=1,MSC SPSIGLOC = SPSIG(IS) CALL WAVEKCG(DEPLOC,SPSIGLOC,WVN,WVC,WVK,WVCG) WKLOC(IS)=WVK END DO CALL MEAN_PARAMETER_LOC(ACLOC,CURTXYLOC,DEPLOC,WKLOC,ISMAX,HS,TM01,TM02,TM10,KLM,WLM) IF (SPPAR(5) .gt. 0) THEN CALL PEAK_PARAMETER_LOC(ACLOC,DEPLOC,ISMAX,FPP,TPP,CPP,WNPP,CGPP,KPP,LPP,PEAKDSPR,PEAKDM,DPEAK,TPPD,KPPD,CGPD,CPPD) TM=TPP DM=PEAKDM ELSE CALL MEAN_DIRECTION_AND_SPREAD_LOC(ACLOC,ISMAX,ETOTS,ETOTC,DM,DSPR) TM=TM01 END IF END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE OPTI_SPECTRAL_SHAPE(SPPAR,ACLOC,LDEBUG,CALLFROM) USE DATAPOOL IMPLICIT NONE REAL(rkind), INTENT(OUT) :: ACLOC(MSC,MDC) REAL(rkind), INTENT(INOUT) :: SPPAR(8) CHARACTER(LEN=*), INTENT(IN) :: CALLFROM LOGICAL, INTENT(IN) :: LDEBUG REAL(rkind) :: HS, TM, DM, TheErr, DeltaPer, Tper REAL(rkind) :: DiffAng, DEG, ADIR REAL(rkind) :: SPPARwork1(8), SPPARwork2(8), SPPARwork(8) integer :: iIter, nbIter, eSign, IS, ID REAL(rkind) :: eSum IF (ABS(SPPAR(5)) .eq. 3) THEN CALL KERNEL_SPECTRAL_SHAPE(SPPAR,ACLOC,LDEBUG,CALLFROM) RETURN END IF SPPARwork1=SPPAR SPPARwork2=SPPAR Tper=SPPAR(2) CALL KERNEL_SPECTRAL_SHAPE(SPPAR,ACLOC,LDEBUG,CALLFROM) CALL COMPUTE_ESTIMATE_PER_DIR_SHAPE(SPPAR, ACLOC, HS, TM, DM) DeltaPer=Tper - TM IF (TM < Tper) THEN eSign=1 ELSE eSign=-1 END IF SPPARwork=SPPAR SPPARwork(2)=SPPAR(2) + DeltaPer DO CALL KERNEL_SPECTRAL_SHAPE(SPPARwork,ACLOC,LDEBUG,CALLFROM) CALL COMPUTE_ESTIMATE_PER_DIR_SHAPE(SPPAR, ACLOC, HS, TM, DM) TheErr=(TM - Tper)*eSign IF (TheErr > 0) THEN EXIT END IF SPPARwork(2)=SPPARwork(2) + DeltaPer END DO IF (eSign .eq. 1) THEN SPPARwork1=SPPAR SPPARwork2=SPPARwork ELSE SPPARwork1=SPPARwork SPPARwork2=SPPAR END IF nbIter=20 iIter=0 DO SPPARwork=0.5_rkind*SPPARwork1 + 0.5_rkind*SPPARwork2 CALL KERNEL_SPECTRAL_SHAPE(SPPARwork,ACLOC,LDEBUG,CALLFROM) CALL COMPUTE_ESTIMATE_PER_DIR_SHAPE(SPPAR, ACLOC, HS, TM, DM) IF (TM > Tper) THEN SPPARwork2=SPPARwork ELSE SPPARwork1=SPPARwork END IF iIter=iIter + 1 IF (iIter > nbIter) THEN EXIT END IF END DO CALL KERNEL_SPECTRAL_SHAPE(SPPARwork,ACLOC,LDEBUG,CALLFROM) CALL COMPUTE_ESTIMATE_PER_DIR_SHAPE(SPPAR, ACLOC, HS, TM, DM) SPPARwork(1)=SPPAR(1)*(SPPAR(1)/HS) CALL KERNEL_SPECTRAL_SHAPE(SPPARwork,ACLOC,LDEBUG,CALLFROM) DO IS=1,MSC eSum=sum(ACLOC(IS,:)) END DO CALL DEG2NAUT (SPPAR(3), DEG, LNAUTIN) ADIR = DEG * DEGRAD DO ID=1,MDC eSum=sum(ACLOC(:,ID)) DiffAng=(360.0_rkind/PI2)*(SPDIR(ID) - ADIR) END DO END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE KERNEL_SPECTRAL_SHAPE(SPPAR,ACLOC,LDEBUG,CALLFROM) USE DATAPOOL IMPLICIT NONE REAL(rkind), INTENT(OUT) :: ACLOC(MSC,MDC) REAL(rkind), INTENT(INOUT) :: SPPAR(8) CHARACTER(LEN=*), INTENT(IN) :: CALLFROM LOGICAL, INTENT(IN) :: LDEBUG INTEGER ID, IS, LSHAPE, ITPER, ISP, ISIGMP REAL(rkind) :: APSHAP, AUX1, AUX2, AUX3, AM0, AM1, AS2, AS3, ETOT, VEC2DEG REAL(rkind) :: COEFF, SYF , MPER, PKPER, DIFPER, EHFR REAL(rkind) :: MS, DEG, ETOTS, ETOTC, FF, CPSHAP, PPSHAP, DM, EAD, DS REAL(rkind) :: RA, SALPHA, SF, SF4, SF5, FPK, FPK4, EFTAIL, CDIR REAL(rkind) :: GAMMA_FUNC, DSPR, AACOS, ADIR, EPTAIL, APTAIL, PPTAIL REAL(rkind) :: OMEG, EFTOT, ETOTT, CTOT, TM1, TPEAK LOGICAL LOGPM, LINCOUT ! SPPARM(1), WBHS: Hs, sign. wave height ! SPPARM(2), WBTP: Wave period given by the user (either peak or mean) ! SPPARM(3), WBDM: average direction ! SPPARM(4), WBDS: directional spread ! SPPARM(5), WBSS: spectral shape (1-4), (1 - Pierson-Moskowitz, 2 - JONSWAP, 3 - BIN, 4 - Gauss) peak (+) or mean frequency (-) ! SPPARM(6), WBDSMS: directional spreading in degree (2) or exponent (1) ! SPPARM(7), WBGAUSS: gaussian width for the gauss spectrum 0.1 ! SPPARM(8), WBPKEN: peak enhancement factor for the JONSWAP spectra 3.3 IF (LDEBUG) THEN WRITE(DBG%FHNDL,*) 'HS PER DIR DPSR SHAPE DEGEXP GAUSS PEAK' WRITE(DBG%FHNDL,'(8F10.4)') SPPAR(8) ENDIF ETOT = 0. EFTOT = 0. ACLOC = 0. IF (SPPAR(1) .LT. THR .OR. SPPAR(2) .LT. THR .OR. SPPAR(4) .LT. THR) THEN ACLOC = 0. RETURN END IF IF (SPPAR(5) .LT. 0) THEN LSHAPE = -INT(SPPAR(5)) LOGPM = .FALSE. ELSE LSHAPE = INT(SPPAR(5)) LOGPM = .TRUE. ENDIF ! PKPER = SPPAR(2) ITPER = 0 IF (LSHAPE.EQ.3) THEN ! select bin closest to given period DIFPER = 1.E10 DO IS = 1, MSC IF (ABS(PKPER - PI2/SPSIG(IS)) .LT. DIFPER) THEN ISP = IS DIFPER = ABS(PKPER - PI2/SPSIG(IS)) END IF ENDDO ENDIF ! 100 FPK = (1./PKPER) FPK4 = FPK**4 IF (LSHAPE.EQ.1) THEN SALPHA = SPPAR(1)**2*FPK4 * 5./16. ELSE IF (LSHAPE.EQ.2) THEN SALPHA = (SPPAR(1)**2 * FPK4) / ((0.06533*(SPPAR(8)**0.8015)+0.13467)*16.) ELSE IF (LSHAPE.EQ.4) THEN AUX1 = SPPAR(1)**2 / ( 16.* SQRT (PI2) * SPPAR(7)) AUX3 = 2._rkind * SPPAR(7)**2 ENDIF ! DO IS = 1, MSC ! IF (LSHAPE.EQ.1) THEN SF = SPSIG(IS) / PI2 SF4 = SF**4 SF5 = SF**5 RA = (SALPHA/SF5)*EXP(-(5.*FPK4)/(4.*SF4))/(PI2*SPSIG(IS)) ACLOC(IS,MDC) = RA ELSE IF (LSHAPE.EQ.2) THEN SF = SPSIG(IS)/(PI2) SF4 = SF**4 SF5 = SF**5 CPSHAP = 1.25_rkind * FPK4 / SF4 IF (CPSHAP.GT.10._rkind) THEN RA = 0._rkind ELSE RA = (SALPHA/SF5) * EXP(-CPSHAP) ENDIF IF (SF .LT. FPK) THEN COEFF = 0.07_rkind ELSE COEFF = 0.09_rkind ENDIF IF (COEFF*FPK .GT. SMALL) THEN APSHAP = 0.5_rkind * ((SF-FPK) / (COEFF*FPK))**2 ELSE APSHAP = ZERO ENDIF IF (APSHAP .GT. 10._rkind) THEN SYF = 1. ELSE PPSHAP = EXP(-APSHAP) SYF = SPPAR(8)**PPSHAP ENDIF RA = SYF*RA/(SPSIG(IS)*PI2) ACLOC(IS,MDC) = RA IF (LDEBUG) WRITE(DBG%FHNDL,*) 'IS LOOP', IS, SF, FPK, SYF, RA ! ELSE IF (LSHAPE .EQ. 3) THEN IF (IS.EQ.ISP) THEN ISBIN = ISP ACLOC(IS,MDC) = ( SPPAR(1)**2 ) / ( 16. * SIGPOW(IS,2) * FRINTF ) ELSE ACLOC(IS,MDC) = 0. END IF ELSE IF (LSHAPE .EQ. 4) THEN AUX2 = ( SPSIG(IS) - ( PI2 / PKPER ) )**2 RA = AUX1 * EXP ( -1. * AUX2 / AUX3 ) / SPSIG(IS) ACLOC(IS,MDC) = RA ELSE CALL WWM_ABORT('Wrong type for frequency shape 1 - 4') ENDIF END DO MPER = 0. IF (.NOT.LOGPM.AND.ITPER.LT.100) THEN ITPER = ITPER + 1 ! calculate average frequency AM0 = 0. AM1 = 0. DO IS = 1, MSC AS2 = ACLOC(IS,MDC) * SIGPOW(IS,2) AS3 = AS2 * SPSIG(IS) AM0 = AM0 + AS2 AM1 = AM1 + AS3 ENDDO ! contribution of tail to total energy density PPTAIL = PTAIL(1) - 1. APTAIL = 1. / (PPTAIL * (1. + PPTAIL * (FRINTH-1.))) AM0 = AM0 * FRINTF + APTAIL * AS2 PPTAIL = PTAIL(1) - 2. EPTAIL = 1. / (PPTAIL * (1. + PPTAIL * (FRINTH-1.))) AM1 = AM1 * FRINTF + EPTAIL * AS3 ! Mean period: IF ( AM1.GT.THR) THEN MPER = PI2 * AM0 / AM1 ELSE MPER = ZERO ENDIF ! write(*,'(I10,4F15.8)') ITPER, MPER, & ! & ABS(MPER-SPPAR(2)) .GT. 0.01*SPPAR(2), & ! & (SPPAR(2) / MPER) * PKPER IF (ABS(MPER-SPPAR(2)) .GT. 0.01*SPPAR(2) .AND. MPER .GT. THR) THEN ! modification suggested by Mauro Sclavo PKPER = (SPPAR(2)/MPER) * PKPER GOTO 100 ELSE PKPER = ZERO ENDIF ELSE IF (ITPER.GE.100) THEN WRITE (STAT%FHNDL,*) 'No convergence calculating the spectrum' FLUSH(STAT%FHNDL) ENDIF CALL DEG2NAUT (SPPAR(3), DEG, LNAUTIN) ADIR = DEG * DEGRAD IF (INT(SPPAR(6)) .EQ. 1) THEN DSPR = PI * SPPAR(4) / 180._rkind MS = MAX (DSPR**(-2) - TWO, 1._rkind) ELSE MS = SPPAR(4) ENDIF IF (MS.LT.12._rkind) THEN CTOT = (2._rkind**MS)*(GAMMA_FUNC(0.5_rkind*MS+1._rkind))**2/(PI*GAMMA_FUNC(MS+1._rkind)) ELSE CTOT = SQRT(0.5_rkind*MS/PI)/(1._rkind-0.25_rkind/MS) ENDIF LINCOUT = .FALSE. DO ID = 1, MDC AACOS = COS(SPDIR(ID) - ADIR) !write(*,*) aacos, spdir(id), adir IF (AACOS .GT. 0._rkind) THEN CDIR = CTOT * MAX (AACOS**MS, THR) IF (.NOT.LCIRD) THEN IF (AACOS .GE. COS(DDIR)) THEN LINCOUT = .TRUE. !WRITE(*,*) 'ERROR', AACOS, COS(DDIR) !STOP 'AVERAGE DIRECTION IS OUT OF SECTOR' ENDIF ENDIF ELSE CDIR = 0._rkind ENDIF DO IS = 1, MSC ACLOC(IS,ID) = CDIR * ACLOC(IS,MDC) !write(*,'(2I10,2F15.8)') is, id, cdir, ACLOC(IS,MDC) ENDDO ENDDO ! ! Finished creating ! ! Integral Parameters of the Input Spectra ... ! AR: 2DO: Optimize the output routine for this task ... ! IF (LDEBUG) THEN ETOT = 0. ETOTC = 0. ETOTS = 0. EFTAIL = 1.0 / (PTAIL(1)-1.0) IF (MSC .GE. 2) THEN 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 ETOTC = ETOTC + EAD * COS(SPDIR(ID)) ETOTS = ETOTS + EAD * SIN(SPDIR(ID)) END DO IF (MSC > 3) THEN EHFR = ACLOC(MSC,ID) * SPSIG(MSC) ETOT = ETOT + DDIR * EHFR * SPSIG(MSC) * EFTAIL ENDIF END DO ELSE DS = SGHIGH - SGLOW DO ID = 1, MDC EAD = ACLOC(1,ID) * DS * DDIR ETOT = ETOT + EAD ETOTC = ETOTC + EAD * COS(SPDIR(ID)) ETOTS = ETOTS + EAD * SIN(SPDIR(ID)) END DO END IF IF (ETOT > THR ) THEN DM = VEC2DEG (ETOTC, ETOTS) CALL DEG2NAUT(DM,DEG,LNAUTOUT) DM = DEG FF = MIN (ONE, SQRT(ETOTC*ETOTC+ETOTS*ETOTS)/ETOT) DSPR = SQRT(2.-2.*FF) * 180./PI ELSE DM = 0. DSPR = 0. END IF ETOTT = 0.0 EFTOT = 0.0 EFTAIL = PTAIL(3) DO ID = 1, MDC DO IS = 1, MSC OMEG = SPSIG(IS) EAD = FRINTF * SIGPOW(IS,2) * ACLOC(IS,ID) ETOTT = ETOTT + EAD EFTOT = EFTOT + EAD * OMEG END DO END DO IF (EFTOT > VERYSMALL) THEN TM1 = PI2*ETOTT/EFTOT ELSE TM1 = 0.0 END IF ETOTT = 0.0 ISIGMP = -1 DO IS = 1, MSC EAD = 0.0 DO ID = 1, MDC EAD = EAD + SPSIG(IS)*ACLOC(IS,ID)*DDIR ENDDO IF (EAD > ETOTT) THEN ETOTT = EAD ISIGMP = IS END IF END DO IF (ISIGMP > 0) THEN TPEAK = 1./(SPSIG(ISIGMP)/PI2) ELSE TPEAK = 0.0 END IF WRITE (STAT%FHNDL,'("+TRACE...",A)') 'GIVEN BOUNDARY SPECTRA AND RECALCULATED WAVE PARAMETERS' WRITE (STAT%FHNDL,'("+TRACE...",A)') 'THE DIFFERENCE IS MOSTLY DUE TO COARSE RESOLUTION IN SPECTRAL SPACE' WRITE (STAT%FHNDL,*) 'GIVEN ', 'HS =', SPPAR(1) WRITE (STAT%FHNDL,*) 'GIVEN ', 'DM =', SPPAR(3) WRITE (STAT%FHNDL,*) 'GIVEN ', 'DSPR =', SPPAR(4) WRITE (STAT%FHNDL,*) 'GIVEN ', 'TM or TP', SPPAR(2) WRITE (STAT%FHNDL,*) 'SIMUL ', 'HS =', 4. * SQRT(ETOT) WRITE (STAT%FHNDL,*) 'SIMUL ', 'DM =', DM WRITE (STAT%FHNDL,*) 'SIMUL ', 'DSPR =', DSPR WRITE (STAT%FHNDL,*) 'SIMUL ', 'TM=', TM1, 'TPEAK=', TPEAK WRITE (STAT%FHNDL,*) 'TOT AC =', SUM(ACLOC) WRITE (STAT%FHNDL,*) SPPAR FLUSH(STAT%FHNDL) END IF END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE SPECTRUM_INT(ACLOC) USE DATAPOOL IMPLICIT NONE REAL(rkind), INTENT(OUT) :: ACLOC(MSC,MDC) REAL(rkind) :: MS(MSC), MS1, ADIR1, DS, EAD REAL(rkind) :: INSPF(WBMSC) REAL(rkind) :: INSPE(WBMSC) REAL(rkind) :: INDIR(WBMSC) REAL(rkind) :: INSPRD(WBMSC) REAL(rkind) :: INMS(WBMSC) REAL(rkind) :: SPCDIR(MSC) INTEGER :: IS, IS2, ID REAL(rkind) :: CTOT(MSC), CDIRT, CDIR(MDC), CTOT1, CDIR1 REAL(rkind) :: DDACOS, DEG, DX, DIFFDX, YINTER REAL(rkind) :: GAMMA_FUNC, ETOT, TM2 REAL(rkind) :: EFTOT, TM1, OMEG, PPTAIL, OMEG2 REAL(rkind) :: RA, ETAIL, EFTAIL REAL(rkind) :: THD(MDC) REAL(rkind) :: DTHD, RTH0 MS1 = WBDS ADIR1 = WBDM PSHAP(1) = 3.3 PSHAP(2) = 0.07 PSHAP(3) = 0.09 PSHAP(4) = 0.01 PSHAP(5) = 2.0 PSHAP(6) = 0.01 ! ! 2do Convert ... set the correct units ... ! DO IS = 1, WBMSC INSPF(IS) = SFRQ(IS,1) * PI2 INDIR(IS) = SDIR(IS,1) INSPRD(IS) = SPRD(IS,1) * DEGRAD INSPE(IS) = SPEG(IS,1,1) / PI2 END DO ETOT = 0.0 DO IS = 2, WBMSC DS = (INSPF(IS) - INSPF(IS-1)) EAD = 0.5*(INSPE(IS) + INSPE(IS-1))*DS ETOT = ETOT + EAD END DO WRITE (STAT%FHNDL,*) 'HS - INPUTSPECTRA - 1', 4.0*SQRT(ETOT) ACLOC = 0. CALL INTERLIN (WBMSC, MSC, INSPF, SPSIG, INSPE, ACLOC(:,1)) DO IS = 1, MSC IF (SPSIG(IS) .GT. INSPF(WBMSC)) THEN WRITE (STAT%FHNDL,*) 'Discrete Frequency is bigger then measured set FRMAX =', INSPF(WBMSC)/PI2 WRITE (STAT%FHNDL,*) 'Setting all Action above the max. measured freq. zero' ACLOC(IS,1) = 0.0 END IF END DO ETOT = 0.0 DO IS = 2, MSC DS = SPSIG(IS) - SPSIG(IS-1) EAD = 0.5*(ACLOC(IS,1) + ACLOC(IS-1,1))*DS ETOT = ETOT + EAD END DO WRITE (STAT%FHNDL,*) 'HS - INPUTSPECTRA - INTERPOLATED', 4.0*SQRT(ETOT) ! ! Convert to Wave Action if nessasary ! IF (LWBAC2EN) THEN DO IS = 1, MSC ACLOC(IS,1) = ACLOC(IS,1) / SPSIG(IS) END DO END IF ETOT = 0.0 DO IS = 2, MSC DS = SPSIG(IS) - SPSIG(IS-1) EAD = 0.5 * (SPSIG(IS)*ACLOC(IS,1)+SPSIG(IS-1)*ACLOC(IS-1,1))*DS ETOT = ETOT + EAD END DO WRITE (STAT%FHNDL,*) 'HS - INPUTSPECTRA - WAVE ACTION', 4.0*SQRT(ETOT) ! ! Convert from nautical to cartesian direction if nessasary and from deg2rad ! DO IS = 1, WBMSC CALL DEG2NAUT (INDIR(IS), DEG, LNAUTIN) INDIR(IS) = DEG END DO ! ! Interpolate Directions in frequency space ! DO IS = 1, MSC DO IS2 = 1, WBMSC - 1 IF (SPSIG(IS) .GT. INSPF(IS2) .AND. SPSIG(IS) .LT. INSPF(IS2+1)) THEN DX = INSPF(IS2+1) - INSPF(IS2) DIFFDX = SPSIG(IS) - INSPF(IS2) CALL INTERDIR( INDIR(IS2), INDIR(IS2+1), DX, DIFFDX, YINTER) SPCDIR(IS) = YINTER IF (SPSIG(IS) .GT. INSPF(WBMSC) ) SPCDIR(IS) = 0.0 END IF END DO IF (SPSIG(IS) .GT. INSPF(WBMSC) ) SPCDIR(IS) = 0.0 END DO CALL INTERLIN (WBMSC, MSC, INSPF, SPSIG, INDIR, SPCDIR) DO IS = 1, MSC IF ( SPSIG(IS) .GT. INSPF(WBMSC) ) SPCDIR(IS) = 0. END DO DO IS = 1, MSC DEG = SPCDIR(IS) * DEGRAD SPCDIR(IS) = DEG END DO IF (LINDSPRDEG) THEN DO IS = 1, WBMSC INMS(IS) = MAX (INSPRD(IS)**(-2) - TWO, ONE) END DO ELSE DO IS = 1, WBMSC INMS(IS) = INSPRD(IS) END DO END IF ! ! Interpolate MS in Frequency Space, if LCUBIC than Cubic Spline Interpolation is used ! MS = 0. CALL INTERLIN (WBMSC, MSC, INSPF, SPSIG, INMS, MS) DO IS = 1, MSC IF ( SPSIG(IS) .GT. INSPF(WBMSC) ) MS(IS) = MS(IS-1) END DO ! ! Construction of a JONSWAP Spectra ! IF (LPARMDIR) THEN IF (MS1.GT.10.) THEN CTOT1 = SQRT(MS1/(2.*PI)) * (1. + 0.25/MS1) ELSE CTOT1 = 2.**MS1 * (GAMMA_FUNC(1.+0.5*MS1))**2 / (PI * GAMMA_FUNC(1.+MS1)) ENDIF DO IS = 1, MSC RA = ACLOC(IS,1) CDIRT = 0. DO ID = 1, MDC DDACOS = COS(SPDIR(ID) - ADIR1) IF (DDACOS .GT. 0.) THEN CDIR1 = CTOT1 * MAX (DDACOS**MS1, THR) ELSE CDIR1 = 0. ENDIF ACLOC(IS,ID) = RA * CDIR1 CDIRT = CDIRT + CDIR1 * DDIR END DO END DO ELSE DO IS = 1, MSC IF (MS(IS).GT.10.) THEN CTOT(IS) = SQRT(MS(IS)/(2.*PI)) * (1. + 0.25/MS(IS)) ELSE CTOT(IS) = 2.**MS(IS) * (GAMMA_FUNC(1.+0.5*MS(IS)))**2.0 / (PI * GAMMA_FUNC(1.+MS(IS))) ENDIF END DO DO IS = 1, MSC RA = ACLOC(IS,1) CDIRT = 0. DO ID = 1, MDC DDACOS = COS(SPDIR(ID) - SPCDIR(IS)) IF (DDACOS .GT. 0.) THEN CDIR(ID) = CTOT(IS) * MAX (DDACOS**MS(IS), ZERO) ELSE CDIR(ID) = 0. ENDIF ACLOC(IS,ID) = RA * CDIR(ID) CDIRT = CDIRT + CDIR(ID) * DDIR ENDDO IF (100. - 1./CDIRT*100. .GT. 1.) THEN WRITE (STAT%FHNDL,*) 100 - 1./CDIRT*100., 'ERROR BIGGER THAN 1% IN THE DIRECTIONAL DISTTRIBUTION' WRITE (STAT%FHNDL,*) 'PLEASE CHECK THE AMOUNG OF DIRECTIONAL BINS AND THE PRESCRIBED DIRECTIONAL SPREADING' END IF ENDDO END IF ETOT = 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 WRITE (STAT%FHNDL,*) 'HS - INPUTSPECTRA - AFTER 2D', 4.0*SQRT(ETOT) ETOT = 0. EFTOT = 0. PPTAIL = PTAIL(1) - 1. ETAIL = 1. / (PPTAIL * (1. + PPTAIL * (FRINTH-1.))) PPTAIL = PTAIL(1) - 2. EFTAIL = 1. / (PPTAIL * (1. + PPTAIL * (FRINTH-1.))) DO ID=1, MDC DO IS = 1, MSC OMEG = SPSIG(IS) EAD = FRINTF * SIGPOW(IS,2) * ACLOC(IS,ID) ETOT = ETOT + EAD EFTOT = EFTOT + EAD * OMEG ENDDO IF (MSC .GT. 3) THEN EAD = SIGPOW(MSC,2) * ACLOC(MSC,ID) ETOT = ETOT + ETAIL * EAD EFTOT = EFTOT + EFTAIL * OMEG * EAD ! WRITE (*,*) ETAIL * EAD, EFTAIL * OMEG * EAD, ACLOC(MSC,ID) ENDIF ENDDO IF (EFTOT.GT.THR) THEN TM1 = PI2 * ETOT / EFTOT ELSE TM1 = 0. ENDIF ETOT = 0. EFTOT = 0. PPTAIL = PTAIL(1) - 1. ETAIL = 1. / (PPTAIL * (1. + PPTAIL * (FRINTH-1.))) PPTAIL = PTAIL(1) - 3. EFTAIL = 1. / (PPTAIL * (1. + PPTAIL * (FRINTH-1.))) DO ID=1, MDC DO IS=1,MSC EAD = SIGPOW(IS,2) * ACLOC(IS,ID) * FRINTF OMEG2 = SIGPOW(IS,2) ETOT = ETOT + EAD EFTOT = EFTOT + EAD * OMEG2 ENDDO IF (MSC .GT. 3) THEN EAD = SIGPOW(MSC,2) * ACLOC(MSC,ID) ETOT = ETOT + ETAIL * EAD EFTOT = EFTOT + EFTAIL * OMEG2 * EAD ENDIF ENDDO IF (EFTOT .GT. THR) THEN TM2 = PI2 * SQRT(ETOT/EFTOT) ELSE TM2 = 0. END IF ETOT = 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 WRITE (STAT%FHNDL,*) 'HS - INPUTSPECTRA - AFTER 2D', 4.0*SQRT(ETOT) WRITE (STAT%FHNDL,*) 'TM01, TM02 & HS', TM1, TM2, 4.0*SQRT(ETOT) !FLUSH(DBG%FHNDL) !FLUSH(STAT%FHNDL) IF (.FALSE.) THEN ! Write WW3 spectra of the input boundary condition ... DTHD=360./MDC RTH0=SPDIR(1)/DDIR DO ID = 1, MDC THD(ID)=DTHD*(RTH0+MyREAL(ID-1)) END DO WRITE (4001,1944) 'WAVEWATCH III SPECTRA', MSC, MDC, 1, 'LAI ET AL' WRITE (4001,1945) (SPSIG(IS)*INVPI2,IS=1,MSC) WRITE (4001,1946) (MOD(2.5*PI-SPDIR(ID),PI2),ID=1,MDC) WRITE (4001,901) 'LAI SPEC', 0., 0., 0., 0., 0., 0., 0. WRITE (4001,902) ((ACLOC(IS,ID)*SPSIG(IS)/PI2*RHOW*G9,IS=1,MSC),ID=1,MDC) END IF 901 FORMAT ('''',A10,'''',2F7.2,F10.1,2(F7.2,F6.1)) 902 FORMAT (7E11.3) 1943 FORMAT ( ' File name : ',A,' (',A,')') 1944 FORMAT ('''',A,'''',1X,3I6,1X,'''',A,'''') 1945 FORMAT (8E10.3) 1946 FORMAT (7E11.3) END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE SET_WAVE_BOUNDARY USE DATAPOOL IMPLICIT NONE INTEGER :: IP, eIdx IF (LBCWA .OR. LBCSP) THEN ! Spectrum or parametric boundary condition IF (LINHOM) THEN DO IP = 1, IWBMNP eIdx = IWBNDLC(IP) AC2(:,:,eIdx) = WBAC(:,:,IP) END DO ELSE DO IP = 1, IWBMNP eIdx = IWBNDLC(IP) AC2(:,:,eIdx) = WBAC(:,:,1) END DO END IF END IF END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE COMPUTE_IFILE_IT(IFILE, IT) USE DATAPOOL IMPLICIT NONE INTEGER, INTENT(OUT) :: IFILE, IT REAL(rkind) :: DTMP INTEGER ITMP DTMP = (MAIN%TMJD-BND_TIME_ALL_FILES(1,1)) * DAY2SEC ITMP = 0 DO IFILE = 1, NUM_NETCDF_FILES_BND ITMP = ITMP + NDT_BND_FILE(IFILE) IF (ITMP .GT. INT(DTMP/SEBO%DELT)) EXIT END DO ITMP = SUM(NDT_BND_FILE(1:IFILE-1)) IT = NINT(DTMP/SEBO%DELT) - ITMP + 1 IF (IT .GT. NDT_BND_FILE(IFILE)) THEN IFILE = IFILE + 1 IT = 1 ENDIF END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE WAVE_BOUNDARY_CONDITION(WBACOUT) USE DATAPOOL IMPLICIT NONE REAL(rkind), INTENT(OUT) :: WBACOUT(MSC,MDC,IWBMNP) INTEGER :: IP !AR: WAVE BOUNDARY ! SPPARM(1): Hs, sign. wave height ! SPPARM(2): Wave period given by user (either peak or mean) ! SPPARM(3): average direction ! SPPARM(4): directional spread ! SPPARM(5): spectral shape (1-4), ! 1 - Pierson-Moskowitz ! 2 - JONSWAP ! 3 - BIN ! 4 - Gauss ! positive peak (+) or mean frequency (-) ! SPPARM(6): directional spreading in degree (1) or exponent (2) ! SPPARM(7): gaussian width for the gauss spectrum 0.1 ! SPPARM(8): peak enhancement factor for the JONSWAP spectra 3. ! ! Count number of active boundary points ... ! IF(LWW3GLOBALOUT) THEN IF (.NOT. ALLOCATED(WW3GLOBAL)) THEN ALLOCATE(WW3GLOBAL(8,MNP), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_bdcons, allocate error 7') END IF END IF IF (LBCWA) THEN ! Parametric Wave Boundary is prescribed WRITE(STAT%FHNDL,'("+TRACE...",A)') 'Parametric Wave Boundary Condition is prescribed' IF (LINHOM) THEN ! Inhomogenous in space IF (LBCSE) THEN ! Unsteady in time SPPARM = 0. IF (IBOUNDFORMAT == 1) THEN ! WWM CALL READWAVEPARWWM ELSE IF (IBOUNDFORMAT == 2) THEN ! FVCOM ... THIS WILL BE REPLACED BY SWAN TYPE BOUNDARY! CALL READWAVEPARFVCOM ELSE IF (IBOUNDFORMAT == 3) THEN ! WW3 #ifdef NCDF CALL READ_NETCDF_WW3_PARAM #else CALL WWM_ABORT('compile with netcdf for IBOUNDFORMAT=3') #endif CALL INTER_STRUCT_BOUNDARY(NDX_BND,NDY_BND,DX_BND,DY_BND,OFFSET_X_BND,OFFSET_Y_BND,SPPARM) IF (LWW3GLOBALOUT) CALL INTER_STRUCT_DOMAIN(NDX_BND,NDY_BND,DX_BND,DY_BND,OFFSET_X_BND,OFFSET_Y_BND,WW3GLOBAL) ELSE IF (IBOUNDFORMAT == 4) THEN ! WWM SPPARM netcdf file ! Print *, 'Before READ_NETCDF_BOUNDARY_SPPARM' #ifdef NCDF CALL READ_NETCDF_BOUNDARY_SPPARM #else CALL WWM_ABORT('compile with netcdf for IBOUNDFORMAT=4') #endif ELSE IF (IBOUNDFORMAT == 5) THEN ! WAM format of waves CALL WWM_ABORT('No possibility of using parametric boundary for IBOUNDFORMAT=5') END IF ELSE ! Steady ... SPPARM = 0. IF (IBOUNDFORMAT == 1) THEN CALL READWAVEPARWWM ELSE IF (IBOUNDFORMAT == 2) THEN CALL READWAVEPARFVCOM ELSE IF (IBOUNDFORMAT == 3) THEN #ifdef NCDF CALL READ_NETCDF_WW3_PARAM #else CALL WWM_ABORT('compile with netcdf for IBOUNDFORMAT=3') #endif CALL INTER_STRUCT_BOUNDARY(NDX_BND,NDY_BND,DX_BND,DY_BND,OFFSET_X_BND,OFFSET_Y_BND,SPPARM) END IF END IF ! LBCSE ... DO IP = 1, IWBMNP CALL SPECTRAL_SHAPE(SPPARM(:,IP),WBACOUT(:,:,IP),.FALSE.,'CALL FROM WB 1', .FALSE.) END DO ELSE ! Homogenous in space IF (IWBMNP .gt. 0) THEN IF (LBCSE) THEN ! Unsteady in time IF (IBOUNDFORMAT == 1) THEN CALL READWAVEPARWWM ELSE IF (IBOUNDFORMAT == 2) THEN CALL READWAVEPARFVCOM END IF CALL SPECTRAL_SHAPE(SPPARM(:,1),WBACOUT(:,:,1), .FALSE.,'CALL FROM WB 3', .FALSE.) ELSE ! Steady in time ... SPPARM = 0. IF (LMONO_IN) THEN SPPARM(1,1) = WBHS * SQRT(2.) ELSE SPPARM(1,1) = WBHS END IF SPPARM(2,1) = WBTP SPPARM(3,1) = WBDM SPPARM(4,1) = WBDS SPPARM(5,1) = WBSS SPPARM(6,1) = WBDSMS SPPARM(7,1) = WBGAUSS SPPARM(8,1) = WBPKEN CALL SPECTRAL_SHAPE(SPPARM(:,1),WBACOUT(:,:,1),.FALSE.,'CALL FROM WB 4', .TRUE.) END IF ! LBCSE END IF END IF ! LINHOM END IF IF (LBCSP) THEN ! Spectrum is prescribed IF (LINHOM) THEN ! The boundary conditions is not homogenous! IF (LBSP1D) CALL WWM_ABORT('No inhomogenous 1d spectra boundary cond. available') IF (IBOUNDFORMAT == 1) THEN ! WWM !CALL READSPEC2D CALL WWM_ABORT('No inhomogenous 2d spectra boundary cond. available in WWM Format') END IF IF (IBOUNDFORMAT == 3) THEN ! WW3 CALL GET_BINARY_WW3_SPECTRA(WBACOUT) END IF IF (IBOUNDFORMAT == 4) THEN ! WWM WBAC netcdf #ifdef NCDF CALL READ_NETCDF_BOUNDARY_WBAC(WBACOUT) #else CALL WWM_ABORT('compile with netcdf for IBOUNDFORMAT=4') #endif END IF IF (IBOUNDFORMAT == 5) THEN ! WWM WBAC netcdf Print *, 'Before call to READ_GRIB_WAM_BOUNDARY_WBAC' #ifdef GRIB_API_ECMWF CALL READ_GRIB_WAM_BOUNDARY_WBAC(WBACOUT) #else CALL WWM_ABORT('compile with GRIB_API for IBOUNDFORMAT=5') #endif END IF IF (IBOUNDFORMAT == 6) THEN ! WW3 netcdf #ifdef NCDF CALL GET_NC_WW3_SPECTRA(WBACOUT) #else CALL WWM_ABORT('compile with netcdf for IBOUNDFORMAT=4') #endif END IF ELSE ! The boundary conditions is homogeneous in space ! IF (LBSP1D) THEN ! 1-D Spectra is prescribed WRITE(STAT%FHNDL,'("+TRACE...",A)') '1d Spectra is given as Wave Boundary Condition' CALL READSPEC1D(LFIRSTREAD) CALL SPECTRUM_INT(WBACOUT(:,:,1)) ELSE IF (LBSP2D) THEN ! 2-D Spectra is prescribed WRITE(STAT%FHNDL,'("+TRACE...",A)') '2d Spectra is given as Wave Boundary Condition' IF (IBOUNDFORMAT == 1) THEN CALL READSPEC2D ELSE IF (IBOUNDFORMAT == 3) THEN CALL GET_BINARY_WW3_SPECTRA(WBACOUT) ELSE IF (IBOUNDFORMAT == 6) THEN CALL GET_NC_WW3_SPECTRA(WBACOUT) END IF CALL SPECTRUM_INT(WBACOUT(:,:,1)) END IF ! LBSP1D .OR. LBSP2D END IF ! LINHOM ENDIF ! LBCSP END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE SET_WAVE_BOUNDARY_CONDITION USE DATAPOOL IMPLICIT NONE CHARACTER(len=29) :: CHR IF (LNANINFCHK) THEN WRITE(DBG%FHNDL,*) ' ENTERING SET_WAVE_BOUNDARY_CONDITION ', SUM(AC2) IF (SUM(AC2) .NE. SUM(AC2)) CALL WWM_ABORT('NAN IN BOUNDARY CONDTITION l.1959') ENDIF IF ( MAIN%TMJD > SEBO%TMJD-1.E-8 .AND. MAIN%TMJD < SEBO%EMJD ) THEN ! Read next time step from boundary file ... IF (LBINTER) THEN CALL WAVE_BOUNDARY_CONDITION(WBACNEW) DSPEC = (WBACNEW-WBACOLD)/SEBO%DELT*MAIN%DELT WBAC = WBACOLD WBACOLD = WBACNEW ELSE ! .NOT. LBINTER CALL WAVE_BOUNDARY_CONDITION(WBAC) END IF ! LBINTER SEBO%TMJD = SEBO%TMJD + SEBO%DELT*SEC2DAY ELSE ! Interpolate in time ... no need to read ... IF (LBINTER) THEN WBAC = WBAC + DSPEC END IF END IF CALL SET_WAVE_BOUNDARY IF (LNANINFCHK) THEN WRITE(DBG%FHNDL,*) ' FINISHED WITH BOUNDARY CONDITION ', SUM(AC2) IF (SUM(AC2) .NE. SUM(AC2)) CALL WWM_ABORT('NAN IN BOUNDARY CONDTITION l.1959') ENDIF END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE INIT_WAVE_BOUNDARY_CONDITION USE DATAPOOL IMPLICIT NONE REAL(rkind) :: DTMP INTEGER :: ITMP LOGICAL :: DoAllocate WRITE(STAT%FHNDL, *) 'Begin of INIT_WAVE_BOUNDARY_CONDITION' !TODO: Makes sure initial condition work also when no wave boundary is set ... IF (IBOUNDFORMAT == 3) THEN IF (LBCSP) THEN ! Spectrum is prescribed CALL INIT_BINARY_WW3_SPECTRA END IF IF (LBCWA) THEN ! Parametric Wave Boundary is prescribed #ifdef NCDF CALL INIT_NETCDF_WW3_WAVEPARAMETER #else CALL WWM_ABORT('Compile with netcdf For WW3 bdcons') #endif END IF END IF IF (IBOUNDFORMAT == 4) THEN #ifdef NCDF CALL INIT_NETCDF_BOUNDARY_WWM #else CALL WWM_ABORT('Compile with netcdf for IBOUNDFORMAT=4') #endif END IF IF (IBOUNDFORMAT == 5) THEN #ifdef GRIB_API_ECMWF CALL INIT_GRIB_WAM_BOUNDARY #else CALL WWM_ABORT('Compile with GRIB_API for IBOUNDFORMAT=5') #endif END IF IF (IBOUNDFORMAT == 6) THEN #ifdef NCDF IF (LBCSP) THEN ! Spectrum is prescribed CALL INIT_NC_WW3_SPECTRA ELSE CALL WWM_ABORT('IBOUNDFORMAT == 6 is for prescribed spectra only') END IF #else CALL WWM_ABORT('Compile with netcdf for IBOUNDFORMAT=6') #endif END IF IF ((IBOUNDFORMAT .eq. 4).or.(IBOUNDFORMAT .eq. 5).or.LEXPORT_BOUC_WW3.or.BOUC_NETCDF_OUT_PARAM.or.BOUC_NETCDF_OUT_SPECTRA) THEN #ifdef MPI_PARALL_GRID CALL SETUP_BOUNDARY_SCATTER_REDUCE_ARRAY #endif #ifdef MPI_PARALL_GRID IF (myrank .eq. rank_boundary) THEN DoAllocate=.TRUE. ELSE DoAllocate=.FALSE. END IF #else DoAllocate=.TRUE. #endif WRITE(STAT%FHNDL, *) 'DoAllocate=', DoAllocate IF (DoAllocate) THEN IF (LBCWA .or. BOUC_NETCDF_OUT_PARAM) THEN IF (.NOT. ALLOCATED(SPPARM_GL)) THEN allocate(SPPARM_GL(8,IWBMNPGL), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_bdcons, allocate error 26') END IF END IF IF (LBCSP .or. LEXPORT_BOUC_WW3 .or. BOUC_NETCDF_OUT_SPECTRA) THEN IF (.NOT. ALLOCATED(WBAC_GL)) THEN allocate(WBAC_GL(MSC,MDC,IWBMNPGL), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_bdcons, allocate error 27') END IF END IF END IF END IF WRITE(STAT%FHNDL, *) 'LEXPORT_BOUC_WW3=', LEXPORT_BOUC_WW3 WRITE(STAT%FHNDL, *) 'IWBMNP=', IWBMNP WRITE(STAT%FHNDL, *) 'allocated(WBAC)=', allocated(WBAC) WRITE(STAT%FHNDL, *) 'allocated(WBAC_GL)=', allocated(WBAC_GL) IF (LBCWA .or. LBCSP) THEN CALL WAVE_BOUNDARY_CONDITION(WBAC) IF (LBINTER) WBACOLD = WBAC END IF END SUBROUTINE !********************************************************************** !* * !********************************************************************** #ifdef NCDF SUBROUTINE INIT_NETCDF_WW3_WAVEPARAMETER USE DATAPOOL USE NETCDF IMPLICIT NONE INTEGER :: IT, IFILE, IVAR, BND_NCID INTEGER :: ILON_ID, ILAT_ID, ITIME_ID, I, J, COUNTER REAL(rkind), ALLOCATABLE :: BND_TIME(:) character (len = *), parameter :: CallFct = "INIT_NETCDF_WW3_WAVEPARAMETER" integer, dimension(nf90_max_var_dims) :: dimIDs INTEGER MAX_NDT_BND_FILE, iProc, eSize INTEGER iVect(4) REAL(rkind) rVect(4) # ifdef MPI_PARALL_GRID IF (MULTIPLE_IN_BOUND .or. (myrank .eq. 0)) THEN # endif CALL TEST_FILE_EXIST_DIE("Missing WW3 boundary file : ", TRIM(WAV%FNAME)) OPEN(WAV%FHNDL,FILE=WAV%FNAME,STATUS='OLD') ! WRITE(STAT%FHNDL,*) WAV%FHNDL, WAV%FNAME, BND%FHNDL, BND%FNAME NUM_NETCDF_FILES_BND = 0 DO READ( WAV%FHNDL, *, IOSTAT = ISTAT ) IF ( ISTAT /= 0 ) EXIT NUM_NETCDF_FILES_BND = NUM_NETCDF_FILES_BND + 1 END DO REWIND(WAV%FHNDL) ! WRITE(STAT%FHNDL,*) 'NUM_NETCDF_FILES_BND=', NUM_NETCDF_FILES_BND NUM_NETCDF_FILES_BND = NUM_NETCDF_FILES_BND / NUM_NETCDF_VAR_TYPES WRITE(STAT%FHNDL,*) 'NUM_NETCDF_FILES_BND=', NUM_NETCDF_FILES_BND WRITE(STAT%FHNDL,*) 'NUM_NETCDF_VAR_TYPES=', NUM_NETCDF_VAR_TYPES ALLOCATE(NETCDF_FILE_NAMES_BND(NUM_NETCDF_FILES_BND,NUM_NETCDF_VAR_TYPES), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_bdcons, allocate error 9') DO IT = 1, NUM_NETCDF_FILES_BND DO IVAR = 1, NUM_NETCDF_VAR_TYPES READ( WAV%FHNDL, *) NETCDF_FILE_NAMES_BND(IT,IVAR) END DO END DO CLOSE (WAV%FHNDL) ! ! four files are read to set up the wave spectra Hs, Tm01, Dir, Sprd ! ALLOCATE(NDT_BND_FILE(NUM_NETCDF_FILES_BND), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_bdcons, allocate error 10') NDT_BND_FILE = 0 ! DO IFILE = 1, NUM_NETCDF_FILES_BND ! WRITE(STAT%FHNDL,'(I10,10X,5A30)') IFILE, NETCDF_FILE_NAMES_BND(IFILE,:) ! END DO ! ! check number of time steps in netcdf file ... it is assumed that all files have the same ammount of time steps ... ! DO IFILE = 1, NUM_NETCDF_FILES_BND WRITE(STAT%FHNDL,*) ifile, TRIM(NETCDF_FILE_NAMES_BND(IFILE,1)) FLUSH(STAT%FHNDL) CALL TEST_FILE_EXIST_DIE("Missing ww3 boundary condition file : ", TRIM(NETCDF_FILE_NAMES_BND(IFILE,1))) ISTAT = NF90_OPEN(TRIM(NETCDF_FILE_NAMES_BND(IFILE,1)), NF90_NOWRITE, BND_NCID) IF (ISTAT /= 0) THEN Print *, 'Error while trying to open netcdf file' Print *, 'FILE=', TRIM(NETCDF_FILE_NAMES_BND(IFILE,1)) Print *, 'One possible error is that the file is NC4 but you linked to NC3' CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 1, ISTAT) END IF ISTAT = nf90_inq_varid(BND_NCID, 'time', ITIME_ID) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 2, ISTAT) ISTAT = NF90_INQUIRE_VARIABLE(BND_NCID, ITIME_ID, dimids = dimids) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 3, ISTAT) ISTAT = nf90_inquire_dimension(BND_NCID, dimIDs(1), len = NDT_BND_FILE(IFILE)) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 4, ISTAT) ! WRITE(STAT%FHNDL,*) IFILE, NDT_BND_FILE(IFILE) END DO ! ! check dimensions in the netcdf ... again it is assumed that this is not changing for all files ... ! ISTAT = nf90_inq_varid(BND_NCID, 'longitude', ILON_ID) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 5, ISTAT) ISTAT = NF90_INQUIRE_VARIABLE(BND_NCID, ILON_ID, dimids = dimIDs) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 6, ISTAT) ISTAT = nf90_inquire_dimension(BND_NCID, dimIDs(1), len = NDX_BND) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 7, ISTAT) ISTAT = nf90_inq_varid(BND_NCID, 'latitude', ILAT_ID) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 8, ISTAT) ISTAT = NF90_INQUIRE_VARIABLE(BND_NCID, ILAT_ID, dimids = dimIDs) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 9, ISTAT) ISTAT = nf90_inquire_dimension(BND_NCID, dimIDs(1), len = NDY_BND) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 10, ISTAT) ! WRITE(STAT%FHNDL,*) 'Number of Gridpoints', NDX_BND, NDY_BND ALLOCATE (COORD_BND_X(NDX_BND), COORD_BND_Y(NDY_BND), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_bdcons, allocate error 11') ! ! read coordinates from files .... ! ISTAT = NF90_GET_VAR(BND_NCID, ILON_ID, COORD_BND_X) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 11, ISTAT) ISTAT = NF90_GET_VAR(BND_NCID, ILAT_ID, COORD_BND_Y) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 12, ISTAT) ! ! estimate offset ... ! OFFSET_X_BND = MINVAL(COORD_BND_X) OFFSET_Y_BND = MINVAL(COORD_BND_Y) ! ! resolution ... ! DX_BND = ABS(MAXVAL(COORD_BND_X) - MINVAL(COORD_BND_X))/(NDX_BND-1) DY_BND = ABS(MAXVAL(COORD_BND_Y) - MINVAL(COORD_BND_Y))/(NDY_BND-1) ! ! close netcdf file ... ! ISTAT = NF90_CLOSE(BND_NCID) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 13, ISTAT) ! ! total number of time steps ... in all files ! NDT_BND_ALL_FILES = 0 ! write(STAT%FHNDL,*) NUM_NETCDF_FILES_BND DO IT = 1, NUM_NETCDF_FILES_BND NDT_BND_ALL_FILES = NDT_BND_ALL_FILES + NDT_BND_FILE(IT) ! write(STAT%FHNDL,*) it, NDT_BND_FILE(it) END DO ! WRITE(STAT%FHNDL,*) NDT_BND_ALL_FILES, NDT_BND_FILE WRITE(STAT%FHNDL, *) 'INIT_NETCDF_WW3_WAVEPARAMETER' MAX_NDT_BND_FILE=MAXVAL(NDT_BND_FILE) ALLOCATE (BND_TIME_ALL_FILES(NUM_NETCDF_FILES_BND,MAX_NDT_BND_FILE), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_bdcons, allocate error 12') BND_TIME_ALL_FILES = ZERO ! ! read all time steps in the proper format and transform in wwm time line ! DO IFILE = 1, NUM_NETCDF_FILES_BND CALL TEST_FILE_EXIST_DIE("Missing ww3 boundary condition file : ", TRIM(NETCDF_FILE_NAMES_BND(IFILE,1))) ISTAT = NF90_OPEN(NETCDF_FILE_NAMES_BND(IFILE,1),NF90_NOWRITE,BND_NCID) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 14, ISTAT) ALLOCATE (BND_TIME(NDT_BND_FILE(IFILE)), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_bdcons, allocate error 13') BND_TIME = ZERO ! MDS: It looks dangerous to use previous id. ISTAT = NF90_GET_VAR(BND_NCID,ITIME_ID,BND_TIME) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 15, ISTAT) DO IT = 1, NDT_BND_FILE(IFILE) BND_TIME_ALL_FILES(IFILE,IT) = BND_TIME(IT) ! CALL CT2MJD('19000101.000000',DTMP1) ! CALL CT2MJD('19900101.000000',DTMP2) ! CALL MJD2CT(DTMP1,chrdate) ! WRITE(*,*) '19000101.000000', DTMP1, chrdate ! CALL MJD2CT(DTMP2,chrdate) ! WRITE(*,*) '19900101.000000', DTMP2, chrdate ! CALL MJD2CT(0.0_rkind,chrdate) ! WRITE(*,*) '00000000.000000', 0.0_rkind, chrdate ! WRITE(*,*) BND_TIME_ALL_FILES(1,1), DT_DIFF_19901900 ! IF (IT == 1 .AND. IFILE ==1) WRITE(*,*) DTMP1, DTMP2, DTMP1+DT_DIFF_19901900 ! IF (IT == 1 .AND. IFILE ==1) WRITE(*,*) IFILE, IT, BND_TIME(IT), chrdate END DO DEALLOCATE(BND_TIME) END DO BND_TIME_ALL_FILES = BND_TIME_ALL_FILES + DT_DIFF_19901900 IF (LWRITE_ALL_WW3_RESULTS) THEN OPEN(3010, FILE = 'sysglobalboundary.dat', STATUS = 'UNKNOWN') WRITE (3010, '(I10)') 0 WRITE (3010, '(I10)') NDX_BND * NDY_BND COUNTER = 0 DO I = 1, NDY_BND DO J = 1, NDX_BND WRITE (3010, '(I10,3F15.4)') COUNTER, OFFSET_X_BND+(J-1)*DX_BND,OFFSET_Y_BND+(I-1)*DY_BND, 0.0 COUNTER = COUNTER + 1 END DO END DO WRITE (3010, *) (NDX_BND-1)*(NDY_BND-1)*2 DO J = 0, NDY_BND-2 DO I = 0, NDX_BND-2 WRITE (3010, '(5I10)') I+J*NDX_BND , NDX_BND+I+J* NDX_BND, NDX_BND+I+1+J*NDX_BND, 0, 0 WRITE (3010, '(5I10)') NDX_BND+I+1+J*NDX_BND, I+1+J*NDX_BND , I+J*NDX_BND , 0, 0 END DO END DO CLOSE(3010) END IF # ifdef MPI_PARALL_GRID END IF # endif # ifdef MPI_PARALL_GRID IF (.NOT. MULTIPLE_IN_BOUND) THEN IF (myrank.eq.0) THEN iVect(1)=NUM_NETCDF_FILES_BND iVect(2)=MAX_NDT_BND_FILE iVect(3)=NDX_BND iVect(4)=NDY_BND DO IPROC=2,nproc CALL MPI_SEND(iVect,4,itype, iProc-1, 811, comm, ierr) END DO rVect(1)=OFFSET_X_BND rVect(2)=OFFSET_Y_BND rVect(3)=DX_BND rVect(4)=DY_BND DO IPROC=2,nproc CALL MPI_SEND(rVect,4,rtype, iProc-1, 820, comm, ierr) END DO eSize=NUM_NETCDF_FILES_BND*MAX_NDT_BND_FILE DO IPROC=2,nproc CALL MPI_SEND(NDT_BND_FILE,NUM_NETCDF_FILES_BND,itype, iProc-1, 812, comm, ierr) CALL MPI_SEND(BND_TIME_ALL_FILES,eSize,rtype, iProc-1, 813, comm, ierr) END DO ELSE CALL MPI_RECV(iVect,4,itype, 0, 811, comm, istatus, ierr) NUM_NETCDF_FILES_BND=iVect(1) MAX_NDT_BND_FILE=iVect(2) NDX_BND=iVect(3) NDY_BND=iVect(4) CALL MPI_RECV(rVect,4,rtype, 0, 820, comm, istatus, ierr) OFFSET_X_BND=rVect(1) OFFSET_Y_BND=rVect(2) DX_BND=rVect(3) DY_BND=rVect(4) ALLOCATE(BND_TIME_ALL_FILES(NUM_NETCDF_FILES_BND,MAX_NDT_BND_FILE), NDT_BND_FILE(NUM_NETCDF_FILES_BND), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_bdcons, allocate error 14') CALL MPI_RECV(NDT_BND_FILE,NUM_NETCDF_FILES_BND,itype, 0, 812, comm, istatus, ierr) eSize=NUM_NETCDF_FILES_BND*MAX_NDT_BND_FILE CALL MPI_RECV(BND_TIME_ALL_FILES,eSize,rtype, 0, 813, comm, istatus, ierr) END IF END IF # endif SEBO%DELT = (BND_TIME_ALL_FILES(1,2) - BND_TIME_ALL_FILES(1,1)) * DAY2SEC ! write(STAT%FHNDL,*) SEBO%DELT, BND_TIME_ALL_FILES(1,2), BND_TIME_ALL_FILES(1,1) ALLOCATE (HS_WW3(NDX_BND,NDY_BND), FP_WW3(NDX_BND,NDY_BND), T02_WW3(NDX_BND,NDY_BND), DSPR_WW3(NDX_BND,NDY_BND), DIR_WW3(NDX_BND,NDY_BND), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_bdcons, allocate error 15') HS_WW3 = 0. FP_WW3 = 0. T02_WW3 = 0. DSPR_WW3 = 0. DIR_WW3 = 0. END SUBROUTINE #endif !********************************************************************** !* * !********************************************************************** #ifdef NCDF SUBROUTINE READ_NETCDF_WW3_IVAR(IFILE, IT, IVAR, VAR_READ) USE DATAPOOL USE NETCDF IMPLICIT NONE integer, intent(in) :: IFILE, IT, IVAR character (len = *), parameter :: CallFct = "READ_NETCDF_WW3_IVAR" CHARACTER(LEN=40) :: EVAR REAL(rkind), intent(out) :: VAR_READ(NDX_BND,NDY_BND) integer :: ITMP(NDX_BND,NDY_BND) REAL(rkind) :: scale_factor integer ncid, var_id IF (IVAR .eq. 3) THEN EVAR='hs' ELSE IF (IVAR .eq. 2) THEN EVAR='fp' ELSE IF (IVAR .eq. 5) THEN EVAR='t02' ELSE IF (IVAR .eq. 4) THEN EVAR='spr' ELSE IF (IVAR .eq. 1) THEN EVAR='dir' ELSE Print *, 'IVAR=', IVAR CALL WWM_ABORT('Wrong IVAR') END IF CALL TEST_FILE_EXIST_DIE("Missing ww3 boundary condition file : ", TRIM(NETCDF_FILE_NAMES_BND(IFILE,IVAR))) ISTAT = NF90_OPEN(NETCDF_FILE_NAMES_BND(IFILE,IVAR),NF90_NOWRITE,ncid) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 1, ISTAT) ISTAT = nf90_inq_varid(ncid, TRIM(EVAR), var_id) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 2, ISTAT) ISTAT = nf90_get_att(ncid, var_id, 'scale_factor', scale_factor) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 3, ISTAT) ISTAT = NF90_GET_VAR(ncid, var_id, ITMP, start = (/ 1, 1, IT /), count = (/ NDX_BND, NDY_BND, 1/)) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 4, ISTAT) VAR_READ = MyREAL(ITMP) * scale_factor ISTAT = nf90_close(ncid) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 5, ISTAT) END SUBROUTINE #endif !********************************************************************** !* * !********************************************************************** #ifdef NCDF SUBROUTINE READ_NETCDF_WW3_SINGLE(IFILE,IT) USE DATAPOOL IMPLICIT NONE INTEGER, INTENT(IN) :: IFILE, IT CALL READ_NETCDF_WW3_IVAR(IFILE, IT, 3, HS_WW3) CALL READ_NETCDF_WW3_IVAR(IFILE, IT, 2, FP_WW3) CALL READ_NETCDF_WW3_IVAR(IFILE, IT, 5, T02_WW3) CALL READ_NETCDF_WW3_IVAR(IFILE, IT, 4, DSPR_WW3) CALL READ_NETCDF_WW3_IVAR(IFILE, IT, 1, DIR_WW3) END SUBROUTINE #endif !********************************************************************** !* * !********************************************************************** #ifdef NCDF SUBROUTINE READ_NETCDF_WW3_PARAM USE DATAPOOL IMPLICIT NONE INTEGER :: IFILE, IT INTEGER :: counter, ip, i, j REAL(rkind), ALLOCATABLE :: U(:), V(:), H(:) REAL(rkind), SAVE :: TIME, scale_factor INTEGER IX, IY, IPROC REAL(rkind), ALLOCATABLE :: ARR_send_recv(:) integer, allocatable :: bnd_rqst(:), bnd_stat(:,:) CALL COMPUTE_IFILE_IT(IFILE, IT) # ifdef MPI_PARALL_GRID IF (MULTIPLE_IN_BOUND) THEN CALL READ_NETCDF_WW3_SINGLE(IFILE,IT) ELSE allocate(ARR_send_recv(5*NDX_BND*NDY_BND)) IF (istat/=0) CALL WWM_ABORT('wwm_bdcons, allocate error 16') IF (myrank .eq. 0) THEN CALL READ_NETCDF_WW3_SINGLE(IFILE,IT) J=0 DO IX=1,NDX_BND DO IY=1,NDY_BND J=J+1 ARR_send_recv(J)=HS_WW3(IX,IY) J=J+1 ARR_send_recv(J)=FP_WW3(IX,IY) J=J+1 ARR_send_recv(J)=T02_WW3(IX,IY) J=J+1 ARR_send_recv(J)=DSPR_WW3(IX,IY) J=J+1 ARR_send_recv(J)=DIR_WW3(IX,IY) END DO END DO allocate(bnd_rqst(nproc-1), bnd_stat(MPI_STATUS_SIZE,nproc-1), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_bdcons, allocate error 17') DO iProc=2,nproc CALL mpi_isend(ARR_send_recv, 5*NDX_BND*NDY_BND, rtype, iProc-1, 2032, comm, bnd_rqst(iProc-1), ierr) END DO IF (nproc > 1) THEN CALL MPI_WAITALL(nproc-1, bnd_rqst, bnd_stat, ierr) END IF deallocate(bnd_rqst, bnd_stat) ELSE CALL MPI_RECV(ARR_send_recv, 5*NDX_BND*NDY_BND, rtype, 0, 2032, comm, istatus, ierr) J=0 DO IX=1,NDX_BND DO IY=1,NDY_BND J=J+1 HS_WW3(IX,IY) = ARR_send_recv(J) J=J+1 FP_WW3(IX,IY) = ARR_send_recv(J) J=J+1 T02_WW3(IX,IY) = ARR_send_recv(J) J=J+1 DSPR_WW3(IX,IY) = ARR_send_recv(J) J=J+1 DIR_WW3(IX,IY) = ARR_send_recv(J) END DO END DO END IF END IF # else CALL READ_NETCDF_WW3_SINGLE(IFILE,IT) # endif IF (LWRITE_WW3_RESULTS) THEN OPEN(3012, FILE = 'ergwiii.bin', FORM = 'UNFORMATTED') ALLOCATE(U(NDX_BND*NDY_BND), V(NDX_BND*NDY_BND), H(NDX_BND*NDY_BND), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_bdcons, allocate error 18') COUNTER = 1 DO J = 1, NDY_BND DO I = 1, NDX_BND U(COUNTER) = HS_WW3(I,J) V(COUNTER) = DIR_WW3(I,J) H(COUNTER) = DSPR_WW3(I,J) COUNTER = COUNTER + 1 END DO END DO TIME = TIME + 1. WRITE(3012) TIME WRITE(3012) (U(IP), V(IP), H(IP), IP = 1, NDX_BND*NDY_BND) DEALLOCATE(U,V,H) END IF END SUBROUTINE #endif !********************************************************************** !* * !********************************************************************** SUBROUTINE SPPARM_INTER_STRUCT(NDX,NDY,DX,DY,OFFSET_X,OFFSET_Y, MNPT, XPT, YPT, VAL, DOPEAK) USE DATAPOOL IMPLICIT NONE INTEGER, INTENT(IN) :: NDX, NDY REAL(rkind), INTENT(IN) :: DX, DY, OFFSET_X, OFFSET_Y INTEGER, INTENT(IN) :: MNPT REAL(rkind), INTENT(IN) :: XPT(MNPT), YPT(MNPT) !JZ: error REAL(rkind), INTENT(OUT) :: VAL(8,MNPT) !IWBMNP) LOGICAL, INTENT(IN) :: DOPEAK REAL(rkind) :: eVect(5) REAL(rkind) :: WX, WX1, WX2, WX3, WX4 REAL(rkind) :: eVAL, HX1, HX2, LEN_X, LEN_Y REAL(rkind) :: DELTA_X, DELTA_Y INTEGER :: IVAR, I, J, IDX1, IDX2 INTEGER :: IP, J_INT, I_INT DO IP = 1, MNPT LEN_X = XPT(IP) - OFFSET_X LEN_Y = YPT(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 ! DO IVAR=1,5 DO I=0,1 DO J=0,1 IDX1=I_INT + I IDX2=J_INT + J IF (IVAR .eq. 1) THEN WX=HS_WW3(IDX1,IDX2) ELSE IF (IVAR .eq. 2) THEN WX=DIR_WW3(IDX1,IDX2) ELSE IF (IVAR .eq. 3) THEN WX=FP_WW3(IDX1,IDX2) ELSE IF (IVAR .eq. 4) THEN WX=T02_WW3(IDX1,IDX2) ELSE IF (IVAR .eq. 5) THEN WX=DSPR_WW3(IDX1,IDX2) ELSE CALL WWM_ABORT('Wrong IVAR') END IF IF ((I .eq. 0).and.(J .eq. 0)) THEN WX1=WX END IF IF ((I .eq. 0).and.(J .eq. 1)) THEN WX2=WX END IF IF ((I .eq. 1).and.(J .eq. 1)) THEN WX3=WX END IF IF ((I .eq. 1).and.(J .eq. 0)) THEN WX4=WX END IF END DO END DO HX1 = WX1 + (WX4-WX1)/DX * DELTA_X HX2 = WX2 + (WX3-WX2)/DX * DELTA_X IF (WX1 .LT. 0. .OR. WX2 .LT. 0. .OR. WX3 .LT. 0. .OR. WX4 .LT. 0. ) THEN eVAL = 0. ELSE eVAL = HX1 + (HX2-HX1)/DY * DELTA_Y ENDIF eVect(IVAR)=eVAL END DO VAL(1,IP)=eVect(1) IF (DOPEAK) THEN IF (eVect(3) .gt. TINY(1.)) THEN VAL(2,IP)=ONE/eVect(3) VAL(5,IP) = 2. END IF ELSE VAL(2,IP)=eVect(4) VAL(5,IP) = -2. END IF VAL(3,IP) = eVect(2) VAL(4,IP) = eVect(5) VAL(6,IP) = 1. VAL(7,IP) = 0.1 VAL(8,IP) = 3.3 ! SPPARM(1): Hs, sign. wave height ! SPPARM(2): Wave period given by user (either peak or mean) ! SPPARM(3): average direction ! SPPARM(4): directional spread ! SPPARM(5): spectral shape (1-4), ! 1 - Pierson-Moskowitz ! 2 - JONSWAP ! 3 - BIN ! 4 - Gauss ! positive peak (+) or mean frequency (-) ! SPPARM(6): directional spreading in degree (1) or exponent (2) ! SPPARM(7): gaussian width for the gauss spectrum 0.1 ! SPPARM(8): peak enhancement factor for the JONSWAP spectra 3. END DO END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE INTER_STRUCT_BOUNDARY(NDX,NDY,DX,DY,OFFSET_X,OFFSET_Y,VAL) USE DATAPOOL IMPLICIT NONE INTEGER, INTENT(IN) :: NDX, NDY REAL(rkind), INTENT(IN) :: DX, DY, OFFSET_X, OFFSET_Y REAL(rkind), INTENT(OUT) :: VAL(8,IWBMNP) REAL(rkind) :: XPT(IWBMNP), YPT(IWBMNP) INTEGER :: IP DO IP=1,IWBMNP XPT(IP)=XP(IWBNDLC(IP)) YPT(IP)=YP(IWBNDLC(IP)) END DO CALL SPPARM_INTER_STRUCT(NDX,NDY,DX,DY,OFFSET_X,OFFSET_Y, IWBMNP, XPT, YPT, VAL, DOPEAK_BOUNDARY) END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE INTER_STRUCT_DOMAIN(NDX,NDY,DX,DY,OFFSET_X,OFFSET_Y,VAL) USE DATAPOOL IMPLICIT NONE INTEGER, INTENT(IN) :: NDX, NDY REAL(rkind), INTENT(IN) :: DX, DY, OFFSET_X, OFFSET_Y !JZ: error REAL(rkind), INTENT(OUT) :: VAL(8,MNP) !IWBMNP) CALL SPPARM_INTER_STRUCT(NDX,NDY,DX,DY,OFFSET_X,OFFSET_Y, MNP, XP, YP, VAL, DOPEAK_BOUNDARY) END SUBROUTINE !********************************************************************** !* * !********************************************************************** !AR: check this for wave boundary nodes .... SUBROUTINE READWAVEPARFVCOM USE DATAPOOL IMPLICIT NONE INTEGER :: IP # ifdef MPI_PARALL_GRID REAL(rkind) :: RTMP INTEGER :: IPP # endif logical :: LOPEN2 ! SPPARM(1): Hs, sign. wave height ! SPPARM(2): Wave period given by user (either peak or mean) ! SPPARM(3): average direction ! SPPARM(4): directional spread ! SPPARM(5): spectral shape (1-4), ! (1 - Pierson-Moskowitz, ! 2 - JONSWAP, ! 3 - BIN, ! 4 - Gauss) ! negative peak (+) ! or mean frequency (-) ! SPPARM(6): directional spreading in degree (1) or exponent (2) ! SPPARM(7): gaussian width for the gauss spectrum 0.1 ! SPPARM(8): peak enhancement factor for the JONSWAP spectra 3.3 SPPARM(4,:) = 20. SPPARM(5,:) = 2. SPPARM(6,:) = 2. SPPARM(7,:) = 0.1 SPPARM(8,:) = 3.3 IF (LINHOM) THEN inquire(WAV%FHNDL,opened=LOPEN2) if(.not.LOPEN2) open(WAV%FHNDL,file=WAV%FNAME,status='old') !write(12,*)'WWM:',WAV%FHNDL,WAV%FNAME READ(WAV%FHNDL,*) END IF # ifdef MPI_PARALL_GRID IPP = 0 IF (LINHOM) THEN DO IP = 1, IWBMNPGL IF(ipgl(IWBNDGL(IP))%rank == myrank) THEN ! IF boundary nodes belong to local domain read values into boundary array IPP = IPP + 1 READ (WAV%FHNDL, *) SPPARM(1,IPP), SPPARM(2,IPP), SPPARM(3,IPP) ELSE READ (WAV%FHNDL, *) RTMP, RTMP, RTMP ! ELSE ... throw them away ... ENDIF END DO ELSE READ (WAV%FHNDL, *) SPPARM(1,1), SPPARM(2,1), SPPARM(3,1) DO IP = 1, IWBMNPGL SPPARM(1:3,IP) = SPPARM(1:3,1) END DO END IF # else IF (LINHOM) THEN DO IP = 1, IWBMNP READ (WAV%FHNDL, *) SPPARM(1,IP), SPPARM(2,IP), SPPARM(3,IP) END DO ELSE READ (WAV%FHNDL, *) SPPARM(1,1), SPPARM(2,1), SPPARM(3,1) DO IP = 1, IWBMNP SPPARM(1:3,IP) = SPPARM(1:3,1) END DO END IF # endif END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE READSPEC1D(LFIRST) USE DATAPOOL IMPLICIT NONE INTEGER :: IP, IS LOGICAL, INTENT(IN) :: LFIRST ! ! Read Spectrum 1-D File ... ! Second Line ... number of frequencies ! Third Line ... number of directions ! IF (LFIRST) THEN CALL TEST_FILE_EXIST_DIE('Missing wave file : ', TRIM(WAV%FNAME)) OPEN(WAV%FHNDL, FILE = TRIM(WAV%FNAME), STATUS = 'OLD') READ (WAV%FHNDL,*) WBMSC READ (WAV%FHNDL,*) WBMDC !WRITE(*,*) WBMSC, WBMDC CALL ALLOC_SPEC_BND DO IS = 1, WBMSC READ (WAV%FHNDL,*) SFRQ(IS,1) END DO END IF IF (LINHOM) THEN DO IP = 1, IWBMNP DO IS = 1, WBMSC READ (WAV%FHNDL, *) SPEG(IS,1,IP), SDIR(IS,1), SPRD(IS,1) END DO END DO ELSE DO IS = 1, WBMSC READ (WAV%FHNDL, *) SPEG(IS,1,1), SDIR(IS,1), SPRD(IS,1) !WRITE(*,*) SPEG(IS,1,1), SDIR(IS,1), SPRD(IS,1) END DO END IF END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE READSPEC2D USE DATAPOOL IMPLICIT NONE END SUBROUTINE !********************************************************************** !* READSPEC2D_WW3INIT1 !* Read the header of a WAVEWATCHIII binary spectral file to get !* dimensions of the spectral grid and number of output locations !* required for dynamic allocation. !* !* Called by GETWW3SPECTRA !* !* Authors: Aron Roland !* Kai Li !* Guillaume Dodet (18/12/2012) !********************************************************************** SUBROUTINE READSPEC2D_BINARY_WW3_INIT_SPEC USE DATAPOOL IMPLICIT NONE CHARACTER(LEN=30) :: GNAME CHARACTER(LEN=21) :: LABEL CALL TEST_FILE_EXIST_DIE('Missing wave file : ', TRIM(WAV%FNAME)) OPEN(WAV%FHNDL,FILE=WAV%FNAME, STATUS='OLD',CONVERT='BIG_ENDIAN',FORM='UNFORMATTED') READ(WAV%FHNDL)LABEL, MSC_WW3,MDC_WW3, NP_WW3, GNAME WRITE(STAT%FHNDL,*)'START READSPEC2D_WW3_INIT_SPEC' WRITE(STAT%FHNDL,*)'LABEL, MSC_WW3,MDC_WW3, NP_WW3, GNAME' WRITE(STAT%FHNDL,*)LABEL, MSC_WW3,MDC_WW3, NP_WW3, GNAME CLOSE(WAV%FHNDL) WRITE(STAT%FHNDL,*)'DIRECTION NUMBER IN WW3 SPECTRUM:',MDC_WW3 WRITE(STAT%FHNDL,*)'FREQUENCY NUMBER IN WW3 SPECTRUM:',MSC_WW3 WRITE(STAT%FHNDL,'("+TRACE...",A)')'DONE READSPEC2D_WW3_INIT_SPEC' END SUBROUTINE !********************************************************************** !* READSPEC2D_WW3INIT2 !* Read the header of a WAVEWATCHIII binary spectral file to get !* frequencies and directions of the spectral grid and starting time !* !* Called by GETWW3SPECTRA !* !* Authors: Aron Roland !* Kai Li !* Guillaume Dodet (18/12/2012) !********************************************************************** SUBROUTINE READSPEC2D_BINARY_WW3_INIT_TIME USE DATAPOOL IMPLICIT NONE REAL :: SPECDMP(MSC_WW3,MDC_WW3) INTEGER :: TMP, IFLAG, IP, TIME(2), IT INTEGER, ALLOCATABLE :: ITIME(:,:) CHARACTER(LEN=30) :: GNAME CHARACTER(LEN=21) :: LABEL CHARACTER(LEN=10) :: PID CHARACTER(LEN=20) :: CTIME1, CTIME2 CHARACTER(LEN=15) :: TIMESTRING REAL :: TMP1(MSC_WW3),TMP2(MDC_WW3) !GD: in ww3 binary file, reals REAL :: TMPR1, TMPR2, TMPR3, TMPR4, TMPR5, TMPR6, TMPR7 WRITE(STAT%FHNDL,*)'START READSPEC2D_BINARY_WW3_INIT_TIME' ALLOCATE(FQ_WW3(MSC_WW3), DR_WW3(MDC_WW3), XP_WW3(NP_WW3), YP_WW3(NP_WW3), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_bdcons, allocate error 19') FQ_WW3=ZERO DR_WW3=ZERO XP_WW3=ZERO YP_WW3=ZERO CALL TEST_FILE_EXIST_DIE('Missing wave file : ', TRIM(WAV%FNAME)) OPEN(WAV%FHNDL, FILE = WAV%FNAME, STATUS = 'OLD',CONVERT='BIG_ENDIAN',FORM='UNFORMATTED') READ(WAV%FHNDL)LABEL,TMP,TMP,TMP,GNAME READ(WAV%FHNDL)TMP1 FQ_WW3 = TMP1 READ(WAV%FHNDL)TMP2 DR_WW3 = TMP2 MAXSTEP_WW3 = 0 DO READ(WAV%FHNDL,IOSTAT=IFLAG) TIME IF (IFLAG .GT. 0) THEN CALL WWM_ABORT('IFLAG incorrectly set 1') ELSE IF (IFLAG .LT. 0) THEN WRITE(STAT%FHNDL,*) 'END OF FILE REACHED AT 1, WHICH IS NICE' EXIT END IF DO IP = 1, NP_WW3 READ(WAV%FHNDL,IOSTAT=IFLAG) PID,TMPR1,TMPR2,TMPR3,TMPR4,TMPR5,TMPR6,TMPR7 ! WRITE(STAT%FHNDL,'(A10,7F15.4)') PID,TMPR1,TMPR2,TMPR3,TMPR4,TMPR5,TMPR6,TMPR7 IF (IFLAG .GT. 0) THEN CALL WWM_ABORT('IFLAG incorrectly set 2') ELSE IF (IFLAG .LT. 0) THEN WRITE(STAT%FHNDL,*) 'END OF FILE REACHED AT 2, WHICH IS NOT GOOD' EXIT END IF READ(WAV%FHNDL,IOSTAT=IFLAG) SPECDMP(:,:) IF (IFLAG .GT. 0) THEN CALL WWM_ABORT('IFLAG incorrectly set 3') ELSE IF (IFLAG .LT. 0) THEN WRITE(STAT%FHNDL,*) 'END OF FILE REACHED AT 3, WHICH IS NOT GOOD' EXIT END IF ENDDO MAXSTEP_WW3 = MAXSTEP_WW3 + 1 IF (MAXSTEP_WW3 == 1) TSTART_WW3 = TIME END DO WRITE(STAT%FHNDL,*) 'NUMBER OF BUOYS', NP_WW3 WRITE(STAT%FHNDL,*) 'NUMBER OF TIME STEPS IN FILE', MAXSTEP_WW3 REWIND(WAV%FHNDL) ALLOCATE(ITIME(MAXSTEP_WW3,2), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_bdcons, allocate error 20') ITIME = 0 READ(WAV%FHNDL)LABEL,TMP,TMP,TMP,GNAME READ(WAV%FHNDL)TMP1 READ(WAV%FHNDL)TMP2 DO IT = 1, MAXSTEP_WW3 READ(WAV%FHNDL,IOSTAT=IFLAG) ITIME(IT,:) DO IP = 1, NP_WW3 READ(WAV%FHNDL,IOSTAT=IFLAG) PID,TMPR1,TMPR2,TMPR3,TMPR4,TMPR5,TMPR6,TMPR7 READ(WAV%FHNDL,IOSTAT=IFLAG) SPECDMP(:,:) ENDDO END DO ALLOCATE (BND_TIME_ALL_FILES(1,MAXSTEP_WW3), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_bdcons, allocate error 21') BND_TIME_ALL_FILES = ZERO DO IT = 1, MAXSTEP_WW3 WRITE(CTIME1,*) ITIME(IT,1) CALL LEADINGZERO(ITIME(IT,2),CTIME2) TIMESTRING = ADJUSTL(TRIM(CTIME1)//'.'//ADJUSTL(CTIME2)) CALL CT2MJD(TIMESTRING, BND_TIME_ALL_FILES(1,IT)) ! WRITE(STAT%FHNDL,*) IT, TIMESTRING, BND_TIME_ALL_FILES(1,IT) END DO DTBOUND_WW3 = (BND_TIME_ALL_FILES(1,2)-BND_TIME_ALL_FILES(1,1))*DAY2SEC SEBO%DELT = DTBOUND_WW3 SEBO%BMJD = BND_TIME_ALL_FILES(1,1) SEBO%EMJD = BND_TIME_ALL_FILES(1,MAXSTEP_WW3) CLOSE(WAV%FHNDL) WRITE(STAT%FHNDL,*)'MIN. FREQ. IN WW3 SPECTRUM:',FQ_WW3(1) WRITE(STAT%FHNDL,*)'MAX. FREQ. IN WW3 SPECTRUM:',FQ_WW3(MSC_WW3) WRITE(STAT%FHNDL,*)'NUMBER OF TIME STEPS',MAXSTEP_WW3 WRITE(STAT%FHNDL,*)'TIME INCREMENT IN SPECTRAL FILE', DTBOUND_WW3 WRITE(STAT%FHNDL,*)'FIRST TIME STEP IN WW3 SPECTRUM FILE:',BND_TIME_ALL_FILES(1,1) WRITE(STAT%FHNDL,*)'BEGING TIME, END TIME and DELT of wave boundary', SEBO%BMJD, SEBO%EMJD, SEBO%DELT WRITE(STAT%FHNDL,*)'BEGING TIME, END TIME and DELT of simulation', MAIN%BMJD, MAIN%EMJD, MAIN%DELT WRITE(STAT%FHNDL,'("+TRACE...",A)') 'DONE READSPEC2D_WW3INIT2' END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE LEADINGZERO(INTIN,STRING) ! AR: some fortran magic with leading zero's ... INTEGER, INTENT(IN) :: INTIN CHARACTER(LEN=6), INTENT(OUT) :: STRING write( string, '(I6)' ) INTIN string = repeat('0',6-len_trim(adjustl(string)))//adjustl(string) END SUBROUTINE !********************************************************************** !* READSPEC2D_WW3 !* Reads spectra in WAVEWATCHIII binary spectral file !* !* Called by GETWW3SPECTRA !* !* Authors: Aron Roland !* Kai Li !* Guillaume Dodet (18/12/2012) !* !* Remarks: GD: Need to be modified for time interpolation. Input !* arguments should include date for interpolation. If !* constant forcing is required, this time can be set to !* zero and only the first step is read. !* Remakrs: AR: At this time the whole file is read but this should be !* replaced by a direct access call to the binary file !* Gulliaume can you do this? It is not that urgent. !********************************************************************** SUBROUTINE READ_SPEC_BINARY_WW3_KERNEL(ISTEP,SPECOUT) USE DATAPOOL IMPLICIT NONE INTEGER, INTENT(IN) :: ISTEP INTEGER :: TMP INTEGER :: IP, IT, TIME(2) REAL(rkind), INTENT(OUT) :: SPECOUT(MSC_WW3,MDC_WW3,NP_WW3) REAL :: SPECOUT_SGLE(MSC_WW3,MDC_WW3) REAL :: FQ_WW3_SNGL(MSC_WW3), DR_WW3_SNGL(MDC_WW3) REAL :: XP_WW3_SGLE(NP_WW3), YP_WW3_SGLE(NP_WW3) REAL :: D(NP_WW3),UA(NP_WW3),UD(NP_WW3),CA(NP_WW3),CD2(NP_WW3) REAL :: m0,m1,m2,df INTEGER :: is,id CHARACTER(LEN=30) :: GNAME CHARACTER(LEN=21) :: LABEL CHARACTER(LEN=10) :: PID(NP_WW3) CALL TEST_FILE_EXIST_DIE('Missing wave file : ', TRIM(WAV%FNAME)) OPEN(WAV%FHNDL, FILE = WAV%FNAME, STATUS = 'OLD',CONVERT='BIG_ENDIAN',FORM='UNFORMATTED') READ(WAV%FHNDL) LABEL,TMP,TMP,TMP,GNAME READ(WAV%FHNDL) FQ_WW3_SNGL READ(WAV%FHNDL) DR_WW3_SNGL IF(LBCSE) THEN ! non-stationary ... DO IT=1,MAXSTEP_WW3 READ(WAV%FHNDL) TIME DO IP = 1, NP_WW3 READ(WAV%FHNDL) PID(IP),YP_WW3_SGLE(IP),XP_WW3_SGLE(IP),D(IP),UA(IP),UD(IP),CA(IP),CD2(IP) ! As if XP and YP would change in time ... well i leave it as it is ... YP_WW3(IP) = YP_WW3_SGLE(IP)*DEGRAD XP_WW3(IP) = XP_WW3_SGLE(IP)*DEGRAD READ(WAV%FHNDL)SPECOUT_SGLE(:,:) ! write(*,*) sum(SPECOUT_SGLE) SPECOUT(:,:,IP) = SPECOUT_SGLE m0 = 0.; m1 = 0.; m2 = 0. DO ID = 1,MDC_WW3-1 DO IS = 1,MSC_WW3-1 DF = FQ_WW3(IS+1)-FQ_WW3(IS) M0 = M0+((SPECOUT_SGLE(IS+1,ID)+SPECOUT_SGLE(IS,ID))/TWO)*DF*DDIR_WW3 M1 = M1+((FQ_WW3(IS+1)-FQ_WW3(IS))/TWO)*((SPECOUT_SGLE(IS+1,ID)+SPECOUT_SGLE(IS,ID))/TWO)*DF*DDIR_WW3 M2 = M2+(((FQ_WW3(IS+1)-FQ_WW3(IS))/TWO)**TWO)*((SPECOUT_SGLE(IS+1,ID)+SPECOUT_SGLE(IS,ID))/TWO)*DF*DDIR_WW3 ! write(*,*) 4*sqrt(m0), m1, m2, df ENDDO ENDDO ENDDO ! IP IF (IT == ISTEP) EXIT ! Read the certain timestep indicated by ISTEP ... ENDDO ! IT ELSE ! stationary ... READ(WAV%FHNDL) TIME DO IP = 1, NP_WW3 READ(WAV%FHNDL)PID(IP),YP_WW3_SGLE(IP),XP_WW3_SGLE(IP),D(IP), UA(IP),UD(IP),CA(IP),CD2(IP) YP_WW3(IP) = YP_WW3_SGLE(IP)*DEGRAD XP_WW3(IP) = XP_WW3_SGLE(IP)*DEGRAD READ(WAV%FHNDL)SPECOUT_SGLE(:,:) SPECOUT(:,:,IP) = SPECOUT_SGLE END DO ENDIF CLOSE(WAV%FHNDL) WRITE(STAT%FHNDL,'("+TRACE...",A)') 'DONE READSPEC2D_WW3' END SUBROUTINE !********************************************************************** !* READ_SPEC_WW3 !* (doing MPI exchanges if MULTIPLE_IN_BOUND = FALSE) !* (otherwise, all node read the same data) !* !* Called by GET_BINARY_WW3_SPECTRA or GET_NC_WW3_SPECTRA !* !* Authors: Mathieu Dutour Sikiric !********************************************************************** SUBROUTINE READ_SPEC_WW3(ISTEP,SPECOUT) USE DATAPOOL IMPLICIT NONE INTEGER, INTENT(IN) :: ISTEP REAL(rkind), INTENT(OUT) :: SPECOUT(MSC_WW3,MDC_WW3,NP_WW3) integer, allocatable :: send_rqst(:) integer, allocatable :: send_stat(:,:) integer siz, iProc # ifndef MPI_PARALL_GRID IF(IBOUNDFORMAT == 3) THEN CALL READ_SPEC_BINARY_WW3_KERNEL(ISTEP,SPECOUT) ELSEIF(IBOUNDFORMAT == 6) THEN CALL READ_SPEC_NC_WW3_KERNEL(ISTEP,SPECOUT) END IF # else IF (MULTIPLE_IN_BOUND) THEN IF(IBOUNDFORMAT == 3) THEN CALL READ_SPEC_BINARY_WW3_KERNEL(ISTEP,SPECOUT) ELSEIF(IBOUNDFORMAT == 6) THEN CALL READ_SPEC_NC_WW3_KERNEL(ISTEP,SPECOUT) END IF ELSE siz=MSC_WW3*MDC_WW3*NP_WW3 IF (myrank .eq. 0) THEN allocate(send_rqst(nproc-1), send_stat(MPI_STATUS_SIZE,nproc-1), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 30') IF(IBOUNDFORMAT == 3) THEN CALL READ_SPEC_BINARY_WW3_KERNEL(ISTEP,SPECOUT) ELSEIF(IBOUNDFORMAT == 6) THEN CALL READ_SPEC_NC_WW3_KERNEL(ISTEP,SPECOUT) END IF DO iProc=2,nproc CALL mpi_isend(SPECOUT, siz, rtype, iProc-1, 2034, comm, send_rqst(iProc-1), ierr) END DO IF (nproc .gt. 1) THEN call mpi_waitall(nproc-1, send_rqst, send_stat,ierr) END IF deallocate(send_rqst, send_stat) ELSE CALL MPI_RECV(SPECOUT, siz, rtype, 0, 2034, comm, istatus, ierr) END IF END IF # endif END SUBROUTINE !********************************************************************** !* GET_BINARY_WW3_SPECTRA !* Read a WAVEWATCHIII binary spectral file and do time and space !* interpolation if required. !* !* Called by WAVE_BOUNDARY_CONDITION !* !* Authors: Aron Roland !* Kai Li !* Guillaume Dodet (18/12/2012) !********************************************************************** SUBROUTINE GET_BINARY_WW3_SPECTRA(WBACOUT) USE DATAPOOL, ONLY: NP_WW3, rkind, DR_WW3, DDIR_WW3, FQ_WW3, FRLOW, LNANINFCHK, DBG, FRHIGH USE DATAPOOL, ONLY: LINHOM, IWBNDLC, XP, YP, XP_WW3, YP_WW3, STAT, MSC, MDC, IWBMNP, MSC_WW3 USE DATAPOOL, ONLY: MDC_WW3 # ifdef MPI_PARALL_GRID USE DATAPOOL, ONLY: XLON, YLAT # endif IMPLICIT NONE REAL(rkind), INTENT(OUT) :: WBACOUT(MSC,MDC,IWBMNP) INTEGER :: IB,IPGL,IBWW3,TIME(2),IS REAL(rkind) :: SPEC_WW3(MSC_WW3,MDC_WW3,NP_WW3),SPEC_WWM(MSC,MDC,NP_WW3) REAL(rkind) :: DIST(NP_WW3),TMP(NP_WW3), INDBWW3(NP_WW3) REAL(rkind) :: SPEC_WW3_TMP(MSC_WW3,MDC_WW3,NP_WW3),SPEC_WW3_UNSORT(MSC_WW3,MDC_WW3,NP_WW3) REAL(rkind) :: JUNK(MDC_WW3),DR_WW3_UNSORT(MDC_WW3),DR_WW3_TMP(MDC_WW3) REAL(rkind) :: XP_WWM,YP_WWM INTEGER :: IFILE, IT WRITE(STAT%FHNDL,'("+TRACE...",A)') 'Begin GETWW3SPECTRA' CALL COMPUTE_IFILE_IT(IFILE, IT) ! ! Read spectra in file ! CALL READ_SPEC_WW3(IT,SPEC_WW3_UNSORT) ! ! Sort directions and carries spectra along (ww3 directions are not ! montonic) ! DO IBWW3 = 1, NP_WW3 DO IS = 1,MSC_WW3 DR_WW3_TMP=DR_WW3 SPEC_WW3_TMP(IS,:,IBWW3) = SPEC_WW3_UNSORT(IS,:,IBWW3) CALL SSORT2 (DR_WW3_TMP,SPEC_WW3_TMP(IS,:,IBWW3),JUNK,MDC_WW3,2) SPEC_WW3(IS,:,IBWW3) = SPEC_WW3_TMP(IS,:,IBWW3) DR_WW3 = DR_WW3_TMP ENDDO DDIR_WW3 = DR_WW3(2) - DR_WW3(1) ! WRITE(STAT%FHNDL,*) 'AFTER SORTING', IBWW3, SUM(SPEC_WW3(:,:,IBWW3)) ENDDO ! ! Interpolate ww3 spectra on wwm frequency grid ! GD: at the moment 360º spanning grids are mandatory ! IF((FQ_WW3(1).GT.FRLOW).OR.(FQ_WW3(MSC_WW3).LT.FRHIGH)) THEN ! WRITE(STAT%FHNDL,*)'WW3 FMIN = ',FQ_WW3(1),'WWM FMIN = ',FRLOW ! WRITE(STAT%FHNDL,*)'WW3 FMAX = ',FQ_WW3(MSC_WW3),'WWM FMAX = ', FRHIGH ! WRITE(STAT%FHNDL,*)'WW3 spectra does not encompass the whole WWM spectra, please carefully check if this makes sense for your simulations' CALL SPECTRALINT(SPEC_WW3,SPEC_WWM) ELSE ! WRITE(STAT%FHNDL,*)'WW3 FMIN = ',FQ_WW3(1),'WWM FMIN = ',FRLOW ! WRITE(STAT%FHNDL,*)'WW3 FMAX = ',FQ_WW3(MSC_WW3),'WWM FMAX = ', FRHIGH CALL SPECTRALINT(SPEC_WW3,SPEC_WWM) ENDIF IF (LNANINFCHK) THEN WRITE(DBG%FHNDL,*) 'SUMS AFTER INTERPOLATION', SUM(SPEC_WW3),SUM(SPEC_WWM) ENDIF ! ! Interpolate ww3 spectra on wwm boundary nodes ! GD: ww3 forcing works until here. Some more debugging is needed ! TMP = 0 IF(LINHOM) THEN !nearest-neighbour interpolation DO IB=1,IWBMNP IPGL = IWBNDLC(IB) # ifdef SCHISM XP_WWM=XLON(IPGL)! * RADDEG YP_WWM=YLAT(IPGL)! * RADDEG # else XP_WWM = XP(IPGL) YP_WWM = YP(IPGL) # endif IF (NP_WW3 .GT. 1) THEN DO IBWW3=1,NP_WW3 DIST(IBWW3)=SQRT((XP_WWM-XP_WW3(IBWW3))**2+(YP_WWM-YP_WW3(IBWW3))**2) INDBWW3(IBWW3)=IBWW3 ENDDO CALL SSORT2 (DIST, INDBWW3, TMP, NP_WW3, 2) CALL SHEPARDINT2D(2, 1./DIST(1:2),MSC,MDC,SPEC_WWM(:,:,INT(INDBWW3(1:2))), WBACOUT(:,:,IB), 1) !WRITE(STAT%FHNDL,'(A20, 2F20.5,3F30.10)') ' AFTER INTERPOLATION ', INDBWW3(1), INDBWW3(2), sum(SPEC_WWM(:,:,INT(INDBWW3(1)))), sum(SPEC_WWM(:,:,INT(INDBWW3(2)))), SUM(WBACOUT(:,:,IB)) ELSE WBACOUT(:,:,IB) = SPEC_WWM(:,:,1) ENDIF ENDDO ! IB ELSE DO IB=1,IWBMNP WBACOUT(:,:,IB) = SPEC_WWM(:,:,1) ENDDO ENDIF WRITE(STAT%FHNDL,'("+TRACE...",A)') 'DONE GETWW3SPECTRA' END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE INIT_BINARY_WW3_SPECTRA USE DATAPOOL IMPLICIT NONE ! ! Read header to get grid dimension, frequencies and directions, ! output locations and first time step ! WRITE(STAT%FHNDL,'("+TRACE...",A)') 'START WITH INIT_BINARY_WW3_SPECTRA' CALL READSPEC2D_BINARY_WW3_INIT_SPEC CALL READSPEC2D_BINARY_WW3_INIT_TIME WRITE(STAT%FHNDL,'("+TRACE...",A)') 'DONE WITH INIT_BINARY_WW3_SPECTRA' END SUBROUTINE !********************************************************************** !* INIT_NC_WW3_SPECTRA !* Read the header of a WAVEWATCHIII netcdf spectral file to get !* dimensions of the spectral grid and number of output locations !* required for dynamic allocation. !* !* Called by INIT_WAVE_BOUNDARY_CONDITION !* !* Authors: Kevin Martins (07/2018) !********************************************************************** #ifdef NCDF SUBROUTINE INIT_NC_WW3_SPECTRA USE DATAPOOL USE NETCDF IMPLICIT NONE INTEGER :: NCID INTEGER :: FREQ_DIM_ID, DIRS_DIM_ID, TIME_DIM_ID, STAT_DIM_ID INTEGER :: FREQ_VAR_ID, DIRS_VAR_ID, TIME_VAR_ID, STAT_LAT_ID, STAT_LON_ID REAL(rkind) :: eJD_WW3, eJD_WWM !------------------------------------ ! Reading grid and station dimensions !------------------------------------ ! Opening netcdf file CALL TEST_FILE_EXIST_DIE("Missing WW3 boundary file : ", TRIM(WAV%FNAME)) ISTAT = NF90_OPEN(WAV%FNAME,NF90_NOWRITE,NCID) CALL GENERIC_NETCDF_ERROR_WWM("ERROR WHILE OPENING NETCDF FILE",1,ISTAT) ! Reading frequency dimension ISTAT = NF90_INQ_DIMID(NCID,'frequency',FREQ_DIM_ID) CALL GENERIC_NETCDF_ERROR_WWM("ERROR WHILE READING FREQUENCY DIMENSION ID",2,ISTAT) ISTAT = NF90_INQUIRE_DIMENSION(NCID,FREQ_DIM_ID,len = MSC_WW3) CALL GENERIC_NETCDF_ERROR_WWM("ERROR WHILE READING FREQUENCY DIMENSION",3,ISTAT) ! Reading direction dimension ISTAT = NF90_INQ_DIMID(NCID,'direction',DIRS_DIM_ID) CALL GENERIC_NETCDF_ERROR_WWM("ERROR WHILE READING DIRECTION DIMENSION ID",4,ISTAT) ISTAT = NF90_INQUIRE_DIMENSION(NCID,DIRS_DIM_ID,len = MDC_WW3) CALL GENERIC_NETCDF_ERROR_WWM("ERROR WHILE READING DIRECTION DIMENSION",5,ISTAT) ! Reading time dimension ISTAT = NF90_INQ_DIMID(NCID,'time',TIME_DIM_ID) CALL GENERIC_NETCDF_ERROR_WWM("ERROR WHILE READING DIRECTION DIMENSION ID",6,ISTAT) ISTAT = NF90_INQUIRE_DIMENSION(NCID,TIME_DIM_ID,len = MAXSTEP_WW3) CALL GENERIC_NETCDF_ERROR_WWM("ERROR WHILE READING DIRECTION DIMENSION",7,ISTAT) ! Reading station dimension ISTAT = NF90_INQ_DIMID(NCID,'station',STAT_DIM_ID) CALL GENERIC_NETCDF_ERROR_WWM("ERROR WHILE READING STATION DIMENSION ID",8,ISTAT) ISTAT = NF90_INQUIRE_DIMENSION(NCID,STAT_DIM_ID,len = NP_WW3) CALL GENERIC_NETCDF_ERROR_WWM("ERROR WHILE READING STATION DIMENSION",9,ISTAT) !------------------------------------ ! Reading variables !------------------------------------ ! Arrays for the WW3 spectra ALLOCATE(FQ_WW3(MSC_WW3),DR_WW3(MDC_WW3),BND_TIME_ALL_FILES(1,MAXSTEP_WW3), XP_WW3(NP_WW3), YP_WW3(NP_WW3),stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_bdcons, WW3 freq/dirs allocate error') ! Reading longitude ISTAT = NF90_INQ_VARID(NCID,'longitude',STAT_LON_ID) CALL GENERIC_NETCDF_ERROR_WWM("ERROR WHILE READING LONGITUDE VARIABLE",10,ISTAT) ISTAT = NF90_GET_VAR(NCID,STAT_LON_ID,XP_WW3,start=(/1,1/),count = (/NP_WW3,1/)) CALL GENERIC_NETCDF_ERROR_WWM("ERROR WHILE READING LONGITUDE",11, ISTAT) ! Reading latitude ISTAT = NF90_INQ_VARID(NCID,'latitude',STAT_LAT_ID) CALL GENERIC_NETCDF_ERROR_WWM("ERROR WHILE READING LATITUDE VARIABLE",12,ISTAT) ISTAT = NF90_GET_VAR(NCID,STAT_LAT_ID,YP_WW3,start=(/1,1/),count = (/NP_WW3,1/)) CALL GENERIC_NETCDF_ERROR_WWM("ERROR WHILE READING LATITUDE",13, ISTAT) ! Convert coordinates to RAD XP_WW3 = XP_WW3*DEGRAD YP_WW3 = YP_WW3*DEGRAD ! Reading frequencies ISTAT = NF90_INQ_VARID(NCID,'frequency',FREQ_VAR_ID) CALL GENERIC_NETCDF_ERROR_WWM("ERROR WHILE READING FREQUENCY VARIABLE",14,ISTAT) ISTAT = NF90_GET_VAR(NCID,FREQ_VAR_ID,FQ_WW3,start=(/1/),count = (/MSC_WW3/)) CALL GENERIC_NETCDF_ERROR_WWM("ERROR WHILE READING FREQUENCIES",15, ISTAT) ! Reading directions ISTAT = NF90_INQ_VARID(NCID,'direction',DIRS_VAR_ID) CALL GENERIC_NETCDF_ERROR_WWM("ERROR WHILE READING DIRECTION VARIABLE",16,ISTAT) ISTAT = NF90_GET_VAR(NCID,DIRS_VAR_ID,DR_WW3,start=(/1/),count = (/MDC_WW3/)) CALL GENERIC_NETCDF_ERROR_WWM("ERROR WHILE READING DIRECTIONS",17, ISTAT) ! Reading time ISTAT = NF90_INQ_VARID(NCID,'time',TIME_VAR_ID) CALL GENERIC_NETCDF_ERROR_WWM("ERROR WHILE READING DIRECTION VARIABLE",18,ISTAT) ISTAT = NF90_GET_VAR(NCID,TIME_VAR_ID,BND_TIME_ALL_FILES,start=(/1/),count = (/MAXSTEP_WW3/)) CALL GENERIC_NETCDF_ERROR_WWM("ERROR WHILE READING DIRECTIONS",19,ISTAT) !------------------------------------ ! Conventions (direction and time) !------------------------------------ ! Convention for directions: care must be taken here as sometimes the directions in WW3 ! are expressed with the origin taken towards the East. This corresponds to the last version of WW3. IF(MAXVAL(DR_WW3)>300) THEN DR_WW3 = MOD(DR_WW3*PI/180,PI2) ! Simple conversion to RAD ENDIF ! Converting time from WW3 convention to WWM convention ! In WW3, time is given in Julian days, with reference from [1990,1,1,0,0,0] CALL DATE2JD(1858,11,17,0,0,0,eJD_WWM) CALL DATE2JD(1990,1,1,0,0,0,eJD_WW3) BND_TIME_ALL_FILES = BND_TIME_ALL_FILES + eJD_WW3 - eJD_WWM DTBOUND_WW3 = (BND_TIME_ALL_FILES(1,2)-BND_TIME_ALL_FILES(1,1))*DAY2SEC SEBO%DELT = DTBOUND_WW3 SEBO%BMJD = BND_TIME_ALL_FILES(1,1) SEBO%EMJD = BND_TIME_ALL_FILES(1,MAXSTEP_WW3) ! Closing the netcdf file ISTAT = NF90_CLOSE(ncid) CALL GENERIC_NETCDF_ERROR_WWM("ERROR WHILE CLOSING NETCDF FILE",8,ISTAT) END SUBROUTINE #endif !********************************************************************** !* READ_SPEC_NC_WW3_KERNEL, Called by READ_SPEC_WW3 !* Reads spectra in WAVEWATCHIII netcdf spectral file !* Authors: Kevin MARTINS !********************************************************************** #ifdef NCDF SUBROUTINE READ_SPEC_NC_WW3_KERNEL(ISTEP,SPECOUT) USE DATAPOOL USE NETCDF IMPLICIT NONE REAL(rkind), INTENT(OUT) :: SPECOUT(MSC_WW3,MDC_WW3,NP_WW3) INTEGER, INTENT(IN) :: ISTEP REAL :: SPEC_WW3(MDC_WW3,MSC_WW3) INTEGER :: IP, NCID, ENERGY_VAR_ID, DD, FF ! Opening netcdf file ISTAT = NF90_OPEN(WAV%FNAME,NF90_NOWRITE,NCID) CALL GENERIC_NETCDF_ERROR_WWM("ERROR WHILE OPENING NETCDF FILE",1,ISTAT) ! Loading the energy density variable ID ISTAT = NF90_INQ_VARID(NCID,'efth',ENERGY_VAR_ID) CALL GENERIC_NETCDF_ERROR_WWM("ERROR WHILE READING SPECTRA DIMENSION",2,ISTAT) ! Reading spectra IF(LBCSE) THEN ! non-stationary ... DO IP = 1,NP_WW3 ISTAT = NF90_GET_VAR(NCID,ENERGY_VAR_ID,SPEC_WW3,start=(/1,1,IP,ISTEP/),count = (/MDC_WW3,MSC_WW3,1,1/)) CALL GENERIC_NETCDF_ERROR_WWM("ERROR WHILE READING SPECTRA",3,ISTAT) ! Re-ordering spectrum (WW3 is [MDC_WW3,MSC_WW3] while WWM is [MSC_WW3,MDC_WW3]) SPECOUT(:,:,IP) = TRANSPOSE(SPEC_WW3) ENDDO ! IP ELSE ! stationary ... DO IP = 1,NP_WW3 ISTAT = NF90_GET_VAR(NCID,ENERGY_VAR_ID,SPEC_WW3,start=(/1,1,IP,1/),count = (/MDC_WW3,MSC_WW3,1,1/)) CALL GENERIC_NETCDF_ERROR_WWM("ERROR WHILE READING SPECTRA",4,ISTAT) ! Re-ordering spectrum (WW3 is [MDC_WW3,MSC_WW3] while WWM is [MSC_WW3,MDC_WW3]) SPECOUT(:,:,IP) = TRANSPOSE(SPEC_WW3) END DO ENDIF ISTAT = NF90_CLOSE(NCID) CALL GENERIC_NETCDF_ERROR_WWM("ERROR WHILE CLOSING NETCDF FILE",5,ISTAT) END SUBROUTINE #endif !********************************************************************** !* GET_NC_WW3_SPECTRA !* Read a WAVEWATCHIII netcdf spectral file and do time and space !* interpolation if required. !* !* Called by WAVE_BOUNDARY_CONDITION !* !* Authors: Kevin Martins (18/12/2012) !********************************************************************** #ifdef NCDF SUBROUTINE GET_NC_WW3_SPECTRA(WBACOUT) USE DATAPOOL, ONLY: NP_WW3, rkind, DR_WW3, DDIR_WW3, FQ_WW3, FRLOW, LNANINFCHK, DBG, FRHIGH USE DATAPOOL, ONLY: LINHOM, IWBNDLC, XP, YP, XP_WW3, YP_WW3, STAT, MSC, MDC, IWBMNP, MSC_WW3 USE DATAPOOL, ONLY: MDC_WW3, SPDIR, PI, PI2 USE NETCDF # ifdef MPI_PARALL_GRID USE DATAPOOL, ONLY: XLON, YLAT # endif IMPLICIT NONE REAL(rkind), INTENT(OUT) :: WBACOUT(MSC,MDC,IWBMNP) INTEGER :: IB,IPGL,IBWW3,TIME(2),IS REAL(rkind) :: SPEC_WW3(MSC_WW3,MDC_WW3,NP_WW3),SPEC_WWM(MSC,MDC,NP_WW3) REAL(rkind) :: DIST(NP_WW3),TMP(NP_WW3), INDBWW3(NP_WW3) REAL(rkind) :: SPEC_WW3_TMP(MSC_WW3,MDC_WW3,NP_WW3),SPEC_WW3_UNSORT(MSC_WW3,MDC_WW3,NP_WW3) REAL(rkind) :: JUNK(MDC_WW3),DR_WW3_UNSORT(MDC_WW3),DR_WW3_TMP(MDC_WW3) REAL(rkind) :: XP_WWM,YP_WWM INTEGER :: IFILE, IT ! Computing time step CALL COMPUTE_IFILE_IT(IFILE, IT) ! Read spectra in file CALL READ_SPEC_WW3(IT,SPEC_WW3_UNSORT) ! Sort directions and carries spectra along (ww3 directions are not monotonic) DO IBWW3 = 1, NP_WW3 DO IS = 1,MSC_WW3 DR_WW3_TMP=DR_WW3 SPEC_WW3_TMP(IS,:,IBWW3) = SPEC_WW3_UNSORT(IS,:,IBWW3) CALL SSORT2 (DR_WW3_TMP,SPEC_WW3_TMP(IS,:,IBWW3),JUNK,MDC_WW3,2) SPEC_WW3(IS,:,IBWW3) = SPEC_WW3_TMP(IS,:,IBWW3) DR_WW3 = DR_WW3_TMP ENDDO ENDDO ! Interpolate ww3 spectra on WWM frequency grid ! GD: at the moment 360º spanning grids are mandatory IF((FQ_WW3(1).GT.FRLOW).OR.(FQ_WW3(MSC_WW3).LT.FRHIGH)) THEN ! WRITE(STAT%FHNDL,*)'WW3 FMIN = ',FQ_WW3(1),'WWM FMIN = ',FRLOW ! WRITE(STAT%FHNDL,*)'WW3 FMAX = ',FQ_WW3(MSC_WW3),'WWM FMAX = ', FRHIGH ! WRITE(STAT%FHNDL,*)'WW3 spectra does not encompass the whole WWM spectra, please carefully check if this makes sense for your simulations' CALL SPECTRALINT(SPEC_WW3,SPEC_WWM) ELSE ! WRITE(STAT%FHNDL,*)'WW3 FMIN = ',FQ_WW3(1),'WWM FMIN = ',FRLOW ! WRITE(STAT%FHNDL,*)'WW3 FMAX = ',FQ_WW3(MSC_WW3),'WWM FMAX = ', FRHIGH CALL SPECTRALINT(SPEC_WW3,SPEC_WWM) ENDIF ! Interpolate ww3 spectra on wwm boundary nodes ! GD: ww3 forcing works until here. Some more debugging is needed (e.g. for LINHOM = F) TMP = 0 IF(LINHOM) THEN !nearest-neighbour interpolation DO IB=1,IWBMNP IPGL = IWBNDLC(IB) # ifdef SCHISM XP_WWM = XLON(IPGL)! * RADDEG YP_WWM = YLAT(IPGL)! * RADDEG # else XP_WWM = XP(IPGL) YP_WWM = YP(IPGL) # endif IF (NP_WW3 .GT. 1) THEN DO IBWW3=1,NP_WW3 DIST(IBWW3)=SQRT((XP_WWM-XP_WW3(IBWW3))**2+(YP_WWM-YP_WW3(IBWW3))**2) INDBWW3(IBWW3)=IBWW3 ENDDO CALL SSORT2 (DIST, INDBWW3, TMP, NP_WW3, 2) CALL SHEPARDINT2D(2, 1./DIST(1:2),MSC,MDC,SPEC_WWM(:,:,INT(INDBWW3(1:2))), WBACOUT(:,:,IB), 1) !WRITE(STAT%FHNDL,'(A20, 2F20.5,3F30.10)') ' AFTER INTERPOLATION ', INDBWW3(1), INDBWW3(2), sum(SPEC_WWM(:,:,INT(INDBWW3(1)))), sum(SPEC_WWM(:,:,INT(INDBWW3(2)))), SUM(WBACOUT(:,:,IB)) ELSE WBACOUT(:,:,IB) = SPEC_WWM(:,:,1) ENDIF ENDDO ! IB ELSE DO IB=1,IWBMNP WBACOUT(:,:,IB) = SPEC_WWM(:,:,1) ENDDO ENDIF WRITE(STAT%FHNDL,'("+TRACE...",A)') 'DONE GETWW3SPECTRA' END SUBROUTINE #endif !********************************************************************** !* * !********************************************************************** SUBROUTINE SHEPARDINT2D(NP,WEIGHT,D1,D2,Z,ZINT,P) USE DATAPOOL, ONLY: rkind IMPLICIT NONE INTEGER, INTENT(IN) :: NP,P,D1,D2 REAL(rkind), INTENT(IN) :: WEIGHT(NP), Z(D1,D2,NP) REAL(rkind), INTENT(OUT) :: ZINT(D1,D2) INTEGER :: IP REAL :: SW SW=0 ZINT=0 DO IP=1,NP SW=SW+WEIGHT(IP)**P ZINT(:,:)=ZINT(:,:)+Z(:,:,IP)*(WEIGHT(IP)**P) ENDDO ZINT=ZINT/SW END SUBROUTINE !********************************************************************** !* SPECTRALINT !* Interpolate spectrum on WWM spectral grid. !* !* Called by GETWW3SPECTRA !* !* Authors: Aron Roland !* Kai Li !* Guillaume Dodet (18/12/2012) !* !* Remarks: GD)The spectral grid needs to be define over 360º and !* the min (resp. max) frequencies need to be smaller !* (resp. higher) than in WWM frequency grid. !********************************************************************** SUBROUTINE SPECTRALINT(SPEC_WW3,SPEC_WWM) USE DATAPOOL IMPLICIT NONE REAL(rkind), INTENT(IN) :: SPEC_WW3(MSC_WW3,MDC_WW3,NP_WW3) REAL(rkind), INTENT(OUT) :: SPEC_WWM(MSC,MDC,NP_WW3) REAL(rkind) :: SPEC_WW3_TMP(MSC_WW3,MDC,NP_WW3) REAL(rkind) :: DF, M0_WW3, M1_WW3, M2_WW3, M0_WWM, M1_WWM, M2_WWM INTEGER :: IP,IS,ID REAL(rkind) :: JACOBIAN(MSC), AM, SM WRITE(STAT%FHNDL,'("+TRACE...",A)') 'ENTERING SPECTRALINT' JACOBIAN = ONE/(SPSIG*PI2)! ENERGY / HZ -> ACTION / RAD SPEC_WW3_TMP = ZERO SPEC_WWM = ZERO IF (LNANINFCHK) THEN WRITE(DBG%FHNDL,'(A20,I10,3F30.2)') 'BEFORE INTERPOLATION', SUM(SPEC_WW3), SUM(SPEC_WW3_TMP), SUM(SPEC_WWM) ENDIF DO IP=1,NP_WW3 DO IS=1,MSC_WW3 CALL INTERLIND(MDC_WW3,MDC,DR_WW3,SPDIR,SPEC_WW3(IS,:,IP),SPEC_WW3_TMP(IS,:,IP)) ENDDO DO ID=1,MDC CALL INTERLIN (MSC_WW3,MSC,FQ_WW3,FR,SPEC_WW3_TMP(:,ID,IP),SPEC_WWM(:,ID,IP)) ENDDO M0_WW3 = ZERO; M1_WW3 = ZERO; M2_WW3 = ZERO DO ID = 1,MDC_WW3 DO IS = 1,MSC_WW3-1 DF = FQ_WW3(IS+1)-FQ_WW3(IS) AM = (SPEC_WW3(IS+1,ID,IP)+SPEC_WW3(IS,ID,IP))/TWO SM = (FQ_WW3(IS+1)+FQ_WW3(IS))/TWO M0_WW3 =M0_WW3+AM*DF*DDIR_WW3 M1_WW3 =M1_WW3+AM*SM*DF*DDIR_WW3 M2_WW3 =M2_WW3+AM*SM**2*DF*DDIR_WW3 ENDDO ENDDO M0_WWM = ZERO; M1_WWM = ZERO; M2_WWM = ZERO DO ID = 1,MDC DO IS = 1,MSC-1 DF = FR(IS+1)-FR(IS) AM = (SPEC_WWM(IS+1,ID,IP)+SPEC_WWM(IS,ID,IP))/TWO SM = (FR(IS+1)+FR(IS))/TWO M0_WWM =M0_WWM+AM*DF*DDIR M1_WWM =M1_WWM+AM*SM*DF*DDIR M2_WWM =M2_WWM+AM*SM**2*DF*DDIR ENDDO ENDDO WRITE(STAT%FHNDL,*) 'POINT NUMBER', IP WRITE(STAT%FHNDL,'(A10,2F20.10,A10,2F20.10)') 'M1 = ',M1_WW3, M1_WWM, 'M2 = ',M2_WW3, M2_WWM END DO IF (LNANINFCHK) THEN WRITE(DBG%FHNDL,'(A20,I10,3F30.2)') 'AFTER INTERPOLATION', IP, SUM(SPEC_WW3), SUM(SPEC_WW3_TMP), SUM(SPEC_WWM) ENDIF ! Do jacobian DO IP = 1, NP_WW3 DO ID = 1, MDC SPEC_WWM(:,ID,IP) = SPEC_WWM(:,ID,IP) * JACOBIAN(:) ! convert to wave action in sigma,theta space END DO END DO WRITE(STAT%FHNDL,*)'CHECKING INTEGRATED PARAMETERS AFTER JACOBIAN' DO IP = 1, NP_WW3 M0_WW3 = ZERO; M1_WW3 = ZERO; M2_WW3 = ZERO DO ID = 1,MDC_WW3 DO IS = 1,MSC_WW3-1 DF = FQ_WW3(IS+1)-FQ_WW3(IS) AM = (SPEC_WW3(IS+1,ID,IP)+SPEC_WW3(IS,ID,IP))/TWO SM = (FQ_WW3(IS+1)+FQ_WW3(IS))/TWO M0_WW3 =M0_WW3+AM*DF*DDIR_WW3 M1_WW3 =M1_WW3+AM*SM*DF*DDIR_WW3 M2_WW3 =M2_WW3+AM*SM**2*DF*DDIR_WW3 ENDDO ENDDO M0_WWM = ZERO; M1_WWM = ZERO; M2_WWM = ZERO DO ID = 1,MDC DO IS = 1,MSC-1 DF = SPSIG(IS+1)-SPSIG(IS) SM = (SPSIG(IS+1)+SPSIG(IS))/TWO AM = (SPEC_WWM(IS+1,ID,IP)+SPEC_WWM(IS,ID,IP))/TWO * SM M0_WWM =M0_WWM+AM*DF*DDIR M1_WWM =M1_WWM+AM*SM*DF*DDIR M2_WWM =M2_WWM+AM*SM**2*DF*DDIR ENDDO ENDDO WRITE(STAT%FHNDL,*) 'POINT NUMBER', IP WRITE(STAT%FHNDL,'(A10,2F20.10,A10,2F20.10)') 'M1 = ',M1_WW3, M1_WWM, 'M2 = ',M2_WW3, M2_WWM END DO IF (LNANINFCHK) THEN WRITE(DBG%FHNDL,'(A20,I10,3F30.2)') 'AFTER JACOBIAN', IP, SUM(SPEC_WW3), SUM(SPEC_WW3_TMP), SUM(SPEC_WWM) ENDIF END SUBROUTINE !********************************************************************** !* INTERLIND !* Interpolate vector on a 2-pi periodic axis (directions in wwm grid). !* !* Called by SPECTRALINT !* !* Authors: Aron Roland !* Kai Li !* Guillaume Dodet (18/12/2012) !* !* Remarks: GD)This routine should be togeher with INTERLIN either here !* or in wwm_aux.F90. !********************************************************************** SUBROUTINE INTERLIND (NX1, NX2, X1, X2, Y1, Y2) USE DATAPOOL, ONLY : THR, rkind,PI2 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) DO I = 1, NX1 - 1 DX1(I) = X1(I+1)-X1(I) END DO DX1(NX1) = X1(1)+PI2-X1(NX1) DO I = 1, NX2 DO J = 1,NX1 - 1 IF (ABS(X2(I)-X1(J)) .LT. THR) THEN Y2(I) = Y1(J) EXIT 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)) EXIT ELSE IF (X2(I) .GT. X1(NX1)) THEN Y2(I) = Y1(NX1)+(Y1(1)-Y1(NX1))/DX1(NX1)*(X2(I)-X1(NX1)) EXIT ELSE IF (X2(I) .LT. X1(1)) THEN Y2(I) = Y1(NX1)+(Y1(1)-Y1(NX1))/DX1(NX1)*(X2(I)+PI2-X1(NX1)) EXIT END IF END DO END DO END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE ALLOC_SPEC_BND() USE DATAPOOL IMPLICIT NONE IF (LINHOM) THEN ALLOCATE (SFRQ(WBMSC,IWBMNP), SPEG(WBMSC,WBMDC,IWBMNP), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_bdcons, allocate error 22') ALLOCATE (SDIR(WBMSC,IWBMNP), SPRD(WBMSC,IWBMNP), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_bdcons, allocate error 23') ELSE ALLOCATE (SFRQ(WBMSC,1), SPEG(WBMSC,WBMDC,1), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_bdcons, allocate error 24') ALLOCATE (SDIR(WBMSC,1), SPRD(WBMSC,1), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_bdcons, allocate error 25') END IF END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE READWAVEPARWW3 USE DATAPOOL IMPLICIT NONE INTEGER :: IP # ifdef MPI_PARALL_GRID INTEGER :: IPP REAL(rkind) :: RTMP # endif ! SPPARM(1): Hs, sign. wave height ! SPPARM(2): Wave period given by user (either peak or mean) ! SPPARM(3): average direction ! SPPARM(4): directional spread ! SPPARM(5): spectral shape (1-4), ! (1 - Pierson-Moskowitz, ! 2 - JONSWAP, ! 3 - BIN, ! 4 - Gauss) ! negative peak (+) ! or mean frequency (-) ! SPPARM(6): directional spreading in degree (1) or exponent (2) ! SPPARM(7): gaussian width for the gauss spectrum 0.1 ! SPPARM(8): peak enhancement factor for the JONSWAP spectra 3.3 SPPARM(4,:) = 20. SPPARM(5,:) = 2. SPPARM(6,:) = 2. SPPARM(7,:) = 0.1 SPPARM(8,:) = 3.3 IF (LINHOM) THEN READ(WAV%FHNDL,*) END IF # ifdef MPI_PARALL_GRID IPP = 0 IF (LINHOM) THEN DO IP = 1, IWBMNPGL IF(ipgl(IWBNDGL(IP))%rank == myrank) THEN ! IF boundary nodes belong to local domain read values into boundary array IPP = IPP + 1 READ (WAV%FHNDL, *) SPPARM(1,IPP), SPPARM(2,IPP), SPPARM(3,IPP) ELSE READ (WAV%FHNDL, *) RTMP, RTMP, RTMP ! ELSE ... throw them away ... ENDIF END DO ELSE READ (WAV%FHNDL, *) SPPARM(1,1), SPPARM(2,1), SPPARM(3,1) DO IP = 1, IWBMNPGL SPPARM(1:3,IP) = SPPARM(1:3,1) END DO END IF # else IF (LINHOM) THEN DO IP = 1, IWBMNP READ (WAV%FHNDL, *) SPPARM(1,IP), SPPARM(2,IP), SPPARM(3,IP) END DO ELSE READ (WAV%FHNDL, *) SPPARM(1,1), SPPARM(2,1), SPPARM(3,1) DO IP = 1, IWBMNP SPPARM(1:3,IP) = SPPARM(1:3,1) END DO END IF # endif END SUBROUTINE !********************************************************************** !* * !********************************************************************** #ifdef NCDF SUBROUTINE WRITE_NETCDF_BOUND_HEADERS_1(FILE_NAME, nbTime, np_write, nbBound, BOUC_PARAM, BOUC_SPEC) USE DATAPOOL USE NETCDF implicit none character(len=140), intent(in) :: FILE_NAME integer, intent(in) :: nbTime, np_write, nbBound logical, intent(in) :: BOUC_PARAM, BOUC_SPEC ! character (len = *), parameter :: UNITS = "units" integer one_dims, two_dims, three_dims, fifteen_dims integer mnp_dims, mne_dims, nfreq_dims, ndir_dims, eight_dims integer iret, var_id, ncid integer ntime_dims, iwbmnpgl_dims character (len = *), parameter :: CallFct="WRITE_NETCDF_BOUND_HEADERS_1" iret = nf90_create(TRIM(FILE_NAME), NF90_CLOBBER, ncid) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 1, iret) iret = nf90_def_dim(ncid, 'one', 1, one_dims) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 2, iret) iret = nf90_def_dim(ncid, 'two', 2, two_dims) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 3, iret) iret = nf90_def_dim(ncid, 'three', 3, three_dims) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 4, iret) iret = nf90_def_dim(ncid, 'fifteen', 15, fifteen_dims) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 5, iret) iret = nf90_def_dim(ncid, 'np_total', np_write, mnp_dims) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 6, iret) iret = nf90_def_dim(ncid, 'nfreq', MSC, nfreq_dims) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 7, iret) iret = nf90_def_dim(ncid, 'ndir', MDC, ndir_dims) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 8, iret) iret = nf90_def_dim(ncid, 'eight', 8, eight_dims) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 9, iret) iret = nf90_def_dim(ncid, 'IWBMNPGL', nbBound, iwbmnpgl_dims) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 10, iret) ! CALL WRITE_PARAM_1(ncid, one_dims) ! CALL WRITE_NETCDF_TIME_HEADER(ncid, nbTime, ntime_dims) ! iret=nf90_def_var(ncid,'IOBP',NF90_INT,(/ mnp_dims /), var_id) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 11, iret) iret=nf90_put_att(ncid,var_id,UNITS,'integer') CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 12, iret) iret=nf90_put_att(ncid,var_id,'description','boundary status of nodes') CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 13, iret) iret=nf90_put_att(ncid,var_id,'case 0','interior point') CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 14, iret) iret=nf90_put_att(ncid,var_id,'case 1','island') CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 15, iret) iret=nf90_put_att(ncid,var_id,'case 2','dirichlet condition') CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 16, iret) iret=nf90_put_att(ncid,var_id,'case 3','neumann condition') CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 17, iret) iret=nf90_put_att(ncid,var_id,'case 4','unknown') CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 18, iret) ! iret=nf90_def_var(ncid,'IWBNDGL',NF90_INT,(/ iwbmnpgl_dims /), var_id) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 19, iret) iret=nf90_put_att(ncid,var_id,'description','indices of boundary nodes') CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 20, iret) ! IF (BOUC_NETCDF_OUT_PARAM) THEN iret=nf90_def_var(ncid,'SPPARM',NF90_OUTTYPE_BOUC,(/ eight_dims, iwbmnpgl_dims, ntime_dims /), var_id) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 21, iret) iret=nf90_put_att(ncid,var_id,'description','Parametric boundary condition') CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 22, iret) END IF ! IF (BOUC_NETCDF_OUT_SPECTRA) THEN iret=nf90_def_var(ncid,'WBAC',NF90_OUTTYPE_BOUC,(/ nfreq_dims, ndir_dims, iwbmnpgl_dims, ntime_dims /), var_id) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 23, iret) iret=nf90_put_att(ncid,var_id,'description','boundary wave action') CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 24, iret) END IF iret=nf90_close(ncid) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 25, iret) END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE WRITE_NETCDF_BOUND_HEADERS_2(FILE_NAME, np_write, IOBPwrite, nbBound, ListBound) USE DATAPOOL USE NETCDF implicit none character(len=140), intent(in) :: FILE_NAME integer, intent(in) :: np_write integer, intent(in) :: IOBPwrite(np_write) integer, intent(in) :: nbBound integer, intent(in) :: ListBound(nbBound) ! integer ncid integer iret, var_id character (len = *), parameter :: CallFct="WRITE_NETCDF_BOUND_HEADERS_2" ! iret=nf90_open(TRIM(FILE_NAME), NF90_WRITE, ncid) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 1, iret) ! CALL WRITE_PARAM_2(ncid) ! iret=nf90_inq_varid(ncid, "IOBP", var_id) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 2, iret) iret=nf90_put_var(ncid, var_id, IOBPwrite, start=(/1/), count =(/np_write/) ) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 3, iret) ! iret=nf90_inq_varid(ncid, "IWBNDGL", var_id) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 4, iret) iret=nf90_put_var(ncid, var_id, ListBound, start=(/1/), count =(/nbBound/) ) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 5, iret) iret=nf90_close(ncid) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 6, iret) END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE INIT_NETCDF_BOUNDARY_OUTPUT USE DATAPOOL IMPLICIT NONE LOGICAL, SAVE :: IsFirstTime = .true. IF (IsFirstTime) THEN IsFirstTime=.FALSE. # ifdef MPI_PARALL_GRID CALL SETUP_BOUNDARY_SCATTER_REDUCE_ARRAY # endif # ifdef MPI_PARALL_GRID IF (myrank .eq. rank_boundary) THEN # endif IF (BOUC_NETCDF_OUT_PARAM .and. (.NOT. allocated(SPPARM_GL))) THEN allocate(SPPARM_GL(8,IWBMNPGL), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_bdcons, allocate error 26') END IF IF (BOUC_NETCDF_OUT_SPECTRA .and. (.NOT. allocated(WBAC_GL))) THEN allocate(WBAC_GL(MSC,MDC,IWBMNPGL), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_bdcons, allocate error 27') END IF # ifdef MPI_PARALL_GRID END IF # endif END IF END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE WRITE_NETCDF_BOUNDARY USE DATAPOOL USE NETCDF implicit none logical, save :: IsInitDone = .FALSE. character(len =256) :: FILE_NAME, PRE_FILE_NAME character (len = *), parameter :: CallFct="WRITE_NETCDF_BOUNDARY" integer iret, ncid, irec_dim, recs_his, var_id integer, save :: ifile = 1 integer LPOS, IP REAL(rkind) :: eTimeDay integer POSITION_BEFORE_POINT, nbTime LPOS=POSITION_BEFORE_POINT(OUT_BOUC % FNAME) IF (OUT_STATION%IDEF.gt.0) THEN WRITE (PRE_FILE_NAME,10) OUT_BOUC % FNAME(1:LPOS),ifile 10 FORMAT (a,'_',i4.4) ELSE WRITE (PRE_FILE_NAME,20) OUT_BOUC % FNAME(1:LPOS) 20 FORMAT (a) ENDIF WRITE (FILE_NAME,30) TRIM(PRE_FILE_NAME) 30 FORMAT (a,'.nc') IF (IsInitDone .eqv. .FALSE.) THEN IsInitDone=.TRUE. nbTime=-1 #ifdef MPI_PARALL_GRID IF (myrank == 0) THEN #endif CALL WRITE_NETCDF_BOUND_HEADERS_1(FILE_NAME, nbTime, np_total, IWBMNPGL, BOUC_NETCDF_OUT_PARAM, BOUC_NETCDF_OUT_SPECTRA) CALL WRITE_NETCDF_BOUND_HEADERS_2(FILE_NAME, np_total, IOBPtotal, IWBMNPGL, IWBNDGL) #ifdef MPI_PARALL_GRID END IF #endif END IF ! ! Getting the needed global arrays ! CALL INIT_NETCDF_BOUNDARY_OUTPUT IF (BOUC_NETCDF_OUT_PARAM .and. LBCWA) THEN CALL REDUCE_BOUNDARY_ARRAY_SPPARM END IF IF (BOUC_NETCDF_OUT_SPECTRA) THEN CALL REDUCE_BOUNDARY_ARRAY_WBAC END IF ! WRITE(STAT%FHNDL,*) 'sum(WBAC)=', sum(WBAC) ! WRITE(STAT%FHNDL,*) 'sum(SPPARM)=', sum(SPPARM) ! WRITE(STAT%FHNDL,*) 'IWBMNP=', IWBMNP ! WRITE(STAT%FHNDL,*) 'IWBMNPGL=', IWBMNPGL ! FLUSH(STAT%FHNDL) ! ! Now writing the boundary data ! IF (myrank .eq. rank_boundary) THEN eTimeDay=MAIN%TMJD iret=nf90_open(TRIM(FILE_NAME), NF90_WRITE, ncid) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 5, iret) iret=nf90_inquire(ncid, unlimitedDimId = irec_dim) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 6, iret) iret=nf90_inquire_dimension(ncid, irec_dim,len = recs_his) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 7, iret) recs_his=recs_his+1 CALL WRITE_NETCDF_TIME(ncid, recs_his, eTimeDay) IF (BOUC_NETCDF_OUT_PARAM .and. LBCWA) THEN iret=nf90_inq_varid(ncid, 'SPPARM', var_id) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 8, iret) IF (NF90_RUNTYPE == NF90_OUTTYPE_BOUC) THEN iret=nf90_put_var(ncid,var_id,SPPARM_GL, start=(/1,1,recs_his/), count = (/8, IWBMNPGL,1/)) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 9, iret) ELSE iret=nf90_put_var(ncid,var_id,SNGL(SPPARM_GL), start=(/1,1,recs_his/), count = (/8, IWBMNPGL,1/)) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 10, iret) ENDIF END IF IF (BOUC_NETCDF_OUT_SPECTRA) THEN ! WRITE(STAT%FHNDL,*) 'sum(WBAC_GL)=', sum(WBAC_GL) ! WRITE(STAT%FHNDL,*) 'sum(SPPARM_GL)=', sum(SPPARM_GL) ! FLUSH(STAT%FHNDL) ! iret=nf90_inq_varid(ncid, 'WBAC', var_id) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 11, iret) IF (NF90_RUNTYPE == NF90_OUTTYPE_BOUC) THEN iret=nf90_put_var(ncid,var_id,WBAC_GL, start=(/1,1,1,recs_his/), count = (/MSC,MDC, IWBMNPGL,1/)) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 12, iret) ELSE iret=nf90_put_var(ncid,var_id,SNGL(WBAC_GL), start=(/1,1,1,recs_his/), count = (/MSC,MDC, IWBMNPGL,1/)) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 13, iret) ENDIF END IF iret=nf90_close(ncid) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 14, iret) END IF IF (OUT_BOUC % IDEF.gt.0) THEN IF (recs_his .eq. OUT_BOUC % IDEF) THEN ifile=ifile+1 IsInitDone = .FALSE. ENDIF ENDIF END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE READ_NETCDF_BOUNDARY_WBAC_SINGLE(IFILE, IT) USE DATAPOOL USE NETCDF IMPLICIT NONE integer, intent(in) :: IFILE, IT character (len = *), parameter :: CallFct="READ_NETCDF_BOUNDARY_WBAC_SINGLE" integer ncid, var_id ISTAT = NF90_OPEN(BOUC_NETCDF_FILE_NAMES(IFILE), NF90_NOWRITE, ncid) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 1, ISTAT) ISTAT = nf90_inq_varid(ncid, 'WBAC', var_id) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 2, ISTAT) ISTAT = NF90_GET_VAR(ncid, var_id, WBAC_GL, start=(/1,1,1,IT/), count = (/MSC,MDC, IWBMNPGL,1/)) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 3, ISTAT) ISTAT = NF90_CLOSE(ncid) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 4, ISTAT) END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE READ_NETCDF_BOUNDARY_WBAC(WBACOUT) USE DATAPOOL IMPLICIT NONE REAL(rkind), INTENT(OUT) :: WBACOUT(MSC,MDC,IWBMNP) integer IFILE, IT integer IP CALL COMPUTE_IFILE_IT_BOUND(IFILE, IT) # ifdef MPI_PARALL_GRID IF (MULTIPLE_IN_GRID) THEN CALL READ_NETCDF_BOUNDARY_WBAC_SINGLE(IFILE, IT) DO IP=1,IWBMNP WBACOUT(:,:,IP)=WBAC_GL(:,:,Indexes_boundary(IP)) END DO ELSE IF (myrank .eq. rank_boundary) THEN CALL READ_NETCDF_BOUNDARY_WBAC_SINGLE(IFILE, IT) END IF CALL SCATTER_BOUNDARY_ARRAY_WBAC(WBACOUT) END IF # else CALL READ_NETCDF_BOUNDARY_WBAC_SINGLE(IFILE, IT) WBACOUT=WBAC_GL # endif END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE READ_NETCDF_BOUNDARY_SPPARM_SINGLE(IFILE, IT) USE DATAPOOL USE NETCDF IMPLICIT NONE integer, intent(in) :: IFILE, IT character (len = *), parameter :: CallFct="READ_NETCDF_BOUNDARY_SPPARM_SINGLE" integer ncid, var_id ! Print *, 'IWBMNPGL=', IWBMNPGL ! Print *, 'IT=', IT ISTAT = NF90_OPEN(BOUC_NETCDF_FILE_NAMES(IFILE), NF90_NOWRITE, ncid) ! Print *, 'step 1' CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 1, ISTAT) ! Print *, 'step 1' ISTAT = nf90_inq_varid(ncid, 'SPPARM', var_id) ! Print *, 'step 2' CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 2, ISTAT) ! Print *, 'step 3' ! Print *, 'allocated(SPPARM_GL)=', allocated(SPPARM_GL) ISTAT = NF90_GET_VAR(ncid, var_id, SPPARM_GL, start=(/1,1,IT/), count = (/8, IWBMNPGL,1/)) ! Print *, 'step 4' CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 3, ISTAT) ! Print *, 'step 5' ISTAT = NF90_CLOSE(ncid) ! Print *, 'step 6' CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 4, ISTAT) ! Print *, 'step 7' END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE READ_NETCDF_BOUNDARY_SPPARM USE DATAPOOL IMPLICIT NONE INTEGER IFILE, IT integer IP CALL COMPUTE_IFILE_IT_BOUND(IFILE, IT) # ifdef MPI_PARALL_GRID IF (MULTIPLE_IN_GRID) THEN CALL READ_NETCDF_BOUNDARY_SPPARM_SINGLE(IFILE, IT) DO IP=1,IWBMNP SPPARM(:,IP)=SPPARM_GL(:,Indexes_boundary(IP)) END DO ELSE IF (myrank .eq. rank_boundary) THEN CALL READ_NETCDF_BOUNDARY_SPPARM_SINGLE(IFILE, IT) END IF CALL SCATTER_BOUNDARY_ARRAY_SPPARM END IF # else CALL READ_NETCDF_BOUNDARY_SPPARM_SINGLE(IFILE, IT) SPPARM=SPPARM_GL # endif END SUBROUTINE !********************************************************************** !* accepted input files for NETCDF_IN_FILE in wwminput.nml * !* NETCDF_IN_FILE = 'FileIn.nc' if FileIn.nc exists * !* or FileIn_0001.nc, ...., FileIn_0002.nc if they exist * !********************************************************************** SUBROUTINE INIT_NETCDF_BOUNDARY_WWM USE NETCDF USE DATAPOOL IMPLICIT NONE integer POSITION_BEFORE_POINT, LPOS character(len=140) FILE_NAME logical LFLIVE INTEGER iFile, jFile, iTime, nbtime_mjd INTEGER dimids(2), varid, fid real(rkind) :: ConvertToDay real(rkind) :: eTimeStart, eTime, dtbound real(rkind), allocatable :: ListTime_mjd(:) character (len=100) :: eStrUnitTime integer idx character (len = *), parameter :: CallFct="INIT_NETCDF_BOUNDAY_WWM" IF (TRIM(NETCDF_IN_FILE) .eq. "unset") THEN CALL WWM_ABORT('NETCDF_IN_FILE must be set for running') END IF INQUIRE( FILE = TRIM(NETCDF_IN_FILE), EXIST = LFLIVE ) ! ! First determination of the file names ! IF (LFLIVE ) THEN NUMBER_BOUC_NETCDF_FILE=1 ALLOCATE(BOUC_NETCDF_FILE_NAMES(1), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_bdcons, allocate error 26') BOUC_NETCDF_FILE_NAMES(1)=TRIM(NETCDF_IN_FILE) ELSE LPOS=POSITION_BEFORE_POINT(OUT_BOUC % FNAME) iFile=0 DO jFile=iFile+1 WRITE (FILE_NAME,10) NETCDF_IN_FILE(1:LPOS),jFile INQUIRE( FILE = TRIM(FILE_NAME), EXIST = LFLIVE ) IF (LFLIVE .eqv. .FALSE.) THEN EXIT END IF iFile=jFile END DO NUMBER_BOUC_NETCDF_FILE=iFile IF (NUMBER_BOUC_NETCDF_FILE .eq. 0) THEN Print *, 'Error in INIT_NETCDF_BOUNDARY_WWM' Print *, 'NETCDF_IN_FILE=', TRIM(NETCDF_IN_FILE) Print *, 'FILE_NAME=', TRIM(FILE_NAME) CALL WWM_ABORT('We did not find the needed boundary files') END IF ALLOCATE(BOUC_NETCDF_FILE_NAMES(NUMBER_BOUC_NETCDF_FILE), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_bdcons, allocate error 27') DO iFile=1,NUMBER_BOUC_NETCDF_FILE WRITE (FILE_NAME,10) NETCDF_IN_FILE(1:LPOS),jFile BOUC_NETCDF_FILE_NAMES(iFile)=TRIM(FILE_NAME) END DO END IF ! ! next reading the times ! BOUND_NB_TIME=0 DO iFile=1,NUMBER_BOUC_NETCDF_FILE ISTAT = nf90_open(TRIM(BOUC_NETCDF_FILE_NAMES(iFile)), nf90_nowrite, fid) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 1, ISTAT) ISTAT = nf90_inq_varid(fid, "ocean_time", varid) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 5, ISTAT) ISTAT = nf90_inquire_variable(fid, varid, dimids=dimids) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 8, ISTAT) ISTAT = nf90_inquire_dimension(fid, dimids(1), len=nbtime_mjd) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 9, ISTAT) ISTAT = nf90_close(fid) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 11, ISTAT) BOUND_NB_TIME = BOUND_NB_TIME + nbtime_mjd END DO ALLOCATE(BOUND_LIST_IFILE(BOUND_NB_TIME), BOUND_LIST_IT(BOUND_NB_TIME), BOUND_LIST_TIME(BOUND_NB_TIME), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_bdcons, allocate error 28') idx=0 DO iFile=1,NUMBER_BOUC_NETCDF_FILE ISTAT = nf90_open(TRIM(BOUC_NETCDF_FILE_NAMES(iFile)), nf90_nowrite, fid) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 1, ISTAT) ISTAT = nf90_inq_varid(fid, "ocean_time", varid) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 5, ISTAT) ISTAT = nf90_get_att(fid, varid, "units", eStrUnitTime) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 7, ISTAT) CALL CF_EXTRACT_TIME(eStrUnitTime, ConvertToDay, eTimeStart) ISTAT = nf90_inquire_variable(fid, varid, dimids=dimids) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 8, ISTAT) ISTAT = nf90_inquire_dimension(fid, dimids(1), len=nbtime_mjd) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 9, ISTAT) allocate(ListTime_mjd(nbtime_mjd), stat=istat) IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 48') ISTAT = nf90_get_var(fid, varid, ListTime_mjd) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 10, ISTAT) ISTAT = nf90_close(fid) CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 11, ISTAT) DO iTime=1,nbtime_mjd idx=idx+1 eTime=ListTime_mjd(iTime)*ConvertToDay + eTimeStart BOUND_LIST_IFILE(idx) = iFile BOUND_LIST_IT(idx) = iTime BOUND_LIST_TIME(idx) = eTime END DO DEALLOCATE(ListTime_mjd) END DO DTBOUND = (BOUND_LIST_TIME(2) - BOUND_LIST_TIME(1))*DAY2SEC SEBO%DELT = DTBOUND SEBO%BMJD = BOUND_LIST_TIME(1) SEBO%EMJD = BOUND_LIST_TIME(BOUND_NB_TIME) RETURN 10 FORMAT (a,'_',i4.4) END SUBROUTINE #endif !********************************************************************** !* * !********************************************************************** SUBROUTINE COMPUTE_IFILE_IT_BOUND(IFILE, IT) USE DATAPOOL IMPLICIT NONE integer, intent(out) :: IFILE, IT REAL(rkind) :: DeltaTime, DeltaT integer iTime DeltaTime=BOUND_LIST_TIME(2) - BOUND_LIST_TIME(1) DeltaT=(MAIN%TMJD - BOUND_LIST_TIME(1))/DeltaTime iTime=NINT(DeltaT) + 1 IFILE=BOUND_LIST_IFILE(iTime) IT=BOUND_LIST_IT(iTime) END SUBROUTINE !********************************************************************** !* * !********************************************************************** SUBROUTINE EXPORT_BOUC_WW3_FORMAT USE DATAPOOL IMPLICIT NONE CHARACTER(LEN=32), PARAMETER :: IDSTRBC='WAVEWATCH III BOUNDARY DATA FILE' CHARACTER(LEN=10), PARAMETER :: VERBPTBC = 'III 1.03 ' integer nbNeumann, nbDirichlet, nbUnknown integer IP LOGICAL, SAVE :: IsFirst = .TRUE. REAL, allocatable :: XBPI(:), YBPI(:), RDBPI(:,:) INTEGER, allocatable :: IPBPI(:,:) REAL(rkind) :: WVK,WVCG,WVKDEP,WVN,WVC,SPSIGLOC real, allocatable :: ABPIO(:) REAL(rkind) :: eCLATS, eCG, DEPLOC, eVal REAL XFR, eTH, FREQ1 INTEGER NBI, idx, IB, I, J, NSPEC_out, IK, ITH, ISP, NK, NTH, IPglob INTEGER TheOut INTEGER TIME2(2) nbDirichlet=0 nbNeumann=0 nbUnknown=0 DO IP=1,np_total IF (IOBPtotal(IP) == 2) THEN nbDirichlet=nbDirichlet+1 END IF IF (IOBPtotal(IP) == 3) THEN nbNeumann=nbNeumann+1 END IF IF (IOBPtotal(IP) == 4) THEN nbUnknown=nbUnknown+1 END IF END DO IF ((nbNeumann .gt. 0).or.(nbUnknown .gt. 0)) THEN # ifdef MPI_PARALL_GRID IF (myrank == 0) THEN # endif IF (IsFirst .eqv. .TRUE.) THEN Print *, 'nbDirichlet=', nbDirichlet Print *, 'nbNeumann=', nbNeumann Print *, 'nbUnknown=', nbUnknown Print *, 'Those points will be put to 0' END IF # ifdef MPI_PARALL_GRID END IF # endif END IF NBI = nbDirichlet IF (NBI .ne. IWBMNPGL) THEN CALL WWM_ABORT('Code inconsistency error') END IF allocate(XBPI(NBI), YBPI(NBI), IPBPI(NBI,4), RDBPI(NBI,4), stat=istat) IF (istat/=0) CALL WWM_ABORT('Error allocate XBPI/YBPI') idx=0 DO IP=1,NP_TOTAL IF (IOBPtotal(IP) == 2) THEN idx=idx+1 XBPI(idx)=MySNGL(XPtotal(IP)) YBPI(idx)=MySNGL(YPtotal(IP)) END IF END DO DO IB=1,NBI IPBPI(IB,1)=IB IPBPI(IB,2:4)=0 END DO RDBPI(:,1) = 1 RDBPI(:,2:4) = 0 CALL REDUCE_BOUNDARY_ARRAY_WBAC # ifdef MPI_PARALL_GRID IF (myrank == 0) THEN # endif TheOut = FHNDL_EXPORT_BOUC_WW3 ! Print *, 'IsFirst=', IsFirst NK = MSC NTH = MDC XFR = MySNGL(SFAC) FREQ1 = MySNGL(FR(1)) eTH = MySNGL(SPDIR(1)) IF (IsFirst .eqv. .TRUE.) THEN ! Print *, 'Before open, case 1' OPEN(TheOut, FILE='nest.ww3', FORM='UNFORMATTED', status='new', action='write') WRITE(TheOut) IDSTRBC, VERBPTBC, NK, NTH, XFR, FREQ1, eTH, NBI WRITE(TheOut) (XBPI(I),I=1,NBI), (YBPI(I),I=1,NBI), & ((IPBPI(I,J),I=1,NBI),J=1,4), & ((RDBPI(I,J),I=1,NBI),J=1,4) ELSE ! Print *, 'Before open, case 2' OPEN(TheOut, FILE='nest.ww3', FORM='UNFORMATTED', status='old', position='append', action='write') END IF ! CALL COMPUTE_TFN(TIME2) NSPEC_out = NK*NTH allocate(ABPIO(NSPEC_out), stat=istat) IF (istat/=0) CALL WWM_ABORT('Error allocate ABPIO') WRITE(TheOut) TIME2, NBI DO IB=1,NBI IPglob=IWBNDGL(IB) DEPLOC = MAX(DMIN,DEPtotal(IPglob)) IF (LSPHE) THEN eCLATS = COS(DEGRAD*YPtotal(IPglob)) ELSE eCLATS = 1 END IF DO IK=1,NK CALL ALL_FROM_TABLE(SPSIGLOC,DEPLOC,WVK,WVCG,WVKDEP,WVN,WVC) eCG = WVCG DO ITH=1,NTH ISP = ITH + (IK-1)*NTH eVal= WBAC_GL(IK,ITH,IB)*eCG/eCLATS ABPIO(ISP) = MySNGL(eVal) END DO END DO WRITE(TheOut) ABPIO END DO deallocate(ABPIO) CLOSE(TheOut) # ifdef MPI_PARALL_GRID END IF # endif deallocate(XBPI, YBPI, IPBPI, RDBPI) IsFirst=.FALSE. END SUBROUTINE !********************************************************************** !* * !**********************************************************************