Skip to content

Instantly share code, notes, and snippets.

@leverich

leverich/tide_fac.f

Created November 27, 2015 19:02
Show Gist options
  • Save leverich/e0834df944d457962a4e to your computer and use it in GitHub Desktop.
Save leverich/e0834df944d457962a4e to your computer and use it in GitHub Desktop.
Download ZIP
tide_fac.f modified to compute node factors and equilibrium arguments for all 37 NOS/NOAA tidal constituents
Raw
tide_fac.f
C PROGRAM TO COMPUTE NODAL FACTORS AND EQUILIBRIUM ARGUEMENTS
C
C
PARAMETER(NCNST=37)
CHARACTER CNAME(NCNST)*8
COMMON /CNSNAM/ CNAME
REAL NODFAC,MONTH
DIMENSION NCON(NCNST)
COMMON /CNST/ NODFAC(NCNST),GRTERM(NCNST),SPEED(NCNST),P(NCNST)
OPEN(UNIT=11,FILE='tide_fac.out',STATUS='UNKNOWN')
WRITE(*,*) 'ENTER LENGTH OF RUN TIME (DAYS)'
READ(*,*) XDAYS
RHRS=XDAYS*24.
WRITE(*,*)' ENTER START TIME - BHR,IDAY,IMO,IYR (IYR e.g. 1992)'
READ(*,*) BHR,IDAY,IMO,IYR
YR=IYR
MONTH=IMO
DAY=IDAY
HRM=BHR+RHRS/2.
WRITE(11,10) BHR,IDAY,IMO,IYR
WRITE(*,10) BHR,IDAY,IMO,IYR
10 FORMAT(' TIDAL FACTORS STARTING: ',
& ' HR-',F5.2,', DAY-',I3,', MONTH-',I3,' YEAR-',I5,/)
WRITE(*,11) XDAYS
WRITE(11,11) XDAYS
11 FORMAT(' FOR A RUN LASTING ',F8.2,' DAYS',//)
C-- DETERMINE THE JULIAN TIME AT BEGINNING AND MIDDLE OF RECORD
DAYJ=DAYJUL(YR,MONTH,DAY)
C-- DETERMINE NODE FACTORS AT MIDDLE OF RECORD
CALL NFACS(YR,DAYJ,HRM)
C-- DETERMINE GREENWICH EQUIL. TERMS AT BEGINNING OF RECORD
CALL GTERMS(YR,DAYJ,BHR,DAYJ,HRM)
NUMCON=37
NCON(1)=1
NCON(2)=2
NCON(3)=3
NCON(4)=4
NCON(5)=5
NCON(6)=6
NCON(7)=7
NCON(8)=8
NCON(9)=9
NCON(10)=10
NCON(11)=11
NCON(12)=12
NCON(13)=13
NCON(14)=14
NCON(15)=15
NCON(16)=16
NCON(17)=17
NCON(18)=18
NCON(19)=19
NCON(20)=20
NCON(21)=21
NCON(22)=22
NCON(23)=23
NCON(24)=24
NCON(25)=25
NCON(26)=26
NCON(27)=27
NCON(28)=28
NCON(29)=29
NCON(30)=30
NCON(31)=31
NCON(32)=32
NCON(33)=33
NCON(34)=34
NCON(35)=35
NCON(36)=36
NCON(37)=37
WRITE(11,*) 'CONST NODE EQ ARG (ref GM)'
WRITE(11,1300)
1300 FORMAT(' NAME FACTOR (DEG) ',//)
DO 20 NC=1,NUMCON
IC=NCON(NC)
C EQUILIBRIUM ARGUEMENT IS REFERENCED TO THE GRENWICH MERIDIAN
WRITE(11,2001) CNAME(IC),NODFAC(IC),GRTERM(IC)
2001 FORMAT(1X,A4,2x,F7.5,4x,F7.2,2x,F7.4)
20 CONTINUE
STOP
END
SUBROUTINE NFACS(YR,DAYJ,HR)
C-- CALCULATES NODE FACTORS FOR CONSTITUENT TIDAL SIGNAL
C-- THE EQUATIONS USED IN THIS ROUTINE COME FROM:
C "MANUAL OF HARMONIC ANALYSIS AND PREDICTION OF TIDES"
C BY PAUL SCHUREMAN, SPECIAL PUBLICATION #98, US COAST
C AND GEODETIC SURVEY, DEPARTMENT OF COMMERCE (1958).
C-- IF DAYM AND HRM CORRESPOND TO MIDYEAR, THEN THIS ROUTINE
C-- RETURNS THE SAME VALUES AS FOUND IN TABLE 14 OF SCHUREMAN.
C---------------------------------------------------------------------
CHARACTER*8 CST(37)
REAL I,N,NU
COMMON/ORBITF/DS,DP,DH,DP1,DN,DI,DNU,DXI,DNUP,DNUP2,DPC
COMMON/ CNST /FNDCST(37),EQCST(37),ACST(37),PCST(37)
COMMON/CNSNAM/CST
C-- CONSTITUENT NAMES:
DATA CST /'M2 ','S2 ','N2 ','K1 ',
* 'M4 ','O1 ','M6 ','MK3 ',
* 'S4 ','MN4 ','NU2 ','S6 ',
* 'MU2 ','2N2 ','OO1 ','LAMBDA2 ',
* 'S1 ','M1 ','J1 ','MM ',
* 'SSA ','SA ','MSF ','MF ',
* 'RHO1 ','Q1 ','T2 ','R2 ',
* '2Q1 ','P1 ','2SM2 ','M3 ',
* 'L2 ','2MK3 ','K2 ','M8 ',
* 'MS4 '/
C-- ORBITAL SPEEDS (DEGREES/HOUR):
DATA ACST/28.9841042,30.0,28.4397295,15.0410686,57.9682084,
*13.9430356,86.9523127,44.0251729,60.0,57.4238337,28.5125831,90.0,
*27.9682084,27.8953548,16.1391017,29.4556253,15.0,14.4966939,
*15.5854433,0.5443747,0.0821373,0.0410686,1.0158958,1.0980331,
*13.4715145,13.3986609,29.9589333,30.0410667,12.8542862,14.9589314,
*31.0158958,43.4761563,29.5284789,42.9271398,30.0821373,
*115.9364169,58.9841042/
C-- NUMBER OF TIDE CYCLES PER DAY PER CONSTITUENT:
DATA PCST/2.,2.,2.,1.,4.,1.,6.,3.,4.,4.,2.,6.,2.,2.,1.,2.,1.,1.,
$1.,0.,0.,0.,0.,0.,1.,1.,2.,2.,1.,1.,2.,3.,2.,3.,2.,8.,4./
PI180=3.14159265/180.
CALL ORBIT(YR,DAYJ,HR)
N=DN*PI180
I=DI*PI180
NU=DNU*PI180
XI=DXI*PI180
P=DP*PI180
PC=DPC*PI180
SINI=SIN(I)
SINI2=SIN(I/2.)
SIN2I=SIN(2.*I)
COSI2=COS(I/2.)
TANI2=TAN(I/2.)
C-- EQUATION 197, SCHUREMAN
QAINV=SQRT(2.310+1.435*COS(2.*PC))
C-- EQUATION 213, SCHUREMAN
RAINV=SQRT(1.-12.*TANI2**2*COS(2.*PC)+36.*TANI2**4)
C-- VARIABLE NAMES REFER TO EQUATION NUMBERS IN SCHUREMAN
EQ73=(2./3.-SINI**2)/.5021
EQ74=SINI**2/.1578
EQ75=SINI*COSI2**2/.37988
EQ76=SIN(2*I)/.7214
EQ77=SINI*SINI2**2/.0164
EQ78=(COSI2**4)/.91544
EQ149=COSI2**6/.8758
EQ207=EQ75*QAINV
EQ215=EQ78*RAINV
EQ227=SQRT(.8965*SIN2I**2+.6001*SIN2I*COS(NU)+.1006)
EQ235=.001+SQRT(19.0444*SINI**4+2.7702*SINI**2*COS(2.*NU)+.0981)
C-- NODE FACTORS FOR 37 CONSTITUENTS:
FNDCST(1)=EQ78
FNDCST(2)=1.0
FNDCST(3)=EQ78
FNDCST(4)=EQ227
FNDCST(5)=FNDCST(1)**2
FNDCST(6)=EQ75
FNDCST(7)=FNDCST(1)**3
FNDCST(8)=FNDCST(1)*FNDCST(4)
FNDCST(9)=1.0
FNDCST(10)=FNDCST(1)**2
FNDCST(11)=EQ78
FNDCST(12)=1.0
FNDCST(13)=EQ78
FNDCST(14)=EQ78
FNDCST(15)=EQ77
FNDCST(16)=EQ78
FNDCST(17)=1.0
C** EQUATION 207 NOT PRODUCING CORRECT ANSWER FOR M1
C**SET NODE FACTOR FOR M1 = 0 UNTIL CAN FURTHER RESEARCH
FNDCST(18)=0.
C FNDCST(18)=EQ207
FNDCST(19)=EQ76
FNDCST(20)=EQ73
FNDCST(21)=1.0
FNDCST(22)=1.0
FNDCST(23)=EQ78
FNDCST(24)=EQ74
FNDCST(25)=EQ75
FNDCST(26)=EQ75
FNDCST(27)=1.0
FNDCST(28)=1.0
FNDCST(29)=EQ75
FNDCST(30)=1.0
FNDCST(31)=EQ78
FNDCST(32)=EQ149
C** EQUATION 215 NOT PRODUCING CORRECT ANSWER FOR L2
C** SET NODE FACTOR FOR L2 = 0 UNTIL CAN FURTHER RESEARCH
FNDCST(33)=0.
C FNDCST(33)=EQ215
FNDCST(34)=FNDCST(1)**2*FNDCST(4)
FNDCST(35)=EQ235
FNDCST(36)=FNDCST(1)**4
FNDCST(37)=EQ78
END
SUBROUTINE GTERMS(YR,DAYJ,HR,DAYM,HRM)
C-- CALCULATES EQUILIBRIUM ARGUMENTS V0+U FOR CONSTITUENT TIDE
C-- THE EQUATIONS USED IN THIS ROUTINE COME FROM:
C "MANUAL OF HARMONIC ANALYSIS AND PREDICTION OF TIDES"
C BY PAUL SCHUREMAN, SPECIAL PUBLICATION #98, US COAST
C AND GEODETIC SURVEY, DEPARTMENT OF COMMERCE (1958).
C-- IF DAYM AND HRM CORRESPOND TO MIDYEAR, THEN THIS ROUTINE
C-- RETURNS THE SAME VALUES AS FOUND IN TABLE 15 OF SCHUREMAN.
C---------------------------------------------------------------------
REAL NU,NUP,NUP2,I
COMMON /ORBITF/DS,DP,DH,DP1,DN,DI,DNU,DXI,DNUP,DNUP2,DPC
COMMON /CNST/ FNDCST(37),EQCST(37),ACST(37),PCST(37)
PI180=3.14159265/180.
C* OBTAINING ORBITAL VALUES AT BEGINNING OF SERIES FOR V0
CALL ORBIT(YR,DAYJ,HR)
S=DS
P=DP
H=DH
P1=DP1
T=ANGLE(180.+HR*(360./24.))
C** OBTAINING ORBITAL VALUES AT MIDDLE OF SERIES FOR U
CALL ORBIT(YR,DAYM,HRM)
NU=DNU
XI=DXI
NUP=DNUP
NUP2=DNUP2
C* SUMMING TERMS TO OBTAIN EQUILIBRIUM ARGUMENTS
EQCST(1)=2.*(T-S+H)+2.*(XI-NU)
EQCST(2)=2.*T
EQCST(3)=2.*(T+H)-3.*S+P+2.*(XI-NU)
EQCST(4)=T+H-90.-NUP
EQCST(5)=4.*(T-S+H)+4.*(XI-NU)
EQCST(6)=T-2.*S+H+90.+2.*XI-NU
EQCST(7)=6.*(T-S+H)+6.*(XI-NU)
EQCST(8)=3.*(T+H)-2.*S-90.+2.*(XI-NU)-NUP
EQCST(9)=4.*T
EQCST(10)=4.*(T+H)-5.*S+P+4.*(XI-NU)
EQCST(11)=2.*T-3.*S+4.*H-P+2.*(XI-NU)
EQCST(12)=6.*T
EQCST(13)=2.*(T+2.*(H-S))+2.*(XI-NU)
EQCST(14)=2.*(T-2.*S+H+P)+2.*(XI-NU)
EQCST(15)=T+2.*S+H-90.-2.*XI-NU
EQCST(16)=2.*T-S+P+180.+2.*(XI-NU)
EQCST(17)=T
I=DI*PI180
PC=DPC*PI180
TOP=(5.*COS(I)-1.)*SIN(PC)
BOTTOM=(7.*COS(I)+1.)*COS(PC)
Q=ARCTAN(TOP,BOTTOM,1)
EQCST(18)=T-S+H-90.+XI-NU+Q
EQCST(19)=T+S+H-P-90.-NU
EQCST(20)=S-P
EQCST(21)=2.*H
EQCST(22)=H
EQCST(23)=2.*(S-H)
EQCST(24)=2.*S-2.*XI
EQCST(25)=T+3.*(H-S)-P+90.+2.*XI-NU
EQCST(26)=T-3.*S+H+P+90.+2.*XI-NU
EQCST(27)=2.*T-H+P1
EQCST(28)=2.*T+H-P1+180.
EQCST(29)=T-4.*S+H+2.*P+90.+2.*XI-NU
EQCST(30)=T-H+90.
EQCST(31)=2.*(T+S-H)+2.*(NU-XI)
EQCST(32)=3.*(T-S+H)+3.*(XI-NU)
R=SIN(2.*PC)/((1./6.)*(1./TAN(.5*I))**2-COS(2.*PC))
R=ATAN(R)/PI180
EQCST(33)=2.*(T+H)-S-P+180.+2.*(XI-NU)-R
EQCST(34)=3.*(T+H)-4.*S+90.+4.*(XI-NU)+NUP
EQCST(35)=2.*(T+H)-2.*NUP2
EQCST(36)=8.*(T-S+H)+8.*(XI-NU)
EQCST(37)=2.*(2.*T-S+H)+2.*(XI-NU)
DO 1 IH=1,37
1 EQCST(IH)=ANGLE(EQCST(IH))
END
SUBROUTINE ORBIT(YR,DAYJ,HR)
C-- DETERMINATION OF PRIMARY AND SECONDARY ORBITAL FUNCTIONS
C-- THE EQUATIONS PROGRAMMED HERE ARE NOT REPRESENTED BY EQUATIONS IN
C SCHUREMAN. THE CODING IN THIS ROUTINE DERIVES FROM A PROGRAM BY
C THE NATIONAL OCEANIC AND ATMOSPHERIC ADMINISTRATION (NOAA).
C HOWEVER, TABULAR VALUES OF THE ORBITAL FUNCTIONS CAN BE FOUND IN
C TABLE 1 OF SCHUREMAN.
C---------------------------------------------------------------------
REAL I,N,NU,NUP,NUP2
COMMON /ORBITF/DS,DP,DH,DP1,DN,DI,DNU,DXI,DNUP,DNUP2,DPC
PI180=3.14159265/180.
X=AINT((YR-1901.)/4.)
DYR=YR-1900.
DDAY=DAYJ+X-1.
C-- DN IS THE MOON'S NODE (CAPITAL N, TABLE 1, SCHUREMAN)
DN=259.1560564-19.328185764*DYR-.0529539336*DDAY-.0022064139*HR
DN=ANGLE(DN)
N=DN*PI180
C-- DP IS THE LUNAR PERIGEE (SMALL P, TABLE 1)
DP=334.3837214+40.66246584*DYR+.111404016*DDAY+.004641834*HR
DP=ANGLE(DP)
P=DP*PI180
I=ACOS(.9136949-.0356926*COS(N))
DI=ANGLE(I/PI180)
NU=ASIN(.0897056*SIN(N)/SIN(I))
DNU=NU/PI180
XI=N-2.*ATAN(.64412*TAN(N/2.))-NU
DXI=XI/PI180
DPC=ANGLE(DP-DXI)
C-- DH IS THE MEAN LONGITUDE OF THE SUN (SMALL H, TABLE 1)
DH=280.1895014-.238724988*DYR+.9856473288*DDAY+.0410686387*HR
DH=ANGLE(DH)
C-- DP1 IS THE SOLAR PERIGEE (SMALL P1, TABLE 1)
DP1=281.2208569+.01717836*DYR+.000047064*DDAY+.000001961*HR
DP1=ANGLE(DP1)
C-- DS IS THE MEAN LONGITUDE OF THE MOON (SMALL S, TABLE 1)
DS=277.0256206+129.38482032*DYR+13.176396768*DDAY+.549016532*HR
DS=ANGLE(DS)
NUP=ATAN(SIN(NU)/(COS(NU)+.334766/SIN(2.*I)))
DNUP=NUP/PI180
NUP2=ATAN(SIN(2.*NU)/(COS(2.*NU)+.0726184/SIN(I)**2))/2.
DNUP2=NUP2/PI180
END
FUNCTION ANGLE(ARG)
C
C*** THIS ROUTINE PLACES AN ANGLE IN 0-360 (+) FORMAT
C
M=-IFIX(ARG/360.)
ANGLE=ARG+FLOAT(M)*360.
IF(ANGLE .LT. 0.) ANGLE=ANGLE+360.
END
FUNCTION ARCTAN(TOP,BOTTOM,KEY)
C** DETERMINE ARCTANGENT AND PLACE IN CORRECT QUADRANT
C IF KEY EQ 0 NO QUADRANT SELECTION MADE
C IF KEY .NE. 0 PROPER QUADRANT IS SELECTED
IF(BOTTOM .NE. 0.0) GO TO 4
IF(TOP) 2,9,3
2 ARCTAN=270.
RETURN
3 ARCTAN=90.
RETURN
4 ARCTAN=ATAN(TOP/BOTTOM)*57.2957795
IF(KEY.EQ.0) RETURN
IF(TOP) 5,5,7
5 IF(BOTTOM) 6,9,8
6 ARCTAN=ARCTAN+180.
RETURN
7 IF(BOTTOM) 6,3,10
8 ARCTAN=ARCTAN+360.
RETURN
9 ARCTAN=0.
10 RETURN
END
FUNCTION DAYJUL(YR,XMONTH,DAY)
C
C*** THIS ROUTINE COMPUTES THE JULIAN DAY (AS A REAL VARIABLE)
C
DIMENSION DAYT(12),DAYS(12)
DATA DAYT/0.,31.,59.,90.,120.,151.,181.,212.,243.,273.,304.,334./
DATA DAYS(1),DAYS(2) /0.,31./
DINC=0.
YRLP=MOD((YR-1900.),4.)
IF(YRLP .EQ. 0.) DINC=1.
DO 1 I=3,12
1 DAYS(I)=DAYT(I)+DINC
DAYJUL=DAYS(IFIX(XMONTH))+DAY
END
@deguc001
Copy link

deguc001 commented Sep 5, 2024

I tried executing the code in Silverfrost FTN95 and it did ask for parameters like run time, hour, day, month, and year, but showed no results in the end. Does it need to be a particular application in order to run the code for the computation? Thanks.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment