cat > aqclib.f <<\EOR C=============================================================================C C International Comprehensive Ocean-Atmosphere Data Set (I-COADS) 10 Sep 2002 C C Filename:level: aqclib:01B Fortran 77 program+shell C C Function: Tools to assist quality control Author: Tom Smith et al. C C=============================================================================C C Software Revision Information (previous version: 12 Aug 2002): Updated C comment structure for Release 2.0 and I-COADS. Changes from level 01A C documented in individual routines. C-----------------------------------------------------------------------3456789 C Software documentation for the (modifiable) library test routine {libtst}, C and for the remaining (invariant) user-interface routines in the library. C C Functionality: This is a library of tools for use with quality control, C as described at the beginning of each subroutine or function. C C Data structures: C The 2-dim spatial domain (180,90) is organized with 180 (90) 2-deg boxes C for each latitude (longitude) band. Box central longitudes are located at: C 1E, 3E, 5E, ..., 359E for i=1,2,3, ..., 180 C box central latitudes are located at: C 89N, 87N, ..., 89S for j=1,2, ..., 90. C C Content: Following are the routines included: C {bchkcl} buddy checking C {bi11s} binomial smoothing/filling (11-pt filter w/ wrap) C {bi3s} binomial smoothing/filling (3-pt filter w/o wrap) C {binomf} binomial filling over the spatial domain C {fmed} function to obtain median C {fslsa} Fourier series estimate C {mask2d} 2-deg landmask (based on COADS Release 1 LLN2F1) C C Machine dependencies: None known. C For more information: See and (electronic documents). C-----------------------------------------------------------------------3456789 c program libtst cc subroutine libtst cc-----temporary code for testing purposes c implicit integer(a-e,g-z) c end C=============================================================================C C WARNING: Code beyond this point should not require any modification. C C=============================================================================C C-----------------------------------------------------------------------3456789 subroutine bchkcl(clim,mask) c-----Climatology buddy checking (using data from the spatial/temporal c neighborhood). The buddy checking is used to remove isolated c extreme values in 12-month climatology. The buddy checking includes c following processes: c (1) Find the min, max, and average of the climatological values of a c cube consisting of 27 2-deg boxes (3 boxes in longitude x 3 boxes c in latitude x 3 boxes in time; ref. Release 1, p. C4. Note: the c true min/max are discarded and the next nearest min/max values used. c (2) Remove isolated values by setting a 2-deg box value to missing if c the cube centered on this 2-deg box contains less than 5 extant c 2-deg box values. c (3) Remove isolated extreme 2-deg box values. For each 2-deg box value, c if (clim(i,j,m)-ave)>(max-min), then the 2-deg box value clim(i,j,m) c is set to missing. c Arguments (using the 2-dim spatial domain, plus month): c Input: real clim(180,90,12): climatology c integer mask(180,90): landmask c Output: real clim(180,90,12): checked climatology c-----Tom Smith, 2001. parameter(im=180,jm=90,maxx=5000) real clim(im,jm,12),wk(im,jm,12),da1(maxx) integer mask(im,jm) do 1 mon=1,12 do 1 j=1,jm do 1 i=1,im wk(i,j,mon)=clim(i,j,mon) 1 continue c do mon=1,12 do 3 j=1,jm j1=j-1 if(j1.lt.1) j1=1 j2=j+1 if(j2.gt.jm)j2=jm do 3 i=1,im if(mask(i,j).eq.0)then clim(i,j,mon)=-999.9 go to 3 endif if(clim(i,j,mon).lt.-990.) go to 3 kn=0 do 2 mm=mon-1,mon+1 md=mm if(mm.lt.1) md=12 if(mm.gt.12)md=1 do 2 jj=j1,j2 do 2 ii=i-1,i+1 id=ii if(ii.lt.1) id=ii+im if(ii.gt.im)id=ii-im if(wk(id,jj,md).gt.-990.)then kn=kn+1 da1(kn)=wk(id,jj,md) if(mm.eq.mon)then if(ii.eq.i)then if(jj.eq.j) kn=kn-1 endif endif endif 2 continue if(kn.lt.5)then clim(i,j,mon)=-999.9 go to 3 endif amin=1.0e9 amax=-amin do k=1,kn if(da1(k).lt.amin)then amin=da1(k) kmin=k endif if(da1(k).gt.amax)then amax=da1(k) kmax=k endif end do s1=0.0 s2=0.0 amin=1.0e9 amax=-amin do k=1,kn if(k.ne.kmin)then if(k.ne.kmax)then s1=s1+da1(k) s2=s2+1.0 if(da1(k).lt.amin) amin=da1(k) if(da1(k).gt.amax) amax=da1(k) endif endif end do avg=s1/s2 civ=(amax-amin) if(abs(clim(i,j,mon)-avg).gt.civ) clim(i,j,mon)=-999.9 3 continue end do return end C-----------------------------------------------------------------------3456789 subroutine bi11s(a1,ntm) c-----Binomial smoothing/filling over the ntm elements of a1 (if sum c or binomial weights used is at least 0.5). An 11-point filter is c used with wrapping, i.e., (a1(ntm+1) = a1(1)). c Arguments: c Input: real a1(ntm): series of ntm values c integer ntm: length of a1 c Output: real a1(ntm): smoothed/filled series c-----Tom Smith, 2001. real a1(*),wk(500),wgt(11) data wgt/0.001,0.010,0.044,0.117,0.205,0.246, 2 0.205,0.117,0.044,0.010,0.001/ do n=1,ntm wk(n)=a1(n) end do do np=ntm+1,ntm+11 n=np-ntm wk(np)=a1(n) end do do n=1,5 k1=n-5 k2=n+5 s1=0.0 s2=0.0 kn=0 do k=k1,k2 kn=kn+1 kd=k if(k.lt.1)kd=k+ntm if(wk(kd).gt.-990.)then s1=s1+wgt(kn)*wk(kd) s2=s2+wgt(kn) endif end do if(s2.ge.0.499) a1(n)=s1/s2 end do c do n=6,ntm k1=n-5 k2=n+5 s1=0.0 s2=0.0 kn=0 do k=k1,k2 kn=kn+1 if(wk(k).gt.-990.)then s1=s1+wgt(kn)*wk(k) s2=s2+wgt(kn) endif end do if(s2.ge.0.499) a1(n)=s1/s2 end do return end C-----------------------------------------------------------------------3456789 subroutine bi3s(a1,ntm) c-----Binomial smoothing/filling over the ntm elements of a1 (if sum of c binomial weights used is at least 0.5). A 3-point filter is used c without wrapping. c Arguments: c Input: real a1(ntm): series of ntm values c integer ntm: length of a1 c Output: real a1(ntm): smoothed/filled series c-----Tom Smith, 2001. real a1(*),wk(500),wgt(3) data wgt/0.25,0.50,0.25/ do n=1,ntm wk(n)=a1(n) end do s1=0.0 s2=0.0 if(wk(1).gt.-990.)then s1=s1+wgt(2)*wk(1) s2=s2+wgt(2) if(wk(2).gt.-990.)then s1=s1+wgt(3)*wk(2) s2=s2+wgt(3) endif a1(1)=s1/s2 endif s1=0.0 s2=0.0 if(wk(ntm).gt.-990.)then s1=s1+wgt(2)*wk(ntm) s2=s2+wgt(2) if(wk(ntm-1).gt.-990.)then s1=s1+wgt(3)*wk(ntm-1) s2=s2+wgt(3) endif a1(ntm)=s1/s2 endif c do n=2,ntm-1 k1=n-1 k2=n+1 s1=0.0 s2=0.0 kn=0 do k=k1,k2 kn=kn+1 if(wk(k).gt.-990.)then s1=s1+wgt(kn)*wk(k) s2=s2+wgt(kn) endif end do if(s2.ge.0.499) a1(n)=s1/s2 end do return end C-----------------------------------------------------------------------3456789 subroutine binomf(a2,mask) c-----Binomial filling over the spatial domain. An 11-point filter is c used with wrapping (a1(ntm+1) = a1(1)) zonally, and a 3-point filter c without wrapping is used meridionally. Existing values are not c replaced. c Arguments (using the 2-dim spatial domain): c Input: real a2(180,90): data array c integer mask(180,90): landmask c Output: real a2(180,90): filled data array c-----Tom Smith, 2001. parameter(im=180,jm=90,jmm1=jm-1) real a2(im,jm),wk(im,jm),wz(11),wm(3) integer mask(im,jm) data wz/0.001,0.010,0.044,0.117,0.205,0.246, 2 0.205,0.117,0.044,0.010,0.001/, wm/0.25,0.50,0.25/ iter=0 10 knu=0 do 1 j=1,jm do 1 i=1,im wk(i,j)=a2(i,j) if(a2(i,j).le.-990.) knu=knu+1 1 continue if(knu.eq.0)then do 2 j=1,jm do 2 i=1,im if(mask(i,j).eq.0) a2(i,j)=-999.9 2 continue return endif iter=iter+1 c do 4 j=2,jmm1 jj1=j-1 jj2=j+1 do 4 i=1,im if(a2(i,j).gt.-990.) go to 4 ii1=i-5 ii2=i+5 s1=0.0 s2=0.0 knj=0 do 3 jj=jj1,jj2 knj=knj+1 kni=0 do 3 ii=ii1,ii2 kni=kni+1 id=ii if(id.lt.1) id=ii+im if(id.gt.im)id=ii-im if(wk(id,jj).gt.-990.)then ww=wz(kni)*wm(knj) s1=s1+ww*wk(id,jj) s2=s2+ww endif 3 continue if(s2.ge.0.1) a2(i,j)=s1/s2 4 continue do i=1,im if(a2(i,2).gt.-990.) a2(i,1)=a2(i,2) if(a2(i,jmm1).gt.-990.)a2(i,jm)=a2(i,jmm1) end do go to 10 end C-----------------------------------------------------------------------3456789 function fmed(ts,ntm) c-----fmed returns the median of the ntm values of ts. ntm can be of c any length greater than or equal to 1. c Arguments (using the 2-dim spatial domain): c Input: real ts(ntm): series of ntm values c integer ntm: the length of ts c Output: real fmed: the median value c-----Tom Smith, 2001. real ts(*) if(ntm.gt.1)then do i=1,ntm-1 do j=i+1,ntm if(ts(j).lt.ts(i))then tmp=ts(i) ts(i)=ts(j) ts(j)=tmp endif end do end do endif if(mod(ntm,2).eq.0)then n1=ntm/2 n2=n1+1 fmed=(ts(n1)+ts(n2))/2.0 else n0=(ntm+1)/2 fmed=ts(n0) endif return end C-----------------------------------------------------------------------3456789 subroutine fslsa(f,m) C-----Fourier series estimate from an annual-cycle series with missing c values. L.S. fitting is used to compute the harmonics. An exact c Fourier series is computed if all values are filled. Missing c values should be assigned a value of -999.9 to be skipped. c The harmonic estimate of f(t) is: c f'(t) = x0 + {sum,k=1,m,[an(k)*sin(k*t*tpi/p)+bn(k)*cos(k*t*tpi/p)]} c Arguments: c Input: real f(nt): series of nt values c integer nt: the length of f c Output: real f(nt): smoothed/filled series c integer m: the number of modes computed c-----Tom Smith, 2001. c-----S. Woodruff, 6 Sep 2002: changed pi to 36 from 9 digits to match c {lmrlib} precision. c-----S. Lubker, 9 Sep 2002: changed to use implicit double precision. implicit double precision(a-h,o-z) parameter(max=200,maxt=5000, & pi=3.14159265358979323846264338327950288d0, & tpi=2*pi,crit=0.05) real f(*),phi(max,maxt)*8 double precision a(max,max),b(max) integer indx(max) mdo=m x0=0.0 c* Check the number of modes, max < nt/2. nchk=12 mmx=5 if(m.gt.mmx)mdo=mmx c* Get the mean, and check to see the maximum number of modes that c* may be supported by the sampling. s1=0.0 s2=0.0 kn=0 do jt=1,12 if(f(jt).gt.-990.)then s1=s1+f(jt) s2=s2+1.0 kn=kn+1 endif end do if(s2.gt.0.0) x0=s1/s2 c if(mod(kn,2).eq.0)then mmx=kn/2 else mmx=(kn-1)/2 endif if(mdo.gt.mmx)then mdo=mmx if(mmx.eq.0) return endif c* Check the number of modes vs max. nn=2*mdo if(nn.gt.max)then c print*, 'Too many modes requested, reducing m.' mdo=max/2 nn=max endif c* Evaluate the function for all times and all modes. do jt=1,12 do k=1,mdo arg=tpi*real(k*jt)/12.0 ks=2*k-1 kc=2*k phi(ks,jt)=dsin(arg) phi(kc,jt)=dcos(arg) end do end do c* Check sampling, truncate modes higher than mode k if sampling is c* not adequate for either the sin or the cos component of mode k+1. mmx=0 do k=1,mdo ks=2*k-1 kc=2*k s1=0.0 s2=0.0 c1=0.0 c2=0.0 do jt=1,12 vs=phi(ks,jt)*phi(ks,jt) vc=phi(kc,jt)*phi(kc,jt) if(f(jt).gt.-990.)then s1=s1+vs c1=c1+vc endif s2=s2+vs c2=c2+vc end do vs=s1/s2 vc=c1/c2 vv=vs if(vc.gt.vs) vv=vc if(vv.lt.crit) go to 10 mmx=mmx+1 end do 10 continue if(mmx.lt.mdo)then if(mmx.eq.0) return mdo=mmx nn=2*mdo endif c* Set up the a and b arrays for the normal equations, using recentered c* data (remove the mean, x0). do i=1,nn do j=1,i sum=0.0 do jt=1,12 if(f(jt).gt.-990.) sum=sum+phi(i,jt)*phi(j,jt) end do a(i,j)=sum a(j,i)=sum end do sum=0.0 do jt=1,12 if(f(jt).gt.-990.) sum=sum+(f(jt)-x0)*phi(i,jt) end do b(i)=sum end do c* Invert a to solve the system Aw=b, return w in b. call ludcmp(a,nn,max,indx,d) call lubksb(a,nn,max,indx,b) c* Fill with the smoothed values and return. do jt=1,12 sum=x0 do k=1,mdo ks=2*k-1 kc=2*k sum=sum+b(ks)*phi(ks,jt)+b(kc)*phi(kc,jt) end do f(jt)=sum end do return end C-----------------------------------------------------------------------3456789 subroutine ludcmp(a,n,np,indx,d) c-----Does an LU decomposition on the matrix A, size nxn. The output is c returned in A. This is Crout's method with partial pivoting, c described in Numerical Recipes, sec 2.3. c-----Tom Smith, 2001. c-----S. Lubker, 9 Sep 2002: changed to use implicit double precision. implicit double precision(a-h,o-z) parameter (nmax=200,tiny=1.0e-20) dimension indx(*),vv(nmax) double precision a(np,np) 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.) pause '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 C-----------------------------------------------------------------------3456789 subroutine lubksb(a,n,np,indx,b) c-----Solves the system Ax=b, where here A is the LU decomposition of c the matrix A that was passed into ludcmp() c-----Tom Smith, 2001. c-----S. Lubker, 9 Sep 2002: changed to use implicit double precision. implicit double precision(a-h,o-z) dimension indx(*) double precision a(np,np),b(*) 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-----------------------------------------------------------------------3456789 subroutine mask2d(mask) c-----Initialize mask, the 2-deg landmask, such that 0=land and 1=ocean. c Reconfigures a version (lacking polar boxes) of the COADS Release c 1 (Slutz et al., 1985, supp. G) 2-deg landmask LLN2F1 to match the c grid used in this library. c-----S. Lubker, 9 September 2002. c Reference: c Slutz, R.J., S.J. Lubker, J.D. Hiscox, S.D. Woodruff, R.L. Jenne, D.H. c Joseph, P.M. Steurer, and J.D. Elms, 1985: Comprehensive Ocean- c Atmosphere Data Set; Release 1. NOAA Environmental Research c Laboratories, Climate Research Program, Boulder, Colo., 268 pp. integer mask(180,90) common /lln2/lln2f1(180,90) save c do 190 j=1,90 do 180 i=1,180 mask(i,j)=1-mod(lln2f1(i,91-j),2) 180 continue 190 continue end C-----------------------------------------------------------------------3456789 block data bdlln2 c-----COADS Release 1 (Slutz et al., 1985, supp. G) 2-deg landmask LLN2F1 c (lacking polar boxes). 0=ocean, 1=land, and 2=coast. c-----S. Lubker, 9 September 2002. c Reference: c Slutz, R.J., S.J. Lubker, J.D. Hiscox, S.D. Woodruff, R.L. Jenne, D.H. c Joseph, P.M. Steurer, and J.D. Elms, 1985: Comprehensive Ocean- c Atmosphere Data Set; Release 1. NOAA Environmental Research c Laboratories, Climate Research Program, Boulder, Colo., 268 pp. common /lln2/lln2f1(180,90) c data ((lln2f1(i,j),i=1,180),j=1,15) +/680*0,5*2,15*0,3*2,150*0,2,3*0,3*2,2*0,2,1,2,3*0,3*2,2*1,2*2,1 +,3*2,4*1,7*2,58*0,2,85*0,2*2,2*0,3*2,5*0,2*2,21*1,2,61*0,2,95*0 +,3*2,18*1,2,55*0,6*2,4*1,2*2,13*0,2,81*0,2*2,15*1,2,38*0,2,13*0 +,3*2,7*1,2*2,0,6*2,58*0,3*2,0,2,3*0,2,3*0,2,2*0,2,0,2*2,0,2,14*0,2 +,11*1,2,48*0,2*2,0,2,0,2*2,0,9*1,2,10*1,2,0,2,0,5*2,1,3*2,47*0,6*2 +,4*0,2,5*0,2*2,0,3*2,11*0,10*1,2,23*0,3*2,2*1,4*2,13*0,2*2,0,2,1 +,2*0,2*2,39*1,11*2,8*0,2,8*1,14*2,2*0,3*2,7*0,2*2,3*0,2,4*0,2,1 +,10*0,8*1,4*2,20*0,2*2,6*1,5*2,4*0,2,1,2*2,6*1,2*2,0,3*2,34*1,2 +,15*1,3*2,6*0,2,22*1,11*2,1,5*2,5*0,1,2*2,8*0,6*1,2*2,23*0,2,3*1,2 +,2*0,4*1,3*0,2,1,2,65*1,2*2,2*0,2*2,3*0,3*2,34*1,2,2*0,2,5*0,2,1 +,2*2,8*0,2,3*1,11*0,3*2,11*0,2*2,2*1,2,2*0,2*2,5*1,2*2,61*1,2*2 +,4*1,2*2,9*0,2*2,33*1,2,11*0,2,9*0,2,2*1,2,25*0,4*1,2,2*0,4*2,0,2 +,60*1,2,4*0,4*2,11*0,3*2,3*1,2,3*0,3*2,22*1,2,8*0,2,1,2*2,12*0,2 +,22*0/ c data ((lln2f1(i,j),i=1,180),j=16,30) +/3*0,2*2,0,2*1,2,3*0,3*2,55*1,2,9*0,2,19*0,3*2,9*0,2*2,18*1,2,8*0 +,2,3*1,2*0,2,38*0,2*2,3*0,2*2,55*1,2,9*0,3*2,33*0,2,18*1,2*2,7*0 +,3*1,2,2*1,2,42*0,56*1,2,10*0,1,2*2,34*0,2,20*1,3*2,2*0,2,7*1,2 +,26*0,2,6*0,5*2,60*1,2*2,0,2,6*0,2,37*0,2,21*1,2,2*0,8*1,3*2,23*0 +,2*2,0,2*2,0,2*2,67*1,47*0,2*2,20*1,2*2,5*1,5*2,28*0,2*2,68*1,2 +,49*0,16*1,3*2,7*1,2,5*0,2*2,26*0,2,15*1,4*2,0,3*1,4*2,41*1,2,49*0 +,2,14*1,2,4*0,2*2,3*1,4*2,32*0,2,4*1,2*2,0,2,6*1,5*0,2,3*1,2,2*0 +,2*2,39*1,2,50*0,2,14*1,3*2,0,2,2*0,4*2,1,2,34*0,1,2,3*0,2*2,2*0,2 +,3*1,2,6*0,2,2*1,2,2*0,2,38*1,2*2,4*0,2,46*0,16*1,2,1,0,1,0,2,0,2 +,2*1,2,31*0,4*2,7*0,2,2*0,1,3*2,0,4*2,0,2,3*1,2,2*0,2,31*1,3*2,1 +,2*2,53*0,16*1,0,1,2,1,2*2,2*1,2,33*0,4*1,13*0,2,10*1,2,2*0,31*1,2 +,4*0,2,54*0,2,22*1,2,34*0,3*1,2,14*0,5*2,5*1,2,2*0,2,31*1,2,3*0 +,2*2,4*0,2,49*0,2,21*1,2,35*0,2*2,0,2,4*1,13*0,41*1,2,3*0,2,4*0,2 +,51*0,2,20*1,36*0,3*2,5*1,2,12*0,41*1,2,61*0,2,17*1,2,36*0,2,10*1 +,2,2*0,2,1,2*2,0,3*2,6*1,2*2,33*1,2,62*0,2*2,10*1,5*2,36*0,2,4*1/ c data ((lln2f1(i,j),i=1,180),j=31,45) +/16*1,0,2,5*1,2,0,2,34*1,2,63*0,2,6*1,2,7*0,2,35*0,2,20*1,2,0,2 +,5*1,2,2*0,2*2,30*1,2,65*0,2,5*1,42*0,2,23*1,2*0,2,5*1,8*0,2,24*1 +,2,67*0,2,4*1,42*0,24*1,2,0,2,6*1,2,1,2*2,4*0,3*2,7*1,2*2,10*1,2*2 +,69*0,2,3*1,41*0,2,25*1,2,0,2,7*1,2,7*0,2,6*1,3*0,2,6*1,2,73*0,2 +,3*1,41*0,26*1,2,2*0,7*1,2,8*0,4*1,2,5*0,2,4*1,2,75*0,2,2*1,2,2*0 +,2*2,36*0,27*1,2,0,2,4*1,2*2,9*0,3*1,2,6*0,2,0,4*1,2,6*0,2,69*0 +,4*2,1,2,36*0,27*1,2,0,2,2*1,2,12*0,2,2*1,9*0,2,4*1,80*0,4*2,34*0 +,28*1,2,0,2,14*0,2,1,2,10*0,2*2,2*1,2,81*0,2,1,34*0,2,28*1,2,2*0,2 +,13*0,2*2,12*0,2*2,83*0,2,5*0,3*2,27*0,31*1,2,13*0,2,103*0,3*2,3*1 +,2*2,23*0,2,6*1,3*2,21*1,2,15*0,2,100*0,2,8*1,2,23*0,2,5*1,3*0,2*2 +,18*1,2,26*0,2*2,6*0,2,82*0,2,9*1,3*2,21*0,2*2,8*0,17*1,2,28*0,2 +,4*0,2,1,2,82*0,2,12*1,2,30*0,16*1,2,27*0,2*2,4*0,2,2*1,2,81*0,2 +,13*1,2,25*0/ c data ((lln2f1(i,j),i=1,180),j=46,60) +/5*0,15*1,30*0,2*2,3*0,2,2*1,8*0,2,73*0,11*1,3*2,2*0,2,28*0,2,13*1 +,2,31*0,2*2,4*0,2,10*0,3*2,69*0,2,9*1,0,2*2,2*1,2,1,3*2,26*0,2 +,12*1,2,32*0,2,16*0,2,2*1,2,66*0,2,20*1,2*2,24*0,2,12*1,2,35*0,2 +,13*0,5*2,66*0,2,21*1,2,24*0,12*1,2,53*0,2,67*0,21*1,25*0,12*1,2 +,121*0,2,19*1,25*0,2,13*1,45*0,3*2,3*0,2,70*0,18*1,2,25*0,2,13*1 +,4*0,2,38*0,2,0,2,1,2,3*0,2*2,69*0,2*2,16*1,26*0,12*1,2*2,2*0,3*2 +,37*0,2,5*1,3*2,1,2,71*0,2,15*1,26*0,2,10*1,2,4*0,2,1,2,36*0,2 +,10*1,2,72*0,14*1,2,27*0,2,9*1,2,4*0,2,1,34*0,2*2,14*1,2,70*0,14*1 +,2,27*0,2,9*1,2,4*0,1,2,33*0,2,16*1,2,70*0,11*1,2*2,29*0,2,8*1,2 +,5*0,2*2,33*0,2,17*1,2,69*0,10*1,2,32*0,8*1,41*0,2,18*1,2,68*0 +,10*1,2,32*0,2,6*1,2,41*0,2,18*1,2,67*0,2,9*1,2,25*0/ c data ((lln2f1(i,j),i=1,180),j=61,75) +/9*0,2,4*1,2,43*0,6*1,3*2,9*1,68*0,2,8*1,2,35*0,2,1,3*2,44*0,1,3*2 +,5*0,3*2,5*1,2,68*0,2,8*1,2,95*0,2,5*1,69*0,2,5*1,2,99*0,2,1,3*2 +,68*0,2,7*1,116*0,2*2,54*0,2,4*1,2,174*0,2,3*1,0,2,104*0,2,12*0,2 +,57*0,3*1,117*0,2*2,58*0,2*1,2,176*0,2,2*1,2,176*0,2,1,2,178*0,1 +,755*0/ c data ((lln2f1(i,j),i=1,180),j=76,90) +/509*0,2,55*0,2,1,2,14*0,6*2,1,2,2*1,4*2,1,15*2,75*0,2,1,2,50*0 +,2*2,12*1,3*2,2*0,2,33*1,4*2,69*0,2,2*1,2,30*0,1,4*2,1,9*2,65*1 +,4*2,59*0,2,1,2*2,2*1,26*0,5*2,84*1,2,45*0,6*2,0,2*2,2*0,2,2*0,3*2 +,2*1,2,21*0,3*2,32*1,2,54*1,2,28*0,7*2,3*1,10*2,19*1,19*0,2*2,91*1 +,2,22*0,4*2,43*1,3*2,11*0,2*2,101*1,12*2,56*1,7*2,917*1/ c end EOR rm aqclib.o make aqclib.o