program trcdata c c --- program to interpolate daSilva COADS data (dimension idmh x jdmh) to a grid c --- of dimension idm x jdm c --- i index (latitude) increases N to S, j index (longitude) increases W to E c c --- this example is for 106 x 198 .5-deg. grid (DAMEE) c --- program interpolates 12 months of data for input field c c **** input parameters are from newtest.f program (also see function c **** 'scaly' below) parameter (idmh=56,jdmh=139) c **** output parameters are for example domain (DAMEE) parameter (idm=106,jdm=198) parameter (minum=1) c real cdata(idmh,jdmh),wkc(idmh,jdmh),cin(jdmh,idmh,2) real mask(idmh,jdmh) common/cclats/clats(idmh) character util(idm*jdm+14)*2,preambl(5)*79,preamb1(5)*79 character util1(idmh*jdmh+14)*2 integer numd(idmh,jdmh) real p(idmh,jdmh,2),hh(idmh,jdmh,2),arrayn(idm,jdm) dimension xyij(idm,jdm,2),coord(2),lrfi(6),lrfo(6) data coord/1.,2./,lrfi/6*0/,lrfo/6*0/ dimension depth(idm,jdm),ip(idm,jdm),plot1(jdmh,idmh), .plot2(jdm,idm),xlat(idm),work(idm,jdm) dimension mon(12) common/conre1/ioffp,spval data spval1/-.03125/ data ii,jj/idm,jdm/ data mon/'JAN.','FEB.','MAR.','APR.','MAY ','JUNE','JULY','AUG.', . 'SEP.','OCT.','NOV.','DEC.'/ c c **** change header for particular data file data hdr/'SST '/ c c **** change clong, clat according to what you chose for newtest.f c --- clong (clat) are W (N) longitude (latitude) limit of COADS data c --- cgridj, cgridi are resolution (deg) of COADS data in long, lat direction data clong/-100./,clat/55./,cgridj/1./,cgridi/1./ c c **** change rlong, xpiv, grid according to output grid (DAMEE here) c --- rlong is W longitude limit of output grid data rlong/-98./ c --- 'grid' is resolution in degrees of output Mercator grid c --- 'xpiv' is i location of equator data rad/57.29578/,pi/3.14159265/,grid/.5/,xpiv/117./ c --- formula relating latitude to equatorial distance on Mercator map alat(dist,gridd)=2.*atan(exp(dist*gridd/rad))-pi/2. c c --- **** read bottom topography for output domain open (unit=9,file='depthdamee.106x198',status='old') read (9,120) preambl,idim,jdim,length,(util(l),l=1,length) 120 format (5(a79/),3i6/(40a2)) close (unit=9) write (*,'(/(1x,a79))') preambl call unpakk(depth,idm,idim,jdim,util,length) do 602 i=1,idm do 602 j=1,jdm ip(i,j)=0 if (depth(i,j).gt.0.) ip(i,j)=1 602 continue c c --- clats contains COADS data latitudes do 1 i=1,idmh 1 clats(i)=clat-cgridi*float(i-1) c call opngks call gsclip(0) ioffp=1 spval=spval1 c --- calculate Mercator latitudes and x,y locations for interpolation do 4 i=1,idm xlat(i)=alat(xpiv-i,grid)*rad ri=(clat-xlat(i))/cgridi+1. do 4 j=1,jdm xlong=rlong+(j-1)*grid rj=(xlong-clong)/cgridj+1. xyij(i,j,1)=amin1(amax1(ri,1.),float(idmh)) 4 xyij(i,j,2)=amin1(amax1(rj,1.),float(jdmh)) c c **** mask file open (unit=20,file='mask.1deg',status='old') read (20,120) preamb1,idim,jdim,length,(util1(l),l=1,length) call unpakk(mask,idmh,idmh,jdmh,util1,length) c iuni=10 c **** 12-month COADS data file open (unit=10,file='sst.1deg',status='old') c --- loop over 12 months do 500 imon=1,12 c --- read in one month of COADS data print 901, imon 901 format(' imon=',i4) c --- store with i index from N to S, j index from W to E read (10,120) preamb1,idim,jdim,length,(util1(l),l=1,length) call unpakk(cdata,idmh,idmh,jdmh,util1,length) do 45 i=1,idmh do 45 j=1,jdmh if (mask(i,j).eq.0.) cdata(i,j)=spval1 45 continue c c --- insert a statement such as the following if plots are not desired for c --- every month if (imon.ne.1.and.imon.ne.7) go to 995 c --- plot COADS data call set(.05,.95,.05,.95,1.,float(jdmh),1.,float(idmh),1) call getset(xxa,xxb,yya,yyb,xa,xb,ya,yb,ltype) call mappos(xxa,xxb,yya,yyb) call maproj('CE',0.,0.,0.) call mapstc('OUTLINE','CO') call mapstl('DOT',.false.) xla1=clats(1) xlo1=clong xla2=clats(idmh) xlo2=clong+(jdmh-1)*cgridj call mapset('CORNERS',xla1,xlo1,xla2,xlo2) call mapsti('LA',0) call mapint call mapdrw call getset(xxa,xxb,yya,yyb,xa,xb,ya,yb,ltype) call set(xxa,xxb,yya,yyb,1.,float(jdmh),1.,float(idmh),1) do 132 i=1,jdmh do 132 j=1,idmh 132 plot1(i,j)=cdata(idmh-j+1,i) call conrec(plot1,jdmh,jdmh,idmh,0.,0.,0.,1,-1,-'125252'o) call pwrit(5.0,float(idmh)+2.,mon(imon),4,2,0,-1) call wtstr(float(jdmh/5),float(idmh)+2.,' - COADS DATA', & 2,0,-1) call pwrit(float(jdmh/5),float(idmh)+2.,hdr,4,2,0,-1) call perim(1,jdmh-1,1,idmh-1) call frame 995 continue c c --- 'embroider' HR data to fill points needed for interpolated grid do 44 i=1,idmh do 44 j=1,jdmh numd(i,j)=1 if (mask(i,j).eq.0.) then numd(i,j)=-1 cdata(i,j)=spval1 endif 44 continue do 824 i=1,idm do 824 j=1,jdm if (ip(i,j).eq.0) go to 824 ih1=xyij(i,j,1) jh1=xyij(i,j,2) do 825 ih=max0(1,ih1-1),min0(idmh,ih1+2) do 825 jh=max0(1,jh1-1),min0(jdmh,jh1+2) if (numd(ih,jh).lt.0) numd(ih,jh)=0 825 continue 824 continue call amend(cdata,numd,idmh,idmh,jdmh,minum,wkc) c --- interpolate data to Mercator grid c --- note 'grdtrns' requires at least 2 levels of data in the vertical. c --- here we make the 2 levels contain the same data, and set the values c --- in array 'p' to 0. and 1., and then interpolate to the value 0.5 do 2 i=1,idmh do 2 j=1,jdmh hh(i,j,1)=cdata(i,j) hh(i,j,2)=cdata(i,j) p(i,j,1)=0. 2 p(i,j,2)=1. c do 910 i=1,idm do 910 j=1,jdm 910 arrayn(i,j)=spval1 call grdtrns(0,p,hh,dumm,idmh,jdmh,2,coord,lrfi,arrayn,arrayn, .idm,jdm,1,0.5,lrfo,xyij) do 922 i=1,idm do 922 j=1,jdm 922 if (ip(i,j).eq.0) arrayn(i,j)=spval1 c call filtr1(arrayn,work,idm,idm,jdm) c c --- insert a statement such as the following if plots are not desired for c --- every month if (imon.ne.1.and.imon.ne.7) go to 996 call getset(xxa,xxb,yya,yyb,xa,xb,ya,yb,ltype) call mappos(xxa,xxb,yya,yyb) call maproj('MERCATOR',0.,0.,0.) call mapstc('OUTLINE','CO') call mapstl('DOT',.false.) xla1=xlat(1) xlo1=rlong xla2=xlat(idm) xlo2=rlong+(jdm-1)*grid call mapset('CORNERS',xla1,xlo1,xla2,xlo2) call mapsti('LA',0) call mapint call mapdrw call getset(xxa,xxb,yya,yyb,xa,xb,ya,yb,ltype) call set(xxa,xxb,yya,yyb,1.,float(jdm),1.,float(idm),1) do 131 i=1,jdm do 131 j=1,idm 131 plot2(i,j)=arrayn(idm-j+1,i) spval=spval1 call conrec(plot2,jdm,jdm,idm,0.,0.,0.,1,-1,-'125252'o) do 822 i=1,idm do 822 j=1,jdm if (ip(i,j).eq.0) arrayn(i,j)=spval1 822 continue call pwrit(5.0,float(idm)+3.,mon(imon),4,2,0,-1) call pwrit(float(jdm/8),float(idm)+3., & 26h - NORTH ATLANTIC (.5) ,26,2,0,-1) call pwrit(float(jdm/8),float(idm)+3.,hdr,4,2,0,-1) call perim(1,jdm-1,1,idm-1) call frame 996 continue c c --- write this month's data set. c c **** output (12 months on Mercator grid for DAMEE) will be on unit 20 iuno=20 call pakk(arrayn,idm,idm,jdm,util,length) write (iuno,120) preambl,idm,jdm,length,(util(l),l=1,length) 100 format (3i6/(40a2)) c 510 continue 500 continue call clsgks end subroutine amend(data,nreps,idm,ii,jj,minum,work) c c --- given an array incompletely filled with grid point values representing c --- grid-square averages of raw observations, use an objective analysis c --- approach to (a) fill in missing grid point values and (b) spatially c --- average grid point values that are based on too few observations. c c --- i n p u t variables: c c --- data : the array of grid point values to be amended and smoothed. c --- nreps : the number of raw observations individual grid point values are c --- based on; a zero indicates that the corresponding value in 'data' c --- is missing; points to be ignored are identified by neg. values. c --- idm : first dimension of 'data' and 'nreps' in the calling program. c --- ii,jj : the number of grid points in i,j direction. c --- minum : the minimum number of observations needed to obtain a 'reliable' c --- grid point average. c --- work : work space, same dimensions as 'data' and 'nreps'. c c --- o u t p u t : the amended and smoothed field is returned in 'data' c parameter (nsrch=20) c --- nsrch = radius of influence (in grid units) searched for data to be used c --- in averaging process. dimension data(idm,1),nreps(idm,1),work(idm,1), . dist(nsrch+1,nsrch+1) data huge/1.e33/ data itest,jtest/7,23/ c c --- objective analysis weight function wfunc(d)=exp(-d) c c --- map scale factors in i,j direction for lateral distance calculation scalx(i,j)=1. ccc scaly(i,j)=1. c do 8 j=1,jj do 8 i=1,ii 8 work(i,j)=data(i,j) c c --- main analysis loop do 1 j=1,jj do 1 i=1,ii c c --- don't modify data points based on more than -minum- observations if (nreps(i,j).lt.0.or.nreps(i,j).ge.minum) go to 1 c c --- sort neighboring points by distance scx=scalx(i,j) scy=scaly(i,j) c do 2 j1=0,nsrch do 2 i1=0,nsrch 2 dist(i1+1,j1+1)=(scx*i1)**2+(scy*j1)**2 c nrp=0 value=0. sum=0. c do 4 m=1,(nsrch+1)**2 c dmin=huge do 3 j1=0,nsrch do 3 i1=0,nsrch if (dmin.gt.dist(i1+1,j1+1)) then imin=i1 jmin=j1 dmin=dist(i1+1,j1+1) end if 3 continue c c --- (imin,jmin) is the closest point not yet used dist(imin+1,jmin+1)=huge c c --- each (imin,jmin) combination yields up to 4 data points. n4=0 do 5 j2=j-jmin,j+jmin,max(1,2*jmin) do 5 i2=i-imin,i+imin,max(1,2*imin) if (i2.lt.1.or.j2.lt.1.or.i2.gt.ii.or.j2.gt.jj) go to 5 if (nreps(i2,j2).gt.0) n4=n4+nreps(i2,j2) 5 continue if (n4.le.0) go to 4 wgt=wfunc(dmin/(9.*(scx*scx+scy*scy))) c do 6 j2=j-jmin,j+jmin,max(1,2*jmin) do 6 i2=i-imin,i+imin,max(1,2*imin) if (i2.lt.1.or.j2.lt.1.or.i2.gt.ii.or.j2.gt.jj) go to 6 n=nreps(i2,j2) if (n.le.0) go to 6 c c --- compute weighted average of original grid point values c --- weights depend on distance and number of reports at each grid point ccc weight=wgt*sqrt(float(n)) weight=wgt*float(n) if (nrp+n4.gt.minum) weight=weight*float(minum-nrp)/float(n4) sum=sum+weight ccc if (iabs(i-itest)+iabs(j-jtest).le.1) write (*,'('' i,j='',2i3, ccc . '' i2,j2='',2i3,'' nreps,wgt,value='',i5,f9.3,1pe13.3)') ccc . i,j,i2,j2,nreps(i2,j2),weight,work(i2,j2) value=value+work(i2,j2)*weight 6 continue c c --- stop once enough reports have been found nrp=nrp+n4 if (nrp.ge.minum) go to 7 4 continue write (*,'('' search around i,j='',2i3,'' yields only'',i3, . '' reports. increase search radius.'')') i,j,nrp cc if (nrp.le.0) stop 'error in amend' c 7 if (sum.ne.0.) data(i,j)=value/sum ccc if (iabs(i-itest)+iabs(j-jtest).le.1) write (*,'('' i,j='',2i3, ccc . '' final value:''f9.3)') i,j,data(i,j) c 1 continue c return end function scaly(i,j) c **** change this to be same as main program parameter (idmh=56) common/cclats/clats(idmh) scaly=cos(clats(i)/57.29578) return end subroutine filtr1(a,b,idm,ii,jj) c c --- a: array to be smoothed c --- b: work space c real a(idm,1),b(idm,1) data wgt1/.5/,wgt2/.25/ data flag/-.03125/ cvmgp(aa,bb,cc)=aa*(.5+sign(.5,cc))+bb*(.5-sign(.5,cc)) cvmgz(aa,bb,cc)=cvmgp(aa,bb,-1.*abs(cc)) c do 3 j=1,jj ja=max0( 1,j-1) jb=min0(jj,j+1) do 3 i=1,ii qja=cvmgz(a(i,j),a(i,ja),a(i,ja)-flag) qjb=cvmgz(a(i,j),a(i,jb),a(i,jb)-flag) 3 b(i,j)=cvmgz(a(i,j),wgt1*a(i,j)+wgt2*(qja+qjb),a(i,j)-flag) c do 4 i=1,ii ia=max0( 1,i-1) ib=min0(ii,i+1) do 4 j=1,jj qia=cvmgz(b(i,j),b(ia,j),b(ia,j)-flag) qib=cvmgz(b(i,j),b(ib,j),b(ib,j)-flag) 4 a(i,j)=cvmgz(b(i,j),wgt1*b(i,j)+wgt2*(qia+qib),b(i,j)-flag) c return end subroutine grdtrns (modev, - qi, fi, dfi, idin, jdin, kdin, pki, lrfi, - po, fo, idou, jdou, kdou, qko, lrfo, xyij) c c c dimension qi(idin,jdin,kdin), fi(idin,jdin,kdin), - dfi(idin,jdin,kdin), - pki(kdin) dimension po(idou,jdou,kdou), fo(idou,jdou,kdou), - qko(kdou) dimension xyij(1) dimension lrfi(1), lrfo(1) dimension fab(2), qab(2), dfab(2) equivalence (fab(1),fija), (fab(2),fijb), - (qab(1),qija), (qab(2),qijb), - (dfab(1),dfa), (dfab(2),dfb) data epsc /1.0e-8/ common /hmtcm1/ xij,nxi,dx,ipth, - yij,myj,dy,jpth, - ipl,ipu c c c statement functions c iflor(x) = int(x) - int(.5-sign(.5,x)) hmtint (px,dx,dfx1,fx1,fx2,dfx2) = fx1 + px*(px*(2.*px-3.)* - (fx1-fx2) + (px-1.)*dx*(px*(dfx1+dfx2) - dfx1) ) c c c c ===================================================================== c ===================================================================== c == == c == == c == subroutine grdtrns == c == == c == == c ===================================================================== c ===================================================================== c c c *** purpose *** c c given volumes of a function f and new vertical coordinate p c on an old coordinate grid (xi,yj,qk), to transform to volumes c of f and the old vertical coordinate q on a new coordinate c grid (xii,yjj,pkk). c c c *** method *** c c transformation to new horizontal grid points on old vertical c coordinate surfaces is via piecewise bicubic quasi-hermite c interpolation -- quasi in the sense that the derivatives d*/dx c and d*/dy are obtained by second order centered finite differ- c ences. vertical transform is then performed to new vertical c coordinate surfaces at the new horizontal grid point using c either linear or hermite cubic interpolation. c c c *** argument list *** c c on input c c qi - a three dimensional volume of new vertical coordinate c data on the old coordinate grid. c fi - a three dimensional volume of function values on the old c coordinate grid. c dfi - a volume of function value vertical derivatives (i.e., of c the fi function), optional, need only be supplied if c hermite interpolation is to be used in the vertical. c idin - first dimension of qi, fi, dfi in the calling program. c jdin - second dimension of qi, fi, dfi in calling program. c kdin - number of vertical levels in the volumes qi, fi, dfi c for this call. must be at least 2. c pki - vector of old vertical coordinate values associated c with vertical levels 1,...,kdin in input volumes. c must be monotone (increasing or decreasing). c modev - vertical interpolation mode flag, c = -1, hermite, and dfi contains df/dq, the derivative c with respect to the new vertical coordinate. c = 0, linear. c = +1, hermite, and dfi contains df/dp, the derivative c with respect to the old vertical coordinate. c lrfi - we define the 'total old coordinate space' to be the c full three dimensional grid in old coordinate space c for which the user has new coord and function data. c on this call to this routine, the data in the input c arrays may entirely span the horizontal dimensions of c the total old coord space, or alternately, the data c may only span a subspace in one or both of the hori- c zontal dimensions. in the latter case, it is assumed c that this call is one of a sequence of calls which c will cycle through the total input space and thus c span the deficient dimension(s). in either case, the c actual data for this call may not fill the user c declared dimensions of qi, fi, dfi -- idin and jdin. c lrfi(1),lrfi(2) - the data for this call, for level *, is in c array(1,1,*) to array(lrfi(1),lrfi(2),*) . c if .le.0, idin and/or jdin taken as defaults. other- c wise must be .ge.4, or an error is diagnosed. c lrfi(3),lrfi(4) - the horizontal first and second dimensions c of the users total old coordinate space. if .le.0, c assumed to be same as lrfi(1),lrfi(2) (or their c defaults idin,jdin ). c lrfi(5),lrfi(6) - horizontal grid point (1,1,*) of the c input volumes corresponds to horizontal grid point c (lrfi(5),lrfi(6),*) in the users total old coord c space. if .le.0, default value(s) 1 are used. c ** note ** the information in lrfi allows us to completely c locate the horizontal grid of the actual data for this call c as a subgrid of the total horizontal grid from which the c user is working. from this, we can tell when an apparent c boundary grid square is actually a boundary square in the c total old coord space, and interpolate with a boundary c formula, or not interpolate, as appropriate. c idou - first dimension of arrays po, fo in calling program. c jdou - second dimension of arrays po, fo in calling program. c kdou - number of vertical levels of new coordinate space to c construct on this call. c qko - vector of new vertical coordinate values of the kdou c output levels desired on this call. must be monotone. c lrfo - similar to lrfi, except applies to that horizontal c portion of the total new coordinate space (output c grid) which is to be constructed on this call. c lrfo(1),lrfo(2) - on this call, construct a lrfo(1) by c lrfo(2) subgrid of total new coord space and store in c array(1,1,*) to array(lrfo(1),lrfo(2),*) of output c arrays. default values of idou and/or jdou are used c if .le.0. c lrfo(3),lrfo(4) - horizontal first and second dimensions of c users total new coord space. if .le.0, default values c lrfo(1) and/or lrfo(2) (or their defaults) are used. c lrfo(5),lrfo(6) - horizontal origin of this-call-subgrid c of total new coord space in horizontal grid of new c coord space. default values 1 and/or 1 if .le.0 . c xyij - array containing old horizontal grid point coords of c new horizontal grid points. the horizontal first and c second dimensions should always exactly match those of c the total new coord space (lrfo(3) and lrfo(4) or c their defaults), or else the array should be loaded in c such a way to mimic this dimensionality. thus for c io=1,..,lrfo(3) and jo=1,..,lrfo(4) we should have c xyij(io,jo,1) = i old grid coord of point c (io,jo,*) in total new grid c xyij(io,jo,2) = j old grid coord of point c (io,jo,*) in total new grid. c on this call we will be examining c xyij( lrfo(5), lrfo(6), *) c to c xyij( lrfo(5)+lrfo(1)-1, lrfo(6)+lrfo(2)-1, *) . c **note** the user may write a 'real function xyij(idx)', c declare it external in the calling program, and delete c the dimension statement above. the single argument c idx will be the linear subscript equivalent of the c three dimensional subscript (io,jo,1) or (io,jo,2). c c on output c c po - three dimensional volume of old vertical coordinate data c p transformed to new coordinate space. c fo - three dimensional volume of function values on the new c coordinate grid. if any grid point of the portion of c the total new coordinate grid to be constructed this c call is exterior to the portion of the old coord grid c supplied for this call, then the corresponding entries c of the output arrays po and fo will be left undefined. c c c ===================================================================== c ===================================================================== c c c c ************************** c * * c **** initialization **** c * * c ************************** c c check input dimensions and pointers, set defaults c kdi = kdin idi = lrfi(1) if (idi.le.0) idi = idin jdi = lrfi(2) if (jdi.le.0) jdi = jdin iit = lrfi(3) if (iit.le.0) iit = idi jit = lrfi(4) if (jit.le.0) jit = jdi iil = max0 (1,lrfi(5)) jil = max0 (1,lrfi(6)) iiu = iil + idi - 1 jiu = jil + jdi - 1 if ( max0(4,min0(idi,idin)) .ne. idi ) go to 9001 if ( max0(4,min0(jdi,jdin)) .ne. jdi ) go to 9002 c c check output dimensions and pointers, set defaults c kdo = kdou ido = lrfo(1) if (ido.le.0) ido = idou jdo = lrfo(2) if (jdo.le.0) jdo = jdou iot = lrfo(3) if (iot.le.0) iot = ido jot = lrfo(4) if (jot.le.0) jot = jdo iol = max0 (1,lrfo(5)) jol = max0 (1,lrfo(6)) iou = iol + ido - 1 jou = jol + jdo - 1 if (ido.gt.idou .or. jdo.gt.jdou) go to 9011 if (iou.gt.iot .or. jou.gt.jot) go to 9012 c c set boundary interpolation and vertical mode flags c ibdl = min0 (iil,2) jbdl = min0 (jil,2) ibdu = min0 (iit-iiu+1, 2) jbdu = min0 (jit-jiu+1, 2) jphmt = modev c c preset other constants c idq = idin ijot = iot*jot idm2 = idi - 2 jdm2 = jdi - 2 idm1 = idi - 1 jdm1 = jdi - 1 xnrm = float(iil) - 1.0 ynrm = float(jil) - 1.0 xmx = float(idi) ymx = float(jdi) c c c ***************************** c * * c *** grid transformation *** c * * c ***************************** c c c ===== output grid horizontal loop ===== c ios = 0 do 7001 ixo = iol,iou ios = ios + 1 jos = 0 do 7000 jyo = jol,jou jos = jos + 1 idx = ixo + (jyo-1)*iot jdy = idx + ijot xi = xyij(idx) - xnrm yj = xyij(jdy) - ynrm if (xi.gt.xmx .or. yj.gt.ymx) go to 7000 nxi = min0 (idm1,iflor(xi)) myj = min0 (jdm1,iflor(yj)) c c compute processing path selectors c ipth = max0 (ibdu*(nxi-idm2), ibdl*min0(0,nxi-2) ) jpth = max0 (jbdu*(myj-jdm2), jbdl*min0(0,myj-2) ) if ( max0(iabs(ipth),iabs(jpth)) .ge. 2 ) go to 7000 dx = xi - float (nxi) dy = yj - float(myj) ipl = max0 (nxi-1, 1) ipu = min0 (nxi+2, idi) c c ===== input vertical levels loop ===== c last = -1 lvi = 0 1900 lvi = lvi + 1 if (lvi.gt.kdi) go to 7000 qijb = qija fijb = fija dfb = dfa qija = qsihmt (qi(1,1,lvi),idq) if (lvi.eq.1) go to 1900 c c ===== output vertical levels loop ===== c lvo = 0 2000 lvo = lvo + 1 if (lvo.gt.kdo) go to 1900 if ( (qko(lvo)-qija) * (qko(lvo)-qijb) .gt. 0.0) go to 2000 dq = qija - qijb if (dq.ne.0.) pq = (qko(lvo)-qijb)/dq ompq = 1.0 - pq lvbrn = max0 (-1, (last+1) - lvi ) if (pq*ompq .gt. epsc) go to 2015 idxc = ifix(lvi - .5 + pq) po(ios,jos,lvo) = pki(idxc) lf1 = lvi - idxc + 1 lf2 = lf1 + lvbrn + 1 if (lf2.ge.3) go to 2004 fab(lf1) = qsihmt (fi(1,1,idxc),idq) 2004 fo(ios,jos,lvo) = fab(lf1) go to (2000,2006,2000,2000), lf2 2006 last = idxc if (jphmt) 2008,2000,2008 2008 dfab(lf1) = qsihmt(dfi(1,1,idxc),idq) go to 2000 2015 pok = pq*pki(lvi) + ompq*pki(lvi-1) last = lvi po(ios,jos,lvo) = pok if (lvbrn) 2025,2030,2035 2025 fijb = qsihmt (fi(1,1,lvi-1),idq) 2030 fija = qsihmt (fi(1,1,lvi),idq) 2035 if (jphmt) 2050,2040,2060 2040 fo(ios,jos,lvo) = pq*fija + ompq*fijb go to 2000 2050 dc = dq go to 2070 2060 dc = pki(lvi) - pki(lvi-1) 2070 pc = pq if (lvbrn) 2075,2080,2085 2075 dfb = qsihmt (dfi(1,1,lvi-1),idq) 2080 dfa = qsihmt (dfi(1,1,lvi),idq) 2085 fo(ios,jos,lvo) = hmtint (pc,dc,dfb,fijb,fija,dfa) go to 2000 7000 continue 7001 continue c c successful completion, normal return c 7007 return c c error conditions c 9001 stop 99001 9002 stop 99002 9011 stop 99011 9012 stop 99012 c end real function qsihmt(fcn,idf) dimension fcn(idf,1) real rw(4) common /hmtcm1/ xij,nxi,dx,ipth, - yij,myj,dy,jpth, - ipl,ipu hmtint (dl,f1,f2,df1,df2) = - f1 + dl*( dl*(2.*dl-3.)*(f1-f2) + - (dl-1.)*(dl*(df1+df2)-df1) ) c iw = 0 do 1050 ir = ipl,ipu iw = iw + 1 fcnl1=fcn(ir,myj-1) fcnl = fcn(ir,myj) fcnu = fcn(ir,myj+1) fcnu1=fcn(ir,myj+2) if (jpth) 1020,1015,1010 1010 drvu = .5*fcn(ir,myj-1) - 2.*fcn(ir,myj) + 1.5*fcn(ir,myj+1) drvl = .5*(fcn(ir,myj+1) - fcn(ir,myj-1)) go to 1025 1015 drvl = .5*(fcn(ir,myj+1) - fcn(ir,myj-1)) drvu = .5*(fcn(ir,myj+2) - fcn(ir,myj)) go to 1025 1020 drvl = -1.5*fcn(ir,myj) + 2.*fcn(ir,myj+1) - .5*fcn(ir,myj+2) drvu = .5*(fcn(ir,myj+2) - fcn(ir,myj)) 1025 continue rw(iw) = hmtint(dy,fcnl,fcnu,drvl,drvu) 1050 continue c jyw = 2 + min0(ipth,0) fcnl = rw(jyw) fcnu = rw(jyw+1) if (ipth) 2020,2015,2010 2010 drvu = .5*rw(jyw-1) - 2.*rw(jyw) + 1.5*rw(jyw+1) drvl = .5*(rw(jyw+1) - rw(jyw-1)) go to 2025 2015 drvl = .5*(rw(jyw+1) - rw(jyw-1)) drvu = .5*(rw(jyw+2) - rw(jyw)) go to 2025 2020 drvl = -1.5*rw(jyw) + 2.*rw(jyw+1) - .5*rw(jyw+2) drvu = .5*(rw(jyw+2) - rw(jyw)) 2025 continue c qsihmt = hmtint(dx,fcnl,fcnu,drvl,drvu) c return end subroutine pakk(array,idim,ii,jj,compac,length) c c --- converts the contents of -array- into an ascii character string which c --- is stored in character*2 array -compac-. compac(1)...compac(7) contain c --- the base value, i.e., the minimum value encountered in -array-. c --- compac(8)...compac(14) contain a scale factor by which the individual c --- 6-bit integers encoded as ascii character pairs in compac(8),... c --- compac(length) must be multiplied before the base value is added c --- during an unpakking operation. base value and scale factor are encoded c --- in e14.7 format. c c --- the printable ascii characters used to encode the integers include c --- the numbers 0...9, upper- and lower-case letters a...z, a...z, plus c --- two additional characters '.' and '/' (total of 64). c c --- a packing operation fills (ii*jj+14) array elements in -compac- which c --- must be dimensioned accordingly in the calling program. the total c --- number of occupied array elements is returned in -length-. in calls to c --- unpack, -length- is treated as input variable. c real array(idim,1) character*2 char,compac(1),comp2(14) character*14 comp14(2) equivalence (comp2,comp14) data nbits/12/ base=1.e22 do 1 i=1,ii do 1 j=1,jj 1 base=amin1(base,array(i,j)) scal=0. do 2 i=1,ii do 2 j=1,jj 2 scal=amax1(scal,array(i,j)-base) scal=scal/float(2**nbits-1) i1=0 i2=0 length=14 do 3 i=1,ii do 3 j=1,jj if (scal.eq.0.) go to 7 numb=(array(i,j)-base)/scal+.5 i1=numb/64 i2=numb-64*i1 c c --- map 6-bit numbers onto character set consisting of numbers c --- 0...9, letters a...z, a...z, and the two characters '.' and '/'. c --- (if mapping into the character range 32...95 -- which includes the c --- characters !"#$%&'()*+,-./:;<=>?@[\]^_ -- is deemed safe, delete c --- the next 6 lines.) ccc if (i1.gt.37) i1=i1+6 ccc if (i1.gt.11) i1=i1+7 ccc i1=i1+14 ccc if (i2.gt.37) i2=i2+6 ccc if (i2.gt.11) i2=i2+7 ccc i2=i2+14 c 7 length=length+1 compac(length)(1:1)=char(i1+32) compac(length)(2:2)=char(i2+32) 100 format (a2) 3 continue write (comp14(1),101) base write (comp14(2),101) scal 101 format (1pe14.7) do 8 l=1,14 8 compac(l)=comp2(l) c return c c entry unpakk(array,idim,ii,jj,compac,length) c do 9 l=1,14 9 comp2(l)=compac(l) read (comp14(1),101) base read (comp14(2),101) scal lngth=14 do 4 i=1,ii do 4 j=1,jj lngth=lngth+1 i1=ichar(compac(lngth)(1:1)) i2=ichar(compac(lngth)(2:2)) c c --- 6-bit numbers are mapped onto character set consisting of numbers c --- 0...9, letters a...z, a...z, and the two characters '.' and '/'. c --- (if mapped into character range 32...95, delete next 6 lines) ccc if (i1.gt.96) i1=i1-6 ccc if (i1.gt.64) i1=i1-7 ccc i1=i1-14 ccc if (i2.gt.96) i2=i2-6 ccc if (i2.gt.64) i2=i2-7 ccc i2=i2-14 c 4 array(i,j)=scal*float(64*(i1-32)+(i2-32))+base if (lngth.ne.length) stop 'unpack' return end