PROGRAM HYCOM_QUADLSQ IMPLICIT NONE C C hycom_quadlsq - Usage: hycom_quadlsq fin1.a fin2.a idm jdm nrec fout.a C C Outputs the quadratic squares linear fit of fin1 to fin2. C Output is s0, s1 and s2, for the fit: C fin2 ~= s0 + s1*fin1 + s2*fin1^2. C C fin*.a is assumed to contain idm*jdm 32-bit IEEE real values for C each array, in standard f77 element order, followed by padding C to a multiple of 4096 32-bit words, but otherwise with no control C bytes/words, and input values of 2.0**100 indicating a data void. C C The output will be a data void if either input file has a data void C at that location. C C this version for "serial" Unix systems. C C Alan J. Wallcraft, Naval Research Laboratory, July 2004. C REAL*4, ALLOCATABLE :: AX(:,:,:),AY(:,:,:) REAL*4 :: PAD(4096) INTEGER IOS,L INTEGER IARGC INTEGER NARG CHARACTER*240 CARG C INTEGER IDM,JDM,NREC,NPAD CHARACTER*240 CFILE1,CFILE2,CFILEO C C READ ARGUMENTS. C NARG = IARGC() C IF (NARG.EQ.6) THEN CALL GETARG(1,CFILE1) CALL GETARG(2,CFILE2) CALL GETARG(3,CARG) READ(CARG,*) IDM CALL GETARG(4,CARG) READ(CARG,*) JDM CALL GETARG(5,CARG) READ(CARG,*) NREC CALL GETARG(6,CFILEO) ELSE WRITE(6,'(3a)') & 'Usage: ', & 'hycom_quadlsq', & ' fin1.a fin2.a idm jdm nrec fout.a' CALL EXIT(1) ENDIF C NPAD = 4096 - MOD(IDM*JDM,4096) IF (NPAD.EQ.4096) THEN NPAD = 0 ENDIF C ALLOCATE( AX(IDM,JDM,NREC), STAT=IOS ) IF (IOS.NE.0) THEN WRITE(6,*) 'Error in hycom_quadlsq: could not allocate ', + IDM*JDM*NREC,' words' CALL EXIT(2) ENDIF ALLOCATE( AY(IDM,JDM,NREC), STAT=IOS ) IF (IOS.NE.0) THEN WRITE(6,*) 'Error in hycom_quadlsq: could not allocate ', + IDM*JDM*NREC,' words' CALL EXIT(2) ENDIF C CALL QUADLSQ(AX,AY,IDM,JDM,NREC,PAD,NPAD, & CFILE1,CFILE2,CFILEO) CALL EXIT(0) END SUBROUTINE QUADLSQ(AX,AY,IDM,JDM,NREC,PAD,NPAD, & CFILE1,CFILE2,CFILEO) IMPLICIT NONE C REAL*4 SPVAL PARAMETER (SPVAL=2.0**100) C CHARACTER*240 CFILE1,CFILE2,CFILEO INTEGER IDM,JDM,NREC,NPAD REAL*4 AX(IDM,JDM,NREC), & AY(IDM,JDM,NREC),PAD(NPAD) C C MOST OF WORK IS DONE HERE. C #ifdef sun INTEGER IR_ISNAN C #endif CHARACTER*18 CASN INTEGER I,J,IOS,IR,K,NRECL REAL X(NREC),Y(NREC),C(3),CMAX(3),CMIN(3) #ifdef CRAY INTEGER*8 IU8,IOS8 #endif C INQUIRE( IOLENGTH=NRECL) AX(:,:,1),PAD #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 = 12 CALL ASNUNIT(IU8,CASN,IOS8) IF (IOS8.NE.0) THEN write(6,*) 'Error: can''t asnunit 12' 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(12,'-F syscall -N ieee',IOS) IF (IOS.NE.0) THEN write(6,*) 'Error: can''t asnunit 12' 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 OPEN(UNIT=11, FILE=CFILE1, FORM='UNFORMATTED', STATUS='OLD', + ACCESS='DIRECT', RECL=NRECL, IOSTAT=IOS) IF (IOS.NE.0) THEN write(6,*) 'Error: can''t open ',TRIM(CFILE1) write(6,*) 'ios = ',ios write(6,*) 'nrecl = ',nrecl CALL EXIT(3) ENDIF OPEN(UNIT=12, FILE=CFILE2, FORM='UNFORMATTED', STATUS='OLD', + ACCESS='DIRECT', RECL=NRECL, IOSTAT=IOS) IF (IOS.NE.0) THEN write(6,*) 'Error: can''t open ',TRIM(CFILE2) write(6,*) 'ios = ',ios write(6,*) 'nrecl = ',nrecl CALL EXIT(3) ENDIF OPEN(UNIT=21, FILE=CFILEO, FORM='UNFORMATTED', STATUS='NEW', + ACCESS='DIRECT', RECL=NRECL, IOSTAT=IOS) IF (IOS.NE.0) THEN write(6,*) 'Error: can''t open ',TRIM(CFILEO) write(6,*) 'ios = ',ios write(6,*) 'nrecl = ',nrecl CALL EXIT(3) ENDIF C DO IR=1,NREC READ(11,REC=IR,IOSTAT=IOS) AX(:,:,IR) #ifdef ENDIAN_IO CALL ENDIAN_SWAP(AX(1,1,IR),IDM*JDM) #endif IF (IOS.NE.0) THEN WRITE(6,*) 'can''t read record ',IR,' of ', & TRIM(CFILE1) CALL EXIT(4) ENDIF READ(12,REC=IR,IOSTAT=IOS) AY(:,:,IR) #ifdef ENDIAN_IO CALL ENDIAN_SWAP(AY(1,1,IR),IDM*JDM) #endif IF (IOS.NE.0) THEN WRITE(6,*) 'can''t read record ',IR,' of ', & TRIM(CFILE2) CALL EXIT(4) ENDIF ENDDO !ir C C QUADRATIC FIT. C CMAX(:) = -SPVAL CMIN(:) = SPVAL DO J= 1,JDM DO I= 1,IDM IF (AX(I,J,1).NE.SPVAL) THEN DO IR= 1,NREC X(IR) = AX(I,J,IR) Y(IR) = AY(I,J,IR) ENDDO !ir CALL QUADFIT(X,Y,NREC,C) DO K= 1,3 AX(I,J,K) = C(K) CMIN(K) = MIN( CMIN(K), C(K) ) CMAX(K) = MAX( CMAX(K), C(K) ) ENDDO !k ELSE AX(I,J,1) = SPVAL AX(I,J,2) = SPVAL AX(I,J,3) = SPVAL ENDIF ENDDO !i ENDDO !j C C OUTPUT THE RESULT. C WRITE(21,REC=1,IOSTAT=IOS) AX(:,:,1) IF (IOS.NE.0) THEN write(6,*) 'Error: can''t write to ',TRIM(CFILEO) write(6,*) 'ios = ',ios write(6,*) 'rec = ',1 CALL EXIT(3) ENDIF WRITE(6,'(a,1p2g16.8)') ' c0: min, max = ',CMIN(1),CMAX(1) C WRITE(21,REC=2,IOSTAT=IOS) AX(:,:,2) IF (IOS.NE.0) THEN write(6,*) 'Error: can''t write to ',TRIM(CFILEO) write(6,*) 'ios = ',ios write(6,*) 'rec = ',2 CALL EXIT(3) ENDIF WRITE(6,'(a,1p2g16.8)') ' c1: min, max = ',CMIN(2),CMAX(2) C WRITE(21,REC=3,IOSTAT=IOS) AX(:,:,3) IF (IOS.NE.0) THEN write(6,*) 'Error: can''t write to ',TRIM(CFILEO) write(6,*) 'ios = ',ios write(6,*) 'rec = ',3 CALL EXIT(3) ENDIF WRITE(6,'(a,1p2g16.8)') ' c2: min, max = ',CMIN(3),CMAX(3) C CLOSE(UNIT=11) CLOSE(UNIT=12) CLOSE(UNIT=21) C RETURN END subroutine quadfit (xd,yd,ndata,c) implicit none c integer ndata real xd(ndata),yd(ndata),c(3) c c Take corresponding data points from the arrays xd and yd and fit c them with the following equation: c c y = c(1) + c(2) * x + c(3) * x**2 c real aa(3,3),fs,suma,sumb,sqrt,rsid integer ipvt(3),ir,is,ij, info c c Input Arguments: c xd - x values for data points c yd - y valuess for data point pairs c ndata - number of data points c c Output Arguements: c c - array containing the three coefficients in the 2nd order c polynomial that provide the "best" fit to the data from a c least squares method. Note that it also temporarily holds c the values of the right hand sides of each Least Squares c equation. c c Other Key Variables: c aa - matrix containing coefficients of the system of equations c generated by the least squares method c suma - a variable that tallies the sum that is needed to generate c each element of aa. c sumb - a variable that tallies the sum that is needed for the right c hand side of each Least Squares equation. c ir - an index that is used to keep track of the equation number c within the system of equations . c is - an index used to track the coefficient number within a given c Least Squares equation. c c c DO loop 55 generates terms for each of the 3 Least Squares Equations c do 55 ir=1,3 c c DO loop 45 generates the right hand side for a given equation c sumb=0. do ij=1,ndata ! loop 45 sumb=sumb+yd(ij)*xd(ij)**(ir-1) enddo c(ir)=sumb c c DO loop 50 generates the coefficients a given equation c do is=1,3 ! loop 50 suma=0. do ij=1,ndata suma=suma+xd(ij)**(is-1)*xd(ij)**(ir-1) enddo aa(ir,is)=suma enddo 55 continue c c Solve the Least Squares Equations c call sgefa(aa,3,3,ipvt,info) call sgesl(aa,3,3,ipvt,c ,0) *c *c Calculate a measure of mean error between the data and the curve *c * rsid=0. * do 65 ir=1,ndata * fs=0. * do 60 is=1,3 * 60 fs=fs+c(is)*xd(ir)**(is-1) * 65 rsid=rsid+(yd(ir)-fs)**2 * rsid=sqrt(rsid/float(ndata-1)) * write(6,2001) rsid * 2001 format(' Fit to 2nd order polynomial has a mean error of', * $ 1p,e12.5) return end subroutine sgefa(a,lda,n,ipvt,info) integer lda,n,ipvt(1),info real a(lda,1) c c sgefa factors a real matrix by gaussian elimination. c c sgefa is usually called by sgeco, but it can be called c directly with a saving in time if rcond is not needed. c (time for sgeco) = (1 + 9/n)*(time for sgefa) . c c on entry c c a real(lda, n) c the matrix to be factored. c c lda integer c the leading dimension of the array a . c c n integer c the order of the matrix a . c c on return c c a an upper triangular matrix and the multipliers c which were used to obtain it. c the factorization can be written a = l*u where c l is a product of permutation and unit lower c triangular matrices and u is upper triangular. c c ipvt integer(n) c an integer vector of pivot indices. c c info integer c = 0 normal value. c = k if u(k,k) .eq. 0.0 . this is not an error c condition for this subroutine, but it does c indicate that sgesl or sgedi will divide by zero c if called. use rcond in sgeco for a reliable c indication of singularity. c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c blas saxpy,sscal,isamax c c internal variables c real t integer isamax,j,k,kp1,l,nm1 c c c gaussian elimination with partial pivoting c info = 0 nm1 = n - 1 if (nm1 .lt. 1) go to 70 do 60 k = 1, nm1 kp1 = k + 1 c c find l = pivot index c l = isamax(n-k+1,a(k,k),1) + k - 1 ipvt(k) = l c c zero pivot implies this column already triangularized c if (a(l,k) .eq. 0.0e0) go to 40 c c interchange if necessary c if (l .eq. k) go to 10 t = a(l,k) a(l,k) = a(k,k) a(k,k) = t 10 continue c c compute multipliers c t = -1.0e0/a(k,k) call sscal(n-k,t,a(k+1,k),1) c c row elimination with column indexing c do 30 j = kp1, n t = a(l,j) if (l .eq. k) go to 20 a(l,j) = a(k,j) a(k,j) = t 20 continue call saxpy(n-k,t,a(k+1,k),1,a(k+1,j),1) 30 continue go to 50 40 continue info = k 50 continue 60 continue 70 continue ipvt(n) = n if (a(n,n) .eq. 0.0e0) info = n return end subroutine sgesl(a,lda,n,ipvt,b,job) integer lda,n,ipvt(1),job real a(lda,1),b(1) c c sgesl solves the real system c a * x = b or trans(a) * x = b c using the factors computed by sgeco or sgefa. c c on entry c c a real(lda, n) c the output from sgeco or sgefa. c c lda integer c the leading dimension of the array a . c c n integer c the order of the matrix a . c c ipvt integer(n) c the pivot vector from sgeco or sgefa. c c b real(n) c the right hand side vector. c c job integer c = 0 to solve a*x = b , c = nonzero to solve trans(a)*x = b where c trans(a) is the transpose. c c on return c c b the solution vector x . c c error condition c c a division by zero will occur if the input factor contains a c zero on the diagonal. technically this indicates singularity c but it is often caused by improper arguments or improper c setting of lda . it will not occur if the subroutines are c called correctly and if sgeco has set rcond .gt. 0.0 c or sgefa has set info .eq. 0 . c c to compute inverse(a) * c where c is a matrix c with p columns c call sgeco(a,lda,n,ipvt,rcond,z) c if (rcond is too small) go to ... c do 10 j = 1, p c call sgesl(a,lda,n,ipvt,c(1,j),0) c 10 continue c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c blas saxpy,sdot c c internal variables c real sdot,t integer k,kb,l,nm1 c nm1 = n - 1 if (job .ne. 0) go to 50 c c job = 0 , solve a * x = b c first solve l*y = b c if (nm1 .lt. 1) go to 30 do 20 k = 1, nm1 l = ipvt(k) t = b(l) if (l .eq. k) go to 10 b(l) = b(k) b(k) = t 10 continue call saxpy(n-k,t,a(k+1,k),1,b(k+1),1) 20 continue 30 continue c c now solve u*x = y c do 40 kb = 1, n k = n + 1 - kb b(k) = b(k)/a(k,k) t = -b(k) call saxpy(k-1,t,a(1,k),1,b(1),1) 40 continue go to 100 50 continue c c job = nonzero, solve trans(a) * x = b c first solve trans(u)*y = b c do 60 k = 1, n t = sdot(k-1,a(1,k),1,b(1),1) b(k) = (b(k) - t)/a(k,k) 60 continue c c now solve trans(l)*x = y c if (nm1 .lt. 1) go to 90 do 80 kb = 1, nm1 k = n - kb b(k) = b(k) + sdot(n-k,a(k+1,k),1,b(k+1),1) l = ipvt(k) if (l .eq. k) go to 70 t = b(l) b(l) = b(k) b(k) = t 70 continue 80 continue 90 continue 100 continue return end integer function isamax(n,sx,incx) c c finds the index of element having max. absolute value. c jack dongarra, linpack, 3/11/78. c modified 3/93 to return if incx .le. 0. c modified 12/3/93, array(1) declarations changed to array(*) c real sx(*),smax integer i,incx,ix,n c isamax = 0 if( n.lt.1 .or. incx.le.0 ) return isamax = 1 if(n.eq.1)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 smax = abs(sx(1)) ix = ix + incx do 10 i = 2,n if(abs(sx(ix)).le.smax) go to 5 isamax = i smax = abs(sx(ix)) 5 ix = ix + incx 10 continue return c c code for increment equal to 1 c 20 smax = abs(sx(1)) do 30 i = 2,n if(abs(sx(i)).le.smax) go to 30 isamax = i smax = abs(sx(i)) 30 continue return end subroutine saxpy(n,sa,sx,incx,sy,incy) c c constant times a vector plus a vector. c uses unrolled loop for increments equal to one. c jack dongarra, linpack, 3/11/78. c modified 12/3/93, array(1) declarations changed to array(*) c real sx(*),sy(*),sa integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if (sa .eq. 0.0) return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n sy(iy) = sy(iy) + sa*sx(ix) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,4) if( m .eq. 0 ) go to 40 do 30 i = 1,m sy(i) = sy(i) + sa*sx(i) 30 continue if( n .lt. 4 ) return 40 mp1 = m + 1 do 50 i = mp1,n,4 sy(i) = sy(i) + sa*sx(i) sy(i + 1) = sy(i + 1) + sa*sx(i + 1) sy(i + 2) = sy(i + 2) + sa*sx(i + 2) sy(i + 3) = sy(i + 3) + sa*sx(i + 3) 50 continue return end real function sdot(n,sx,incx,sy,incy) c c forms the dot product of two vectors. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c modified 12/3/93, array(1) declarations changed to array(*) c real sx(*),sy(*),stemp integer i,incx,incy,ix,iy,m,mp1,n c stemp = 0.0e0 sdot = 0.0e0 if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n stemp = stemp + sx(ix)*sy(iy) ix = ix + incx iy = iy + incy 10 continue sdot = stemp return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m stemp = stemp + sx(i)*sy(i) 30 continue if( n .lt. 5 ) go to 60 40 mp1 = m + 1 do 50 i = mp1,n,5 stemp = stemp + sx(i)*sy(i) + sx(i + 1)*sy(i + 1) + * sx(i + 2)*sy(i + 2) + sx(i + 3)*sy(i + 3) + sx(i + 4)*sy(i + 4) 50 continue 60 sdot = stemp return end subroutine sscal(n,sa,sx,incx) c c scales a vector by a constant. c uses unrolled loops for increment equal to 1. c jack dongarra, linpack, 3/11/78. c modified 3/93 to return if incx .le. 0. c modified 12/3/93, array(1) declarations changed to array(*) c real sa,sx(*) integer i,incx,m,mp1,n,nincx c if( n.le.0 .or. incx.le.0 )return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c nincx = n*incx do 10 i = 1,nincx,incx sx(i) = sa*sx(i) 10 continue return c c code for increment equal to 1 c c c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m sx(i) = sa*sx(i) 30 continue if( n .lt. 5 ) return 40 mp1 = m + 1 do 50 i = mp1,n,5 sx(i) = sa*sx(i) sx(i + 1) = sa*sx(i + 1) sx(i + 2) = sa*sx(i + 2) sx(i + 3) = sa*sx(i + 3) sx(i + 4) = sa*sx(i + 4) 50 continue return end