PROGRAM HYCOM_TIDAL_RI2AP IMPLICIT NONE !DAN====================================================================== C C hycom_tidal_ri2ap - C C Usage: hycom_tidal_ri2ap ftideRIin.a ftideAPout.a [grid.a] C Usage: hycom_tidal_ri2ap180 ftideRIin.a ftideAPout.a [grid.a] C C Purpose: To convert the Real Field, Imag Field output from the C hycom_tidal_foreman program to C Amplitude and Phase fields for each Tidal Mode C The values for the Phase lie between 0 and 360 degrees C or between -180 and 180 degrees C C ftideRIin.a contains two fields for the real and imaginary parts C of each tidal mode generated by hycom_tidal_foreman. C it follows the OSU (TPXO) convention for the Imaginary component: C phase = atan2(-Im,Re). C C ftideAPout.a will contain two fields for the Amplitude and Phase of C each tidal mode generated by hycom_tidal_foreman. C for ri2ap the phase will be in the range 0 to 360 and C for ri2ap180 the phase will be in the range -180 to 180. C C grid.a is a hycom grid file, default regional.grid.a. Note that C the corresponding grid.b must also exist. C this version is for "serial" systems. C C Daniel R. Moore (QNA) and Alan J. Wallcraft (NRL), September 2012. C C========================================================================= REAL*4, ALLOCATABLE :: TideR(:,:),TideI(:,:),TideA(:,:),TideP(:,:) REAL*4 :: PAD(4096),SPVAL,Pi, & Tmax,Tmin,Pmin,Pmax,Phase DATA Pi/3.141592654/ CHARACTER*240 CARG PARAMETER (SPVAL=2.0**100) C LOGICAL L180 INTEGER IDM,JDM,I,J,K,L,NPAD,NRECL,IP_test,JP_test,NARG INTEGER IOS,IREC,IGG,IARGC CHARACTER*240 CFILE_IN,CFILE_OUT,CFILE_GRID,CFILEB CHARACTER*6 CVARIN C C READ ARGUMENTS. C CALL GETARG(0,CARG) L = LEN_TRIM(CARG) * WRITE(6,"(4a)") TRIM(CARG),' "',CARG(L-2:L),'"' L180 = CARG(L-2:L).EQ."180" NARG = IARGC() C IF((NARG-2)*(NARG-3)*(NARG-5).NE.0)THEN WRITE(6,*)'2 or 3 arguments expected!, got ',NARG WRITE(6,*) & 'Usage: hycom_RI2AP ftideRIin.a ftideAPout.a [grid.a]' CALL EXIT(1) ENDIF c c First 2 arguments are common to both possible argument numbers: 2 or 3 c CALL GETARG(1,CFILE_IN) CALL GETARG(2,CFILE_OUT) C C Process Variant in argument numbers (NARG = 3 !) C C NARG = 3 regional.grid.a format file (and the associated .b file !) C specifying the geometry of the region C IF(NARG.EQ.2)THEN CFILE_GRID='regional.grid.a' ELSE CALL GETARG(3,CFILE_GRID) ENDIF C C Undocumented debug option C Usage: hycom_tidal_ri2ap ftideRIin.a ftideAPout.a grid.a itest jtest C Usage: hycom_tidal_ri2ap180 ftideRIin.a ftideAPout.a grid.a itest jtest C IF(NARG.EQ.5)THEN CALL GETARG(4,CARG) READ(CARG,*) IP_test CALL GETARG(5,CARG) READ(CARG,*) JP_test ELSE IP_test = 0 JP_test = 0 ENDIF c------------------------------------------------------------- c Test arguments c print *,'Input File = ',TRIM(CFILE_IN) print *,'grid.a = ',TRIM(CFILE_GRID) print *,'Output File = ',TRIM(CFILE_OUT) C READ(5,*)IGG C----------------------------------------------------------------------------- C Set up to output the tidal modes conversion at a specific point: C (IP_trest,JP_Test) C Skip the generation of this output if NTest_Out <= 0 ! C IF(IP_test.gt.0)write(6,*)'Tidal Mode Conversion to be ', + 'printed at:(',IP_test,',',JP_test,')' c---------------------------------------------------------------- C Get IDM and JDM from regional.grid.n C---------------------------------------------------------------- C C GET IDM,JDM FROM regional.grid.b. C CFILEB = CFILE_GRID(1:LEN_TRIM(CFILE_GRID)-1) // 'b' WRITE(6,'(a,a)')' Grid data file = ',TRIM(CFILEB) C OPEN(UNIT=11,FILE=CFILEB,FORM='FORMATTED', & STATUS='OLD',ACTION='READ') C READ( 11,*) IDM,CVARIN IF (CVARIN.NE.'idm ') THEN WRITE(6,*) 'hycom_RI2AP: bad region.grid.b file ', & CFILEB(1:LEN_TRIM(CFILEB)) CALL EXIT(2) ENDIF READ( 11,*) JDM,CVARIN IF (CVARIN.NE.'jdm ') THEN WRITE(6,*) 'hycom_RI2AP: bad region.grid.b file ', & CFILEB(1:LEN_TRIM(CFILEB)) CALL EXIT(2) ENDIF C CLOSE(UNIT=11) write(6,116)IDM,JDM,CFILE_GRID(1:LEN_TRIM(CFILE_GRID)-1) // 'b' 116 format ( & i5,4x,'''idm '' = longitudinal array size'/ & i5,4x,'''jdm '' = latitudinal array size'/ & a70) C-------------------------------------------------------------------------------- C Allocate Arrays to hold TideR,TideI,TideA,TideP c ALLOCATE( TideR(IDM,JDM), STAT=IOS ) IF (IOS.NE.0) THEN WRITE(6,*) 'Error in hycom_RI2AP: could not allocate ', + IDM*JDM,' words for TideR' CALL EXIT(2) ENDIF write(6,*)'Array TideR(IDM,JDM) allocated' ALLOCATE( TideI(IDM,JDM), STAT=IOS ) IF (IOS.NE.0) THEN WRITE(6,*) 'Error in hycom_RI2AP: could not allocate ', + IDM*JDM,' words for TideI' CALL EXIT(2) ENDIF write(6,*)'Array TideI(IDM,JDM) allocated' ALLOCATE( TideA(IDM,JDM), STAT=IOS ) IF (IOS.NE.0) THEN WRITE(6,*) 'Error in hycom_RI2AP: could not allocate ', + IDM*JDM,' words for TideA' CALL EXIT(2) ENDIF write(6,*)'Array TideA(IDM,JDM) allocated' ALLOCATE( TideP(IDM,JDM), STAT=IOS ) IF (IOS.NE.0) THEN WRITE(6,*) 'Error in hycom_RI2AP: could not allocate ', + IDM*JDM,' words for TideP' CALL EXIT(2) ENDIF write(6,*)'Array TideP(IDM,JDM) allocated' C---------------------------------------------------------------- C Determine Padding to read in a Field as a single record. C NPAD = 4096 - MOD(IDM*JDM,4096) IF (NPAD.EQ.4096) THEN NPAD = 0 INQUIRE(IOLENGTH=NRECL) TideR ELSE INQUIRE(IOLENGTH=NRECL) TideR,PAD(1:NPAD) ENDIF write(6,'(a,i5,i9)') 'npad,nrecl =',npad,nrecl C C read(5,*)IGG #ifdef CRAY #ifdef t3e IF (MOD(NRECL,4096).EQ.0) THEN WRITE(CASN,8000) NRECL/4096 8000 FORMAT('-F cachea:',I4.4,':1:0') IU8 = 11 CALL ASNUNIT(IU8,CASN,IOS8) IF (IOS8.NE.0) THEN write(6,*) 'Error: can''t asnunit 11' write(6,*) 'ios = ',ios8 write(6,*) 'casn = ',casn CALL EXIT(5) ENDIF IU8 = 21 CALL ASNUNIT(IU8,CASN,IOS8) IF (IOS8.NE.0) THEN write(6,*) 'Error: can''t asnunit 21' write(6,*) 'ios = ',ios8 write(6,*) 'casn = ',casn CALL EXIT(5) ENDIF ENDIF #else CALL ASNUNIT(11,'-F syscall -N ieee',IOS) IF (IOS.NE.0) THEN write(6,*) 'Error: can''t asnunit 11' write(6,*) 'ios = ',ios CALL EXIT(5) ENDIF CALL ASNUNIT(21,'-F syscall -N ieee',IOS) IF (IOS.NE.0) THEN write(6,*) 'Error: can''t asnunit 21' write(6,*) 'ios = ',ios CALL EXIT(5) ENDIF #endif #endif C READ(5,*)IGG C====================================================================== C Open Input File C OPEN(UNIT=11, FILE=CFILE_IN, FORM='UNFORMATTED', STATUS='OLD', + ACCESS='DIRECT', RECL=NRECL, IOSTAT=IOS) IF (IOS.NE.0) THEN write(6,*) 'Error: can''t open ',TRIM(CFILE_IN) write(6,*) 'ios = ',ios write(6,*) 'nrecl = ',nrecl CALL EXIT(3) ENDIF C------------------------------------------------------------------------ C Open OUTPUT File C------------------------------------------------------------------------ C OUTPUT FILE TIDAL(IDM,JDM,N2MODES) C C First Open the .a and .b files C OPEN(UNIT=21, FILE=CFILE_OUT, FORM='UNFORMATTED', STATUS='NEW', + ACCESS='DIRECT', RECL=NRECL, IOSTAT=IOS) IF (IOS.NE.0) THEN write(6,*) 'Error: can''t open ',TRIM(CFILE_OUT) write(6,*) 'ios = ',ios write(6,*) 'nrecl = ',nrecl CALL EXIT(3) ENDIF WRITE(6,*)'Output TidePAout.a File Opened,IOS =',IOS C read(5,*)IGG C CFILEB = CFILE_OUT(1:LEN_TRIM(CFILE_OUT)-1) // 'b' OPEN(UNIT=22,FILE=CFILEB,FORM='FORMATTED', & STATUS='NEW',ACTION='WRITE',IOSTAT=IOS) WRITE(6,*)'TidaPAout.b File Opened,IOS =',IOS IF (IOS.NE.0) THEN write(6,*) 'Error: can''t open ',TRIM(CFILEB) write(6,*) 'ios = ',ios CALL EXIT(3) ENDIF WRITE(6,*)'TidePAout.b File Opened' C read(5,*)IGG C C========================================================================= C Loop Through Pairs of Tidal Mode Real and Imag Fields until end C IREC=0 100 CONTINUE IREC=IREC+2 READ(11,REC=IREC-1,IOSTAT=IOS,ERR=200) TideR write(6,*)'Array ',IREC,' read, NRECL=',NRECL #ifdef ENDIAN_IO CALL ENDIAN_SWAP(TideR,IDM*JDM) #endif IF (IOS.NE.0) THEN WRITE(6,*) 'can''t read TideR on ',TRIM(CFILE_IN) c CALL EXIT(4) GO TO 200 ENDIF c READ(11,REC=IREC,IOSTAT=IOS,ERR=200) TideI write(6,*)'Array TideI read, NRECL=',NRECL #ifdef ENDIAN_IO CALL ENDIAN_SWAP(TideI,IDM*JDM) #endif IF (IOS.NE.0) THEN WRITE(6,*) 'can''t read TideI on ',TRIM(CFILE_IN) c CALL EXIT(4) GO TO 200 ENDIF C C Now convert Re Im to Amp Phase C Tmin= 1.e10 Tmax=-1.e10 Pmin= 1.e10 Pmax=-1.e10 DO J=1,JDM DO I=1,IDM IF(TideR(I,J).eq.SPVAL)THEN TideA(I,J)=SPVAL TideP(I,J)=SPVAL ELSE TideA(I,J)=sqrt(TideR(I,J)**2+TideI(I,J)**2) Tmax=max(Tmax,TideA(I,J)) Tmin=min(Tmin,TideA(I,J)) Phase=180.*ATAN2(-TideI(I,J),TideR(I,J))/PI !OSU convention IF(.not.L180 .and. Phase.lt.0.)Phase=Phase+360. TideP(I,J)=Phase Pmax=max(Pmax,TideP(I,J)) Pmin=min(Pmin,TideP(I,J)) ENDIF ENDDO ENDDO C C Now Write Out Amp and Phase Fields C WRITE(21,REC=IREC-1,IOSTAT=IOS)TideA,(PAD(I),I=1,NPAD) WRITE(21,REC=IREC ,IOSTAT=IOS)TideP,(PAD(I),I=1,NPAD) WRITE(22,62)IREC/2,Tmin,Tmax WRITE( 6,62)IREC/2,Tmin,Tmax WRITE(22,63)IREC/2,Pmin,Pmax WRITE( 6,63)IREC/2,Pmin,Pmax C----------------------------------------------------------------- C Print our Values at Test Point c IF(IP_test.gt.0)then WRITE(6,64)IP_test,JP_test,IREC/2,TideR(IP_test,JP_test), + TideI(IP_test,JP_test),TideA(IP_test,JP_test), + TideP(IP_test,JP_test) end if cC========================================================================= C Print Out a trace of the Field C * DO J=2,JDM,8 * WRITE(6,'(20F6.2)')(TideA(I,J),I=2,IDM,9) * END DO * DO J=2,JDM,8 * WRITE(6,'(20F6.1)')(TideP(I,J),I=2,IDM,9) * END DO * READ(5,*)IGG C Loop- to next field c GO TO 100 62 FORMAT('TIDE',I2.2,'Am min, max=',2g15.7) 63 FORMAT('TIDE',I2.2,'Ph min, max=',2g15.7) 64 FORMAT(' At point: (',i4,',',i4,'), Re Im Amp Ph for mode ',I2, + ' = ',4g15.6) C#################################################################### c c 200 CONTINUE WRITE(6,*)IREC/2-1,' Tidal Re, Im fields converted to Amp, Phase' close(21) close(22) CALL EXIT(0) END