Skip to content

Commit 9ed3c88

Browse files
AliS-NoaaAliS-Noaa
authored
Merge branch 'NOAA-EMC:develop' into outp-access
2 parents bb97240 + d6160ab commit 9ed3c88

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

75 files changed

+9698
-9712
lines changed

‎model/src/constants.F90‎

Lines changed: 126 additions & 128 deletions
Original file line numberDiff line numberDiff line change
@@ -34,6 +34,7 @@ MODULE CONSTANTS
3434
!/ 20-Jan-2017 : Add parameters for ESMF ( version 6.02 )
3535
!/ 01-Mar-2018 : Add UNDEF parameter ( version 6.02 )
3636
!/ 05-Jun-2018 : Add PDLIB parameters ( version 6.04 )
37+
!/ 04-Jul-2025 : Remove labelled statements ( version X.XX )
3738
!/
3839
!/ Copyright 2009-2012 National Weather Service (NWS),
3940
!/ National Oceanic and Atmospheric Administration. All rights
@@ -275,134 +276,131 @@ SUBROUTINE KZEONE(X, Y, RE0, IM0, RE1, IM1)
275276
! ABSCISSAS AND THE WEIGHT FACTORS USED IN THE GAUSS-
276277
! HERMITE QUADRATURE.
277278
R2 = X*X + Y*Y
278-
IF (R2.GE.1.96D2) GO TO 50
279-
IF (R2.GE.1.849D1) GO TO 30
280-
! THIS SECTION CALCULATES THE FUNCTIONS USING THE SERIES
281-
! EXPANSIONS
282-
X2 = X/2.0D0
283-
Y2 = Y/2.0D0
284-
P1 = X2*X2
285-
P2 = Y2*Y2
286-
T1 = -(DLOG(P1+P2)/2.0D0+0.5772156649015329D0)
287-
! THE CONSTANT IN THE PRECEDING STATEMENT IS EULER*S
288-
! CONSTANT
289-
T2 = -DATAN2(Y,X)
290-
X2 = P1 - P2
291-
Y2 = X*Y2
292-
RTERM = 1.0D0
293-
ITERM = 0.0D0
294-
RE0 = T1
295-
IM0 = T2
296-
T1 = T1 + 0.5D0
297-
RE1 = T1
298-
IM1 = T2
299-
P2 = DSQRT(R2)
300-
L = 2.106D0*P2 + 4.4D0
301-
IF (P2.LT.8.0D-1) L = 2.129D0*P2 + 4.0D0
302-
DO N=1,INT(L)
303-
P1 = N
304-
P2 = N*N
305-
R1 = RTERM
306-
RTERM = (R1*X2-ITERM*Y2)/P2
307-
ITERM = (R1*Y2+ITERM*X2)/P2
308-
T1 = T1 + 0.5D0/P1
309-
RE0 = RE0 + T1*RTERM - T2*ITERM
310-
IM0 = IM0 + T1*ITERM + T2*RTERM
311-
P1 = P1 + 1.0D0
312-
T1 = T1 + 0.5D0/P1
313-
RE1 = RE1 + (T1*RTERM-T2*ITERM)/P1
314-
IM1 = IM1 + (T1*ITERM+T2*RTERM)/P1
315-
END DO
316-
R1 = X/R2 - 0.5D0*(X*RE1-Y*IM1)
317-
R2 = -Y/R2 - 0.5D0*(X*IM1+Y*RE1)
318-
P1 = DEXP(X)
319-
RE0 = P1*RE0
320-
IM0 = P1*IM0
321-
RE1 = P1*R1
322-
IM1 = P1*R2
323-
RETURN
324-
! THIS SECTION CALCULATES THE FUNCTIONS USING THE INTEGRAL
325-
! REPRESENTATION, EQN 3, EVALUATED WITH 15 POINT GAUSS-
326-
! HERMITE QUADRATURE
327-
30 X2 = 2.0D0*X
328-
Y2 = 2.0D0*Y
329-
R1 = Y2*Y2
330-
P1 = DSQRT(X2*X2+R1)
331-
P2 = DSQRT(P1+X2)
332-
T1 = EXSQ(1)/(2.0D0*P1)
333-
RE0 = T1*P2
334-
IM0 = T1/P2
335-
RE1 = 0.0D0
336-
IM1 = 0.0D0
337-
DO N=2,8
338-
T2 = X2 + TSQ(N)
339-
P1 = DSQRT(T2*T2+R1)
340-
P2 = DSQRT(P1+T2)
341-
T1 = EXSQ(N)/P1
342-
RE0 = RE0 + T1*P2
343-
IM0 = IM0 + T1/P2
344-
T1 = EXSQ(N)*TSQ(N)
345-
RE1 = RE1 + T1*P2
346-
IM1 = IM1 + T1/P2
347-
END DO
348-
T2 = -Y2*IM0
349-
RE1 = RE1/R2
350-
R2 = Y2*IM1/R2
351-
RTERM = 1.41421356237309D0*DCOS(Y)
352-
ITERM = -1.41421356237309D0*DSIN(Y)
353-
! THE CONSTANT IN THE PREVIOUS STATEMENTS IS,OF COURSE,
354-
! SQRT(2.0).
355-
IM0 = RE0*ITERM + T2*RTERM
356-
RE0 = RE0*RTERM - T2*ITERM
357-
T1 = RE1*RTERM - R2*ITERM
358-
T2 = RE1*ITERM + R2*RTERM
359-
RE1 = T1*X + T2*Y
360-
IM1 = -T1*Y + T2*X
361-
RETURN
362-
! THIS SECTION CALCULATES THE FUNCTIONS USING THE
363-
! ASYMPTOTIC EXPANSIONS
364-
50 RTERM = 1.0D0
365-
ITERM = 0.0D0
366-
RE0 = 1.0D0
367-
IM0 = 0.0D0
368-
RE1 = 1.0D0
369-
IM1 = 0.0D0
370-
P1 = 8.0D0*R2
371-
P2 = DSQRT(R2)
372-
L = 3.91D0+8.12D1/P2
373-
R1 = 1.0D0
374-
R2 = 1.0D0
375-
M = -8
376-
K = 3
377-
DO N=1,INT(L)
378-
M = M + 8
379-
K = K - M
380-
R1 = FLOAT(K-4)*R1
381-
R2 = FLOAT(K)*R2
382-
T1 = FLOAT(N)*P1
383-
T2 = RTERM
384-
RTERM = (T2*X+ITERM*Y)/T1
385-
ITERM = (-T2*Y+ITERM*X)/T1
386-
RE0 = RE0 + R1*RTERM
387-
IM0 = IM0 + R1*ITERM
388-
RE1 = RE1 + R2*RTERM
389-
IM1 = IM1 + R2*ITERM
390-
END DO
391-
T1 = DSQRT(P2+X)
392-
T2 = -Y/T1
393-
P1 = 8.86226925452758D-1/P2
394-
! THIS CONSTANT IS SQRT(PI)/2.0, WITH PI=3.14159...
395-
RTERM = P1*DCOS(Y)
396-
ITERM = -P1*DSIN(Y)
397-
R1 = RE0*RTERM - IM0*ITERM
398-
R2 = RE0*ITERM + IM0*RTERM
399-
RE0 = T1*R1 - T2*R2
400-
IM0 = T1*R2 + T2*R1
401-
R1 = RE1*RTERM - IM1*ITERM
402-
R2 = RE1*ITERM + IM1*RTERM
403-
RE1 = T1*R1 - T2*R2
404-
IM1 = T1*R2 + T2*R1
405-
RETURN
279+
IF (R2.GE.1.96D2) THEN
280+
! CALCULATE THE FUNCTIONS USING THE ASYMPTOTIC EXPANSIONS
281+
RTERM = 1.0D0
282+
ITERM = 0.0D0
283+
RE0 = 1.0D0
284+
IM0 = 0.0D0
285+
RE1 = 1.0D0
286+
IM1 = 0.0D0
287+
P1 = 8.0D0*R2
288+
P2 = DSQRT(R2)
289+
L = 3.91D0+8.12D1/P2
290+
R1 = 1.0D0
291+
R2 = 1.0D0
292+
M = -8
293+
K = 3
294+
DO N=1,INT(L)
295+
M = M + 8
296+
K = K - M
297+
R1 = FLOAT(K-4)*R1
298+
R2 = FLOAT(K)*R2
299+
T1 = FLOAT(N)*P1
300+
T2 = RTERM
301+
RTERM = (T2*X+ITERM*Y)/T1
302+
ITERM = (-T2*Y+ITERM*X)/T1
303+
RE0 = RE0 + R1*RTERM
304+
IM0 = IM0 + R1*ITERM
305+
RE1 = RE1 + R2*RTERM
306+
IM1 = IM1 + R2*ITERM
307+
END DO
308+
T1 = DSQRT(P2+X)
309+
T2 = -Y/T1
310+
P1 = 8.86226925452758D-1/P2
311+
! THIS CONSTANT IS SQRT(PI)/2.0, WITH PI=3.14159...
312+
RTERM = P1*DCOS(Y)
313+
ITERM = -P1*DSIN(Y)
314+
R1 = RE0*RTERM - IM0*ITERM
315+
R2 = RE0*ITERM + IM0*RTERM
316+
RE0 = T1*R1 - T2*R2
317+
IM0 = T1*R2 + T2*R1
318+
R1 = RE1*RTERM - IM1*ITERM
319+
R2 = RE1*ITERM + IM1*RTERM
320+
RE1 = T1*R1 - T2*R2
321+
IM1 = T1*R2 + T2*R1
322+
ELSE IF (R2.GE.1.849D1) THEN
323+
! CALCULATE THE FUNCTIONS USING THE INTEGRAL
324+
! REPRESENTATION, EQN 3, EVALUATED WITH 15 POINT GAUSS-
325+
! HERMITE QUADRATURE
326+
X2 = 2.0D0*X
327+
Y2 = 2.0D0*Y
328+
R1 = Y2*Y2
329+
P1 = DSQRT(X2*X2+R1)
330+
P2 = DSQRT(P1+X2)
331+
T1 = EXSQ(1)/(2.0D0*P1)
332+
RE0 = T1*P2
333+
IM0 = T1/P2
334+
RE1 = 0.0D0
335+
IM1 = 0.0D0
336+
DO N=2,8
337+
T2 = X2 + TSQ(N)
338+
P1 = DSQRT(T2*T2+R1)
339+
P2 = DSQRT(P1+T2)
340+
T1 = EXSQ(N)/P1
341+
RE0 = RE0 + T1*P2
342+
IM0 = IM0 + T1/P2
343+
T1 = EXSQ(N)*TSQ(N)
344+
RE1 = RE1 + T1*P2
345+
IM1 = IM1 + T1/P2
346+
END DO
347+
T2 = -Y2*IM0
348+
RE1 = RE1/R2
349+
R2 = Y2*IM1/R2
350+
RTERM = 1.41421356237309D0*DCOS(Y)
351+
ITERM = -1.41421356237309D0*DSIN(Y)
352+
! THE CONSTANT IN THE PREVIOUS STATEMENTS IS,OF COURSE,
353+
! SQRT(2.0).
354+
IM0 = RE0*ITERM + T2*RTERM
355+
RE0 = RE0*RTERM - T2*ITERM
356+
T1 = RE1*RTERM - R2*ITERM
357+
T2 = RE1*ITERM + R2*RTERM
358+
RE1 = T1*X + T2*Y
359+
IM1 = -T1*Y + T2*X
360+
ELSE
361+
! CALCULATE THE FUNCTIONS USING THE SERIES EXPANSIONS
362+
X2 = X/2.0D0
363+
Y2 = Y/2.0D0
364+
P1 = X2*X2
365+
P2 = Y2*Y2
366+
T1 = -(DLOG(P1+P2)/2.0D0+0.5772156649015329D0)
367+
! THE CONSTANT IN THE PRECEDING STATEMENT IS EULER*S
368+
! CONSTANT
369+
T2 = -DATAN2(Y,X)
370+
X2 = P1 - P2
371+
Y2 = X*Y2
372+
RTERM = 1.0D0
373+
ITERM = 0.0D0
374+
RE0 = T1
375+
IM0 = T2
376+
T1 = T1 + 0.5D0
377+
RE1 = T1
378+
IM1 = T2
379+
P2 = DSQRT(R2)
380+
L = 2.106D0*P2 + 4.4D0
381+
IF (P2.LT.8.0D-1) L = 2.129D0*P2 + 4.0D0
382+
DO N=1,INT(L)
383+
P1 = N
384+
P2 = N*N
385+
R1 = RTERM
386+
RTERM = (R1*X2-ITERM*Y2)/P2
387+
ITERM = (R1*Y2+ITERM*X2)/P2
388+
T1 = T1 + 0.5D0/P1
389+
RE0 = RE0 + T1*RTERM - T2*ITERM
390+
IM0 = IM0 + T1*ITERM + T2*RTERM
391+
P1 = P1 + 1.0D0
392+
T1 = T1 + 0.5D0/P1
393+
RE1 = RE1 + (T1*RTERM-T2*ITERM)/P1
394+
IM1 = IM1 + (T1*ITERM+T2*RTERM)/P1
395+
END DO
396+
R1 = X/R2 - 0.5D0*(X*RE1-Y*IM1)
397+
R2 = -Y/R2 - 0.5D0*(X*IM1+Y*RE1)
398+
P1 = DEXP(X)
399+
RE0 = P1*RE0
400+
IM0 = P1*IM0
401+
RE1 = P1*R1
402+
IM1 = P1*R2
403+
END IF
406404
END SUBROUTINE KZEONE
407405

408406
! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

0 commit comments

Comments
 (0)