Routine weight_gen uses splib routine splat. As a test, if

the lapack_gen.F and splat.F routines from splib are made local files, the dependency on splib can be removed from the makefile. Fixes ufs-community#886.
GeorgeGayno-NOAA · Jan 18, 2024 · 631a278 · 631a278
1 parent 35670c8
commit 631a278
Show file tree

Hide file tree

Showing 3 changed files with 310 additions and 1 deletion.
diff --git a/sorc/weight_gen.fd/CMakeLists.txt b/sorc/weight_gen.fd/CMakeLists.txt
@@ -1,5 +1,7 @@
 list(APPEND fortran_src
   scrip.F90
+  lapack_gen.F
+  splat.F
 )
 
 if(CMAKE_Fortran_COMPILER_ID MATCHES "^(Intel)$")
@@ -12,7 +14,6 @@ set(exe_name weight_gen)
 add_executable(${exe_name} ${fortran_src})
 target_link_libraries(
   ${exe_name}
-  sp::sp_d
   NetCDF::NetCDF_Fortran)
 
 install(TARGETS ${exe_name})

diff --git a/sorc/weight_gen.fd/lapack_gen.F b/sorc/weight_gen.fd/lapack_gen.F
@@ -0,0 +1,116 @@
+C> @file
+C> @brief Two Numerical Recipes routines for matrix inversion From Numerical Recipes.
+C>      
+C> ### Program History Log
+C> Date | Programmer | Comments
+C> -----|------------|---------
+C> 2012-11-05 | E.Mirvis | separated this generic LU from the splat.F
+
+C> Solves a system of linear equations, follows call to ludcmp().
+C>
+C> @param A
+C> @param N
+C> @param NP
+C> @param INDX
+C> @param B
+      SUBROUTINE LUBKSB(A,N,NP,INDX,B)
+      REAL A(NP,NP),B(N)
+      INTEGER INDX(N)
+      II=0  
+      DO 12 I=1,N
+        LL=INDX(I)
+        SUM=B(LL)
+        B(LL)=B(I)
+        IF (II.NE.0)THEN
+          DO 11 J=II,I-1
+            SUM=SUM-A(I,J)*B(J)
+   11     CONTINUE
+        ELSE IF (SUM.NE.0.) THEN
+          II=I
+        ENDIF
+        B(I)=SUM
+   12 CONTINUE
+      DO 14 I=N,1,-1  
+        SUM=B(I)
+        IF(I.LT.N)THEN
+          DO 13 J=I+1,N
+            SUM=SUM-A(I,J)*B(J)
+   13 CONTINUE
+        ENDIF
+        B(I)=SUM/A(I,I)
+   14 CONTINUE
+      RETURN
+      END
+
+C> Replaces an NxN matrix a with the LU decomposition.
+C>
+C> @param A
+C> @param N
+C> @param NP
+C> @param INDX
+      SUBROUTINE LUDCMP(A,N,NP,INDX)
+C      PARAMETER (NMAX=400,TINY=1.0E-20)
+      PARAMETER (TINY=1.0E-20)
+C==EM==^^^
+C
+      REAL A(NP,NP),VV(N),D
+C     REAL A(NP,NP),VV(NMAX),D
+C==EM==^^^
+      INTEGER INDX(N)
+      D=1.  
+      DO 12 I=1,N
+        AAMAX=0.
+        DO 11 J=1,N
+          IF (ABS(A(I,J)).GT.AAMAX) AAMAX=ABS(A(I,J))
+   11   CONTINUE
+        IF (AAMAX.EQ.0.) print *, 'SINGULAR MATRIX.'
+        VV(I)=1./AAMAX
+   12 CONTINUE
+      DO 19 J=1,N
+        IF (J.GT.1) THEN
+          DO 14 I=1,J-1
+            SUM=A(I,J)
+            IF (I.GT.1)THEN
+              DO 13 K=1,I-1
+                SUM=SUM-A(I,K)*A(K,J)
+   13         CONTINUE
+              A(I,J)=SUM
+            ENDIF
+   14     CONTINUE
+        ENDIF
+        AAMAX=0.
+        DO 16 I=J,N
+          SUM=A(I,J)
+          IF (J.GT.1)THEN
+            DO 15 K=1,J-1
+              SUM=SUM-A(I,K)*A(K,J)
+   15       CONTINUE  
+            A(I,J)=SUM
+          ENDIF
+          DUM=VV(I)*ABS(SUM)
+          IF (DUM.GE.AAMAX) THEN
+            IMAX=I
+            AAMAX=DUM
+          ENDIF
+   16   CONTINUE
+        IF (J.NE.IMAX)THEN
+          DO 17 K=1,N
+            DUM=A(IMAX,K)
+            A(IMAX,K)=A(J,K)
+            A(J,K)=DUM
+   17     CONTINUE
+          D=-D
+          VV(IMAX)=VV(J)
+        ENDIF
+        INDX(J)=IMAX
+        IF(J.NE.N)THEN
+          IF(A(J,J).EQ.0.)A(J,J)=TINY
+          DUM=1./A(J,J)
+          DO 18 I=J+1,N
+            A(I,J)=A(I,J)*DUM
+   18     CONTINUE
+        ENDIF
+   19 CONTINUE
+      IF(A(N,N).EQ.0.)A(N,N)=TINY
+      RETURN
+      END
diff --git a/sorc/weight_gen.fd/splat.F b/sorc/weight_gen.fd/splat.F
@@ -0,0 +1,192 @@
+C> @file
+C> @brief Computes cosines of colatitude and Gaussian weights
+C> for sets of latitudes.
+C>
+C> ### Program History Log
+C> Date | Programmer | Comments
+C> -----|------------|---------
+C> 96-02-20 | Iredell | Initial.
+C> 97-10-20 | Iredell | Adjust precision.
+C> 98-06-11 | Iredell | Generalize precision using FORTRAN 90 intrinsic.
+C> 1998-12-03 | Iredell | Generalize precision further.
+C> 1998-12-03 | Iredell | Uses AIX ESSL BLAS calls.
+C> 2009-12-27 | D. Stark | Updated to switch between ESSL calls on an AIX platform, and Numerical Recipies calls elsewise.
+C> 2010-12-30 | Slovacek | Update alignment so preprocessor does not cause compilation failure.
+C> 2012-09-01 | E. Mirvis & M.Iredell | Merging & debugging linux errors of _d and _8 using generic LU factorization.
+C> 2012-11-05 | E. Mirvis | Generic FFTPACK and LU lapack were removed.
+C>
+C> @author Iredell @date 96-02-20
+
+C> Computes cosines of colatitude and Gaussian weights
+C> for one of the following specific global sets of latitudes.
+C> - Gaussian latitudes (IDRT=4)
+C> - Equally-spaced latitudes including poles (IDRT=0)
+C> - Equally-spaced latitudes excluding poles (IDRT=256)
+C>
+C> The Gaussian latitudes are located at the zeroes of the
+C> Legendre polynomial of the given order. These latitudes
+C> are efficient for reversible transforms from spectral space.
+C> (About twice as many equally-spaced latitudes are needed.)
+C> The weights for the equally-spaced latitudes are based on
+C> Ellsaesser (JAM,1966). (No weight is given the pole point.)
+C> Note that when analyzing grid to spectral in latitude pairs,
+C> if an equator point exists, its weight should be halved.
+C> This version invokes the ibm essl matrix solver.
+C>
+C> @param[in] IDRT grid identifier
+C>  - 4 for Gaussian grid
+C>  - 0 for equally-spaced grid including poles
+C>  - 256 for equally-spaced grid excluding poles
+C> @param[in] JMAX number of latitudes
+C> @param[out] SLAT sines of latitude
+C> @param[out] WLAT Gaussian weights
+C>
+C> @author Iredell @date 96-02-20
+      SUBROUTINE SPLAT(IDRT,JMAX,SLAT,WLAT)
+      REAL SLAT(JMAX),WLAT(JMAX)
+      INTEGER,PARAMETER:: KD=SELECTED_REAL_KIND(15,45)
+      REAL(KIND=KD):: PK(JMAX/2),PKM1(JMAX/2),PKM2(JMAX/2)
+      REAL(KIND=KD):: SLATD(JMAX/2),SP,SPMAX,EPS=10.*EPSILON(SP)
+      PARAMETER(JZ=50)
+      REAL BZ(JZ)
+      DATA BZ        / 2.4048255577,  5.5200781103,
+     $  8.6537279129, 11.7915344391, 14.9309177086, 18.0710639679,
+     $ 21.2116366299, 24.3524715308, 27.4934791320, 30.6346064684,
+     $ 33.7758202136, 36.9170983537, 40.0584257646, 43.1997917132,
+     $ 46.3411883717, 49.4826098974, 52.6240518411, 55.7655107550,
+     $ 58.9069839261, 62.0484691902, 65.1899648002, 68.3314693299,
+     $ 71.4729816036, 74.6145006437, 77.7560256304, 80.8975558711,
+     $ 84.0390907769, 87.1806298436, 90.3221726372, 93.4637187819,
+     $ 96.6052679510, 99.7468198587, 102.888374254, 106.029930916,
+     $ 109.171489649, 112.313050280, 115.454612653, 118.596176630,
+     $ 121.737742088, 124.879308913, 128.020877005, 131.162446275,
+     $ 134.304016638, 137.445588020, 140.587160352, 143.728733573,
+     $ 146.870307625, 150.011882457, 153.153458019, 156.295034268 /
+      REAL:: DLT,D1=1.
+      REAL AWORK((JMAX+1)/2,((JMAX+1)/2)),BWORK(((JMAX+1)/2))
+      INTEGER:: JHE,JHO,J0=0
+      INTEGER IPVT((JMAX+1)/2)
+      PARAMETER(PI=3.14159265358979,C=(1.-(2./PI)**2)*0.25)
+C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+C  GAUSSIAN LATITUDES
+      IF(IDRT.EQ.4) THEN
+        JH=JMAX/2
+        JHE=(JMAX+1)/2
+        R=1./SQRT((JMAX+0.5)**2+C)
+        DO J=1,MIN(JH,JZ)
+          SLATD(J)=COS(BZ(J)*R)
+        ENDDO
+        DO J=JZ+1,JH
+          SLATD(J)=COS((BZ(JZ)+(J-JZ)*PI)*R)
+        ENDDO
+        SPMAX=1.
+        DO WHILE(SPMAX.GT.EPS)
+          SPMAX=0.
+          DO J=1,JH
+            PKM1(J)=1.
+            PK(J)=SLATD(J)
+          ENDDO
+          DO N=2,JMAX
+            DO J=1,JH
+              PKM2(J)=PKM1(J)
+              PKM1(J)=PK(J)
+              PK(J)=((2*N-1)*SLATD(J)*PKM1(J)-(N-1)*PKM2(J))/N
+            ENDDO
+          ENDDO
+          DO J=1,JH
+            SP=PK(J)*(1.-SLATD(J)**2)/(JMAX*(PKM1(J)-SLATD(J)*PK(J)))
+            SLATD(J)=SLATD(J)-SP
+            SPMAX=MAX(SPMAX,ABS(SP))
+          ENDDO
+        ENDDO
+CDIR$ IVDEP
+        DO J=1,JH
+          SLAT(J)=SLATD(J)
+          WLAT(J)=(2.*(1.-SLATD(J)**2))/(JMAX*PKM1(J))**2
+          SLAT(JMAX+1-J)=-SLAT(J)
+          WLAT(JMAX+1-J)=WLAT(J)
+        ENDDO
+        IF(JHE.GT.JH) THEN
+          SLAT(JHE)=0.
+          WLAT(JHE)=2./JMAX**2
+          DO N=2,JMAX,2
+            WLAT(JHE)=WLAT(JHE)*N**2/(N-1)**2
+          ENDDO
+        ENDIF
+C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+C  EQUALLY-SPACED LATITUDES INCLUDING POLES
+      ELSEIF(IDRT.EQ.0) THEN
+        JH=JMAX/2
+        JHE=(JMAX+1)/2
+        JHO=JHE-1
+        DLT=PI/(JMAX-1)
+        SLAT(1)=1.
+        DO J=2,JH
+          SLAT(J)=COS((J-1)*DLT)
+        ENDDO
+        DO JS=1,JHO
+          DO J=1,JHO
+            AWORK(JS,J)=COS(2*(JS-1)*J*DLT)
+          ENDDO
+        ENDDO
+        DO JS=1,JHO
+          BWORK(JS)=-D1/(4*(JS-1)**2-1)
+        ENDDO
+
+        call ludcmp(awork,jho,jhe,ipvt)
+        call lubksb(awork,jho,jhe,ipvt,bwork)
+
+        WLAT(1)=0.
+        DO J=1,JHO
+          WLAT(J+1)=BWORK(J)
+        ENDDO
+CDIR$ IVDEP
+        DO J=1,JH
+           print *, j, jmax, JMAX+1-J
+          SLAT(JMAX+1-J)=-SLAT(J)
+          WLAT(JMAX+1-J)=WLAT(J)
+        ENDDO
+        IF(JHE.GT.JH) THEN
+          SLAT(JHE)=0.
+          WLAT(JHE)=2.*WLAT(JHE)
+        ENDIF
+C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+C  EQUALLY-SPACED LATITUDES EXCLUDING POLES
+      ELSEIF(IDRT.EQ.256) THEN
+        JH=JMAX/2
+        JHE=(JMAX+1)/2
+        JHO=JHE
+        DLT=PI/JMAX
+        SLAT(1)=1.
+        DO J=1,JH
+          SLAT(J)=COS((J-0.5)*DLT)
+        ENDDO
+        DO JS=1,JHO
+          DO J=1,JHO
+            AWORK(JS,J)=COS(2*(JS-1)*(J-0.5)*DLT)
+          ENDDO
+        ENDDO
+        DO JS=1,JHO
+          BWORK(JS)=-D1/(4*(JS-1)**2-1)
+        ENDDO
+
+        call ludcmp(awork,jho,jhe,ipvt,d)
+        call lubksb(awork,jho,jhe,ipvt,bwork)
+
+        WLAT(1)=0.
+        DO J=1,JHO
+          WLAT(J)=BWORK(J)
+        ENDDO
+CDIR$ IVDEP
+        DO J=1,JH
+          SLAT(JMAX+1-J)=-SLAT(J)
+          WLAT(JMAX+1-J)=WLAT(J)
+        ENDDO
+        IF(JHE.GT.JH) THEN
+          SLAT(JHE)=0.
+          WLAT(JHE)=2.*WLAT(JHE)
+        ENDIF
+      ENDIF
+C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+      RETURN
+      END