#include "cppdefs.h" ! !************************************************************************ SUBROUTINE swan_reader (ng, first, tile) !************************************************************************ ! #if defined INWAVE_MODEL & defined INWAVE_SWAN_COUPLING ! !svn $Id: swan_reader.F 1336 2008-01-24 02:45:56Z jcwarner $ ! LAST CHANGE: mai 12/28/2010 ! !======================================================================! ! ! ! This routine reads the output spectra from swan.................... ! ! ! !======================================================================! ! USE inwave_iounits USE mod_inwave_swan USE mod_inwave_vars USE mod_parallel USE mod_scalars USE Tr3dbc_mod # ifdef DISTRIBUTE USE mp_exchange_mod # endif ! implicit none ! ! Imported variable declarations. ! integer, intent(in) :: ng, first, tile ! ! Local variable declarations. ! character(6) :: rtext character(72) :: rline real(r8) :: factor, exc, m0, frep, cff, DTTIME integer :: nfreq, ndir, switch integer :: i, j, d, inp integer :: nfreqsw, ndirsw integer :: minpos, readerr integer :: INTTIM(6) integer :: LBi, UBi, LBj, UBj integer,allocatable :: posang(:) real(r8), pointer :: temp(:), findline(:) real(r8), pointer :: tempA(:,:) real(r8), pointer :: sdsw1(:,:) real(r8), pointer :: theta1(:) real(r8), parameter :: rad = pi/180.0_r8 real(r8), parameter :: orad = 180.0_r8/pi # include "set_bounds.h" ! ! Determine array lower and upper bounds in the I- and J-directions. ! LBi=BOUNDS(ng)%LBi(tile) UBi=BOUNDS(ng)%UBi(tile) LBj=BOUNDS(ng)%LBj(tile) UBj=BOUNDS(ng)%UBj(tile) ! !----------------------------------------------------------------------- ! Open the spectral file !----------------------------------------------------------------------- ! readerr=0 inp=44 open(inp,file=IWSWNname(ng),form='formatted',status='old') switch=0 ! !----------------------------------------------------------------------- ! Read file until RFREQ or AFREQ is found !----------------------------------------------------------------------- ! DO while (switch==0) read(inp,'(a)')rtext IF (rtext == 'RFREQ ') THEN switch = 1 ELSEIF (rtext == 'AFREQ ') THEN switch = 2 END IF ENDDO ! !----------------------------------------------------------------------- ! Read nfreq and f !----------------------------------------------------------------------- ! read(inp,*)nfreqsw IF (first.eq.1) allocate(WAVES(ng)%fsw(nfreqsw)) DO i=1,nfreqsw read(inp,*) WAVES(ng)%fsw(i) ENDDO WAVES(ng)%nfreqsw=nfreqsw ! !----------------------------------------------------------------------- ! Convert to absolute frequencies !----------------------------------------------------------------------- ! IF (switch == 1) THEN WAVES(ng)%fsw = WAVES(ng)%fsw ELSE WAVES(ng)%fsw = WAVES(ng)%fsw ENDIF ! !----------------------------------------------------------------------- ! Read CDIR or NDIR !----------------------------------------------------------------------- ! read(inp,'(a)')rtext IF (rtext == 'NDIR ') THEN switch = 1 ELSEIF (rtext == 'CDIR ') THEN switch = 2 ELSE stop ENDIF ! !----------------------------------------------------------------------- ! Read ndir, theta !----------------------------------------------------------------------- ! read(inp,*) ndirsw allocate(theta1(ndirsw)) IF (first.eq.1) allocate(WAVES(ng)%theta(ndirsw)) DO i=1,ndirsw read(inp,*)theta1(i) IF(theta1(i).lt.0.0_r8)THEN theta1(i)=(360.0_r8+theta1(i)) ENDIF theta1(i)=theta1(i)*rad ENDDO WAVES(ng)%ndirsw=ndirsw WAVES(ng)%dang=abs(theta1(2)-theta1(1)) minpos=minloc(theta1, DIM = 1) allocate(posang(ndirsw)) ! If swan was run on a full circle or just a sector. IF (INT(2.0_r8*pi/WAVES(ng)%dang).eq.ndirsw) THEN WAVES(ng)%Swancircle=1 ELSE WAVES(ng)%Swancircle=0 END IF IF((theta1(minpos+1)-theta1(minpos)).lt. & & (theta1(minpos-1)-theta1(minpos)))THEN DO i=1,ndirsw posang(i)=minpos+(i-1) IF(posang(i).gt.ndirsw)THEN posang(i)=posang(i)-ndirsw ENDIF ENDDO ELSE DO i=1,ndirsw posang(i)=minpos-(i-1) IF(posang(i).lt.1)THEN posang(i)=ndirsw+posang(i) ENDIF ENDDO ENDIF ! !----------------------------------------------------------------------- ! Find exception value !----------------------------------------------------------------------- ! i=0 DO WHILE (i==0) READ (inp,'(a)',ERR=215,END=240) rline IF(rline(41:49).eq.'exception') THEN read(rline(1:15),*) exc i=1 END IF END DO ! !----------------------------------------------------------------------- ! Get time !----------------------------------------------------------------------- ! i=0 DO WHILE (i==0) READ (inp,'(a)',ERR=215,END=240) rline IF(rline(41:44).eq.'date') THEN read(rline(1:15),*) WAVES(ng)%SpecTimeIso(1) read(rline(1:4),'(I4)') INTTIM(1) read(rline(5:6),'(I4)') INTTIM(2) read(rline(7:8),'(I4)') INTTIM(3) read(rline(10:11),'(I4)') INTTIM(4) read(rline(12:13),'(I4)') INTTIM(5) read(rline(14:15),'(I4)') INTTIM(6) i=1 WAVES(ng)%SpecTime(1)=DTTIME(INTTIM) ! Set dt SpecTime to a large valule in case the spec2d file ! has only 1 time instance. IF (first.eq.1) THEN WAVES(ng)%SpecTime(2)=WAVES(ng)%SpecTime(1) WAVES(ng)%SpecTimedt=99999999.0_r8 END IF END IF END DO ! !----------------------------------------------------------------------- ! Find FACTOR keyword !----------------------------------------------------------------------- ! i=0 DO WHILE (i==0) read(inp,'(a)')rtext IF (rtext == 'FACTOR') THEN i=1 ELSEIF (rtext == 'ZERO ') THEN IF (Master) THEN WRITE(stdout,*) 'Zero energy density input for this point' END IF stop ELSEIF (rtext == 'NODATA') THEN IF (Master) THEN WRITE(stdout,*) 'SWAN file has no data for this point' END IF stop ENDIF END DO read(inp,*)factor ! !----------------------------------------------------------------------- ! Read S_array !----------------------------------------------------------------------- ! allocate(sdsw1(nfreqsw,ndirsw)) IF (first.eq.1) allocate(WAVES(ng)%SDSW(nfreqsw,ndirsw)) DO i=1,nfreqsw read(inp,*)sdsw1(i,:) END DO where (sdsw1 == exc) sdsw1=0.0_r8 ! !----------------------------------------------------------------------- ! If first, then read the next time and close file for now. ! If not first, then read until the next times step and get data. !----------------------------------------------------------------------- ! j=0 DO WHILE (j==0) READ (inp,'(a)',ERR=215,END=245) rline IF(rline(41:44).eq.'date') THEN read(rline(1:15),*) WAVES(ng)%SpecTimeIso(2) read(rline(1:4),'(I4)') INTTIM(1) read(rline(5:6),'(I4)') INTTIM(2) read(rline(7:8),'(I4)') INTTIM(3) read(rline(10:11),'(I4)') INTTIM(4) read(rline(12:13),'(I4)') INTTIM(5) read(rline(14:15),'(I4)') INTTIM(6) WAVES(ng)%SpecTime(2)=DTTIME(INTTIM) IF (first.eq.1) THEN ! Compute the dt SpecTime if we do get a second time stamp. WAVES(ng)%SpecTimedt=WAVES(ng)%SpecTime(2)- & & WAVES(ng)%SpecTime(1) j=1 IF (Master) THEN WRITE(stdout,*) 'Using 2D spec data at ', & & WAVES(ng)%SpecTimeIso(1) END IF ELSE ! If not first, then keep reading until we get the next time stamp. IF ((time(ng)-dstart*24.0_r8*3600.0_r8+dt(ng)).gt. & & (WAVES(ng)%SpecTime(2)-WAVES(ng)%SpecTime(1))) THEN read(inp,'(a)')rtext IF (rtext == 'ZERO ') THEN IF (Master) THEN WRITE(stdout,*) 'Zero energy density input' END IF stop ELSEIF (rtext == 'NODATA') THEN IF (Master) THEN WRITE(stdout,*) 'SWAN file has no data' END IF stop ENDIF read(inp,*)factor DO i=1,nfreqsw read(inp,*)sdsw1(i,:) END DO where (sdsw1 == exc) sdsw1=0.0_r8 IF (Master) THEN WRITE(stdout,*) 'Using 2D spec data at ', & & WAVES(ng)%SpecTimeIso(2) END IF ELSE j=1 END IF END IF END IF END DO 245 CONTINUE close(inp) ! !----------------------------------------------------------------------- ! Convert to m2/Hz/rad !----------------------------------------------------------------------- ! DO j=1,ndirsw DO i=1,nfreqsw sdsw1(i,j)=sdsw1(i,j)*factor*orad END DO END DO ! !----------------------------------------------------------------------- ! Convert from energy density to variance density and order !----------------------------------------------------------------------- ! DO j=1,ndirsw WAVES(ng)%theta(j)=theta1(posang(j)) DO i=1,nfreqsw WAVES(ng)%SDSW(i,j)=sdsw1(i,posang(j))/(rho0*g) ENDDO ENDDO deallocate(sdsw1) deallocate(theta1) deallocate(posang) ! !----------------------------------------------------------------------- ! Integrate the directional spectra over the directions !----------------------------------------------------------------------- ! IF (first.eq.1) allocate(WAVES(ng)%Sf(nfreqsw)) WAVES(ng)%Sf = sum(WAVES(ng)%SDSW, DIM = 2)*WAVES(ng)%dang ! !----------------------------------------------------------------------- ! Find main wave direction !----------------------------------------------------------------------- ! allocate (temp(ndirsw)) ! !----------------------------------------------------------------------- ! Integrate the directional spectra over the frequencies !----------------------------------------------------------------------- ! temp=sum(WAVES(ng)%SDSW, DIM = 1) i=maxval(maxloc(temp)) WAVES(ng)%mainang=WAVES(ng)%theta(i) IF (Master) THEN WRITE(stdout,*) 'Peak dir from ', WAVES(ng)%mainang*180.0_r8/pi END IF deallocate(temp) ! !----------------------------------------------------------------------- ! Calculate zero-order moment !----------------------------------------------------------------------- ! allocate (temp(nfreqsw+1)) temp(1)=0.0_r8 DO i=2,nfreqsw temp(i)=0.5_r8*WAVES(ng)%fsw(i-1)+ & & 0.5_r8*WAVES(ng)%fsw(i) END DO temp(nfreqsw+1)=WAVES(ng)%fsw(nfreqsw) m0=0.0_r8 DO i=1,nfreqsw m0=m0+WAVES(ng)%Sf(i)*(temp(i+1)-temp(i)) END DO deallocate (temp) WAVES(ng)%hm0gew=4.004_r8*sqrt(m0) IF (Master) THEN WRITE(stdout,*) 'Waves hm0 is ', WAVES(ng)%hm0gew END IF ! !----------------------------------------------------------------------- ! Compute the representative period !----------------------------------------------------------------------- ! call tpdcalc(WAVES(ng)%Sf,WAVES(ng)%fsw,frep) ! WAVEG(ng)%Trep=1.0_r8/frep WAVES(ng)%fp=frep IF (Master) THEN WRITE(stdout,*) 'Representative period is ', WAVEG(ng)%Trep END IF ! GO TO 250 215 IF (Master) THEN WRITE (stdout,*) 'swan_reader.F Error reading swan 2d spec file' END IF ! exit_flag=4 240 IF (Master) THEN WRITE(stdout,*) 'END of 2d spec file' END IF 250 CONTINUE #endif RETURN END SUBROUTINE swan_reader !******************************************************************* FUNCTION DTTIME (INTTIM) ! * ! Borrowed from SWAN ocpmix.F DTTIME * !******************************************************************* ! USE mod_scalars IMPLICIT NONE real(r8) :: DTTIME ! INTEGER INTTIM(6) ! INTEGER IDYMON(12), IYEAR, IYRM1, IDNOW, I, II ! LOGICAL LEAPYR, LOGREF ! DATA IDYMON /31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/ ! IYEAR = INTTIM(1) IYRM1 = IYEAR-1 LEAPYR=(MOD(IYEAR,4).EQ.0.AND.MOD(IYEAR,100).NE.0).OR. & & MOD(IYEAR,400).EQ.0 IDNOW=0 DO I = 1,INTTIM(2)-1 IDNOW=IDNOW+IDYMON(I) END DO IDNOW=IDNOW+INTTIM(3) IF (LEAPYR.AND.INTTIM(2).GT.2) IDNOW=IDNOW+1 IDNOW = IDNOW + IYEAR*365 + IYRM1/4 - IYRM1/100 + IYRM1/400 + 1 IF (IYEAR.EQ.0) IDNOW=IDNOW-1 ! IF (.NOT.LOGREF) THEN ! REFDAY = IDNOW ! LOGREF = .TRUE. ! DTTIME = 0. ! ELSE ! DTTIME = REAL(IDNOW-REFDAY) * 24.*3600. DTTIME = REAL(IDNOW) * 24.*3600. ! ENDIF DTTIME = DTTIME + 3600.*REAL(INTTIM(4)) + 60.*REAL(INTTIM(5)) + & & REAL(INTTIM(6)) RETURN END