C subroutine tides(time,x,mx,zp,mz,tnms,shtms,wout) dimension wout(5,mx,mz),x(mx),zp(mz),tnms(mz),shtms(mz) dimension hs(mz),us(mz),vs(mz),ts(mz),ws(mz),xint(mz) c Warning: this subroutine assumes dx= HHHH/20 km (0 to HHHH km; c mx=21) and dz=.2 (-14.0 to 1.0; mz=76). If these values are changed c externally, they must be changed here also. c u=positive eastward (m/s), v=positive northward (m/s) c w=positive upward (s-1), zp=pressure scale height (-14.0 to 1.0) HHHH=8000. pi=acos(-1.) hwave=2.*pi/HHHH c Get height vectors of u,v,T for heights zp(i) at t=local time call tide18v(zp,mz,time,us,vs,ts) C Calculate perturbation height do i=1,mz tt=ts(i)/tnms(i) hs(i)=tt/(1.+tt) hs(i)=-shtms(i)*hs(i) end do C Integrate up (Simpson's rule) call qsf(.2,hs,hs,mz) C Calculate vertical velocity from continuity equation k=mz do i=1,mz xint(i)=vs(k)*exp(-zp(k)) k=k-1 end do C Integrate down (negative dz accounted for in eq.) call qsf(.2,xint,xint,mz) k=mz do i=1,mz ws(k)=(hwave/1000.)*(xint(i)*exp(zp(k))+vs(mz)*exp(zp(k)-zp(mz))) k=k-1 end do do 3 i=1,mz do 4 j=1,mx xx=x(j)/1.0e+05 flat=-cos(xx*hwave) slat=sin(xx*hwave) wout(1,j,i)=us(i)*flat wout(2,j,i)=vs(i)*flat wout(3,j,i)=ws(i)*slat wout(4,j,i)=ts(i)*flat wout(5,j,i)=hs(i)*flat 4 continue 3 continue return end subroutine terms (diag,sdiag,sigma,del) c real diag,sdiag,sigma,del c c coded by alan kaylor cline c from fitpack -- january 26, 1987 c a curve and surface fitting package c a product of pleasant valley software c 8603 altus cove, austin, texas 78759, usa c c this subroutine computes the diagonal and superdiagonal c terms of the tridiagonal linear system associated with c spline under tension interpolation. c c on input-- c c sigma contains the tension factor. c c and c c del contains the step size. c c on output-- c c sigma*del*cosh(sigma*del) - sinh(sigma*del) c diag = del*--------------------------------------------. c (sigma*del)**2 * sinh(sigma*del) c c sinh(sigma*del) - sigma*del c sdiag = del*----------------------------------. c (sigma*del)**2 * sinh(sigma*del) c c and c c sigma and del are unaltered. c c this subroutine references package module snhcsh. c c----------------------------------------------------------- c if (sigma .ne. 0.) go to 1 diag = del/3. sdiag = del/6. return 1 sigdel = sigma*del call snhcsh (sinhm,coshm,sigdel,0) denom = sigma*sigdel*(1.+sinhm) diag = (coshm-sinhm)/denom sdiag = sinhm/denom return end subroutine curv1 (n,x,y,slp1,slpn,islpsw,yp,temp, * sigma,ierr) c integer n,islpsw,ierr real x(n),y(n),slp1,slpn,yp(n),temp(n),sigma c c coded by alan kaylor cline c from fitpack -- january 26, 1987 c a curve and surface fitting package c a product of pleasant valley software c 8603 altus cove, austin, texas 78759, usa c c this subroutine determines the parameters necessary to c compute an interpolatory spline under tension through c a sequence of functional values. the slopes at the two c ends of the curve may be specified or omitted. for actual c computation of points on the curve it is necessary to call c the function curv2. c c on input-- c c n is the number of values to be interpolated (n.ge.2). c c x is an array of the n increasing abscissae of the c functional values. c c y is an array of the n ordinates of the values, (i. e. c y(k) is the functional value corresponding to x(k) ). c c slp1 and slpn contain the desired values for the first c derivative of the curve at x(1) and x(n), respectively. c the user may omit values for either or both of these c parameters and signal this with islpsw. c c islpsw contains a switch indicating which slope data c should be used and which should be estimated by this c subroutine, c = 0 if slp1 and slpn are to be used, c = 1 if slp1 is to be used but not slpn, c = 2 if slpn is to be used but not slp1, c = 3 if both slp1 and slpn are to be estimated c internally. c c yp is an array of length at least n. c c temp is an array of length at least n which is used for c scratch storage. c c and c c sigma contains the tension factor. this value indicates c the curviness desired. if abs(sigma) is nearly zero c (e.g. .001) the resulting curve is approximately a c cubic spline. if abs(sigma) is large (e.g. 50.) the c resulting curve is nearly a polygonal line. if sigma c equals zero a cubic spline results. a standard value c for sigma is approximately 1. in absolute value. c c on output-- c c yp contains the values of the second derivative of the c curve at the given nodes. c c ierr contains an error flag, c = 0 for normal return, c = 1 if n is less than 2, c = 2 if x-values are not strictly increasing. c c and c c n, x, y, slp1, slpn, islpsw and sigma are unaltered. c c this subroutine references package modules ceez, terms, c and snhcsh. c c----------------------------------------------------------- c nm1 = n-1 np1 = n+1 ierr = 0 if (n .le. 1) go to 8 if (x(n) .le. x(1)) go to 9 c c denormalize tension factor c sigmap = abs(sigma)*float(n-1)/(x(n)-x(1)) c c approximate end slopes c if (islpsw .ge. 2) go to 1 slpp1 = slp1 go to 2 1 delx1 = x(2)-x(1) delx2 = delx1+delx1 if (n .gt. 2) delx2 = x(3)-x(1) if (delx1 .le. 0. .or. delx2 .le. delx1) go to 9 call ceez (delx1,delx2,sigmap,c1,c2,c3,n) slpp1 = c1*y(1)+c2*y(2) if (n .gt. 2) slpp1 = slpp1+c3*y(3) 2 if (islpsw .eq. 1 .or. islpsw .eq. 3) go to 3 slppn = slpn go to 4 3 delxn = x(n)-x(nm1) delxnm = delxn+delxn if (n .gt. 2) delxnm = x(n)-x(n-2) if (delxn .le. 0. .or. delxnm .le. delxn) go to 9 call ceez (-delxn,-delxnm,sigmap,c1,c2,c3,n) slppn = c1*y(n)+c2*y(nm1) if (n .gt. 2) slppn = slppn+c3*y(n-2) c c set up right hand side and tridiagonal system for yp and c perform forward elimination c 4 delx1 = x(2)-x(1) if (delx1 .le. 0.) go to 9 dx1 = (y(2)-y(1))/delx1 call terms (diag1,sdiag1,sigmap,delx1) yp(1) = (dx1-slpp1)/diag1 temp(1) = sdiag1/diag1 if (n .eq. 2) go to 6 do 5 i = 2,nm1 delx2 = x(i+1)-x(i) if (delx2 .le. 0.) go to 9 dx2 = (y(i+1)-y(i))/delx2 call terms (diag2,sdiag2,sigmap,delx2) diag = diag1+diag2-sdiag1*temp(i-1) yp(i) = (dx2-dx1-sdiag1*yp(i-1))/diag temp(i) = sdiag2/diag dx1 = dx2 diag1 = diag2 5 sdiag1 = sdiag2 6 diag = diag1-sdiag1*temp(nm1) yp(n) = (slppn-dx1-sdiag1*yp(nm1))/diag c c perform back substitution c do 7 i = 2,n ibak = np1-i 7 yp(ibak) = yp(ibak)-temp(ibak)*yp(ibak+1) return c c too few points c 8 ierr = 1 return c c x-values not strictly increasing c 9 ierr = 2 return end function curv2 (t,n,x,y,yp,sigma) c integer n real t,x(n),y(n),yp(n),sigma c c coded by alan kaylor cline c from fitpack -- january 26, 1987 c a curve and surface fitting package c a product of pleasant valley software c 8603 altus cove, austin, texas 78759, usa c c this function interpolates a curve at a given point c using a spline under tension. the subroutine curv1 should c be called earlier to determine certain necessary c parameters. c c on input-- c c t contains a real value to be mapped onto the interpo- c lating curve. c c n contains the number of points which were specified to c determine the curve. c c x and y are arrays containing the abscissae and c ordinates, respectively, of the specified points. c c yp is an array of second derivative values of the curve c at the nodes. c c and c c sigma contains the tension factor (its sign is ignored). c c the parameters n, x, y, yp, and sigma should be input c unaltered from the output of curv1. c c on output-- c c curv2 contains the interpolated value. c c none of the input parameters are altered. c c this function references package modules intrvl and c snhcsh. c c----------------------------------------------------------- c c determine interval c im1 = intrvl(t,x,n) i = im1+1 c c denormalize tension factor c sigmap = abs(sigma)*float(n-1)/(x(n)-x(1)) c c set up and perform interpolation c del1 = t-x(im1) del2 = x(i)-t dels = x(i)-x(im1) sum = (y(i)*del1+y(im1)*del2)/dels if (sigmap .ne. 0.) go to 1 curv2 = sum-del1*del2*(yp(i)*(del1+dels)+yp(im1)* * (del2+dels))/(6.*dels) return 1 sigdel = sigmap*dels call snhcsh (ss,dummy,sigdel,-1) call snhcsh (s1,dummy,sigmap*del1,-1) call snhcsh (s2,dummy,sigmap*del2,-1) curv2 = sum+(yp(i)*del1*(s1-ss)+yp(im1)*del2*(s2-ss))/ * (sigdel*sigmap*(1.+ss)) return end subroutine ceez (del1,del2,sigma,c1,c2,c3,n) c real del1,del2,sigma,c1,c2,c3 c c coded by alan kaylor cline c from fitpack -- january 26, 1987 c a curve and surface fitting package c a product of pleasant valley software c 8603 altus cove, austin, texas 78759, usa c c this subroutine determines the coefficients c1, c2, and c3 c used to determine endpoint slopes. specifically, if c function values y1, y2, and y3 are given at points x1, x2, c and x3, respectively, the quantity c1*y1 + c2*y2 + c3*y3 c is the value of the derivative at x1 of a spline under c tension (with tension factor sigma) passing through the c three points and having third derivative equal to zero at c x1. optionally, only two values, c1 and c2 are determined. c c on input-- c c del1 is x2-x1 (.gt. 0.). c c del2 is x3-x1 (.gt. 0.). if n .eq. 2, this parameter is c ignored. c c sigma is the tension factor. c c and c c n is a switch indicating the number of coefficients to c be returned. if n .eq. 2 only two coefficients are c returned. otherwise all three are returned. c c on output-- c c c1, c2, and c3 contain the coefficients. c c none of the input parameters are altered. c c this subroutine references package module snhcsh. c c----------------------------------------------------------- c if (n .eq. 2) go to 2 if (sigma .ne. 0.) go to 1 del = del2-del1 c c tension .eq. 0. c c1 = -(del1+del2)/(del1*del2) c2 = del2/(del1*del) c3 = -del1/(del2*del) return c c tension .ne. 0. c 1 call snhcsh (dummy,coshm1,sigma*del1,1) call snhcsh (dummy,coshm2,sigma*del2,1) delp = sigma*(del2+del1)/2. delm = sigma*(del2-del1)/2. call snhcsh (sinhmp,dummy,delp,-1) call snhcsh (sinhmm,dummy,delm,-1) denom = coshm1*(del2-del1)-2.*del1*delp*delm* * (1.+sinhmp)*(1.+sinhmm) c1 = 2.*delp*delm*(1.+sinhmp)*(1.+sinhmm)/denom c2 = -coshm2/denom c3 = coshm1/denom return c c two coefficients c 2 c1 = -1./del1 c2 = -c1 return end subroutine snhcsh (sinhm,coshm,x,isw) c integer isw real sinhm,coshm,x c c coded by alan kaylor cline c from fitpack -- january 26, 1987 c a curve and surface fitting package c a product of pleasant valley software c 8603 altus cove, austin, texas 78759, usa c c this subroutine returns approximations to c sinhm(x) = sinh(x)/x-1 c coshm(x) = cosh(x)-1 c and c coshmm(x) = (cosh(x)-1-x*x/2)/(x*x) c with relative error less than 4.0e-14. c c on input-- c c x contains the value of the independent variable. c c isw indicates the function desired c = -1 if only sinhm is desired, c = 0 if both sinhm and coshm are desired, c = 1 if only coshm is desired, c = 2 if only coshmm is desired, c = 3 if both sinhm and coshmm are desired. c c on output-- c c sinhm contains the value of sinhm(x) if isw .le. 0 or c isw .eq. 3 (sinhm is unaltered if isw .eq.1 or isw .eq. c 2). c c coshm contains the value of coshm(x) if isw .eq. 0 or c isw .eq. 1 and contains the value of coshmm(x) if isw c .ge. 2 (coshm is unaltered if isw .eq. -1). c c and c c x and isw are unaltered. c c----------------------------------------------------------- c data sp14/.227581660976348e-7/, * sp13/.612189863171694e-5/, * sp12/.715314759211209e-3/, * sp11/.398088289992973e-1/, * sq12/.206382701413725e-3/, * sq11/-.611470260009508e-1/, * sq10/.599999999999986e+1/ data sp25/.129094158037272e-9/, * sp24/.473731823101666e-7/, * sp23/.849213455598455e-5/, * sp22/.833264803327242e-3/, * sp21/.425024142813226e-1/, * sq22/.106008515744821e-3/, * sq21/-.449855169512505e-1/, * sq20/.600000000268619e+1/ data sp35/.155193945864942e-9/, * sp34/.511529451668737e-7/, * sp33/.884775635776784e-5/, * sp32/.850447617691392e-3/, * sp31/.428888148791777e-1/, * sq32/.933128831061610e-4/, * sq31/-.426677570538507e-1/, * sq30/.600000145086489e+1/ data sp45/.188070632058331e-9/, * sp44/.545792817714192e-7/, * sp43/.920119535795222e-5/, * sp42/.866559391672985e-3/, * sp41/.432535234960858e-1/, * sq42/.824891748820670e-4/, * sq41/-.404938841672262e-1/, * sq40/.600005006283834e+1/ data cp5/.552200614584744e-9/, * cp4/.181666923620944e-6/, * cp3/.270540125846525e-4/, * cp2/.206270719503934e-2/, * cp1/.744437205569040e-1/, * cq2/.514609638642689e-4/, * cq1/-.177792255528382e-1/, * cq0/.200000000000000e+1/ data zp4/.664418805876835e-8/, * zp3/.218274535686385e-5/, * zp2/.324851059327161e-3/, * zp1/.244515150174258e-1/, * zq2/.616165782306621e-3/, * zq1/-.213163639579425e0/, * zq0/.240000000000000e+2/ c ax = abs(x) if (isw .ge. 0) go to 5 c c sinhm approximation c if (ax .gt. 3.9) go to 2 xs = ax*ax if (ax .gt. 2.2) go to 1 c c sinhm approximation on (0.,2.2) c sinhm = xs*((((sp14*xs+sp13)*xs+sp12)*xs+sp11)*xs+1.)/ . ((sq12*xs+sq11)*xs+sq10) return c c sinhm approximation on (2.2,3.9) c 1 sinhm = xs*(((((sp25*xs+sp24)*xs+sp23)*xs+sp22)*xs+sp21) . *xs+1.)/((sq22*xs+sq21)*xs+sq20) return 2 if (ax .gt. 5.1) go to 3 c c sinhm approximation on (3.9,5.1) c xs = ax*ax sinhm = xs*(((((sp35*xs+sp34)*xs+sp33)*xs+sp32)*xs+sp31) . *xs+1.)/((sq32*xs+sq31)*xs+sq30) return 3 if (ax .gt. 6.1) go to 4 c c sinhm approximation on (5.1,6.1) c xs = ax*ax sinhm = xs*(((((sp45*xs+sp44)*xs+sp43)*xs+sp42)*xs+sp41) . *xs+1.)/((sq42*xs+sq41)*xs+sq40) return c c sinhm approximation above 6.1 c 4 expx = exp(ax) sinhm = (expx-1./expx)/(ax+ax)-1. return c c coshm and (possibly) sinhm approximation c 5 if (isw .ge. 2) go to 7 if (ax .gt. 2.2) go to 6 xs = ax*ax coshm = xs*(((((cp5*xs+cp4)*xs+cp3)*xs+cp2)*xs+cp1) . *xs+1.)/((cq2*xs+cq1)*xs+cq0) if (isw .eq. 0) sinhm = xs*((((sp14*xs+sp13)*xs+sp12) . *xs+sp11)*xs+1.)/((sq12*xs+sq11)*xs+sq10) return 6 expx = exp(ax) coshm = (expx+1./expx)/2.-1. if (isw .eq. 0) sinhm = (expx-1./expx)/(ax+ax)-1. return c c coshmm and (possibly) sinhm approximation c 7 xs = ax*ax if (ax .gt. 2.2) go to 8 coshm = xs*((((zp4*xs+zp3)*xs+zp2)*xs+zp1)*xs+1.)/ . ((zq2*xs+zq1)*xs+zq0) if (isw .eq. 3) sinhm = xs*((((sp14*xs+sp13)*xs+sp12) . *xs+sp11)*xs+1.)/((sq12*xs+sq11)*xs+sq10) return 8 expx = exp(ax) coshm = ((expx+1./expx-xs)/2.-1.)/xs if (isw .eq. 3) sinhm = (expx-1./expx)/(ax+ax)-1. return end C SUBROUTINE QSF C C PURPOSE C TO COMPUTE THE VECTOR OF INTEGRAL VALUES FOR A GIVEN C EQUIDISTANT TABLE OF FUNCTION VALUES. C C USAGE C CALL QSF (H,Y,Z,NDIM) C C DESCRIPTION OF PARAMETERS C H - THE INCREMENT OF ARGUMENT VALUES. C Y - THE INPUT VECTOR OF FUNCTION VALUES. C Z - THE RESULTING VECTOR OF INTEGRAL VALUES. Z MAY BE C IDENTICAL WITH Y. C NDIM - THE DIMENSION OF VECTORS Y AND Z. C C REMARKS C NO ACTION IN CASE NDIM LESS THAN 3. C C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED C NONE C C METHOD C BEGINNING WITH Z(1)=0, EVALUATION OF VECTOR Z IS DONE BY C MEANS OF SIMPSONS RULE TOGETHER WITH NEWTONS 3/8 RULE OR A C COMBINATION OF THESE TWO RULES. TRUNCATION ERROR IS OF C ORDER H**5 (I.E. FOURTH ORDER METHOD). ONLY IN CASE NDIM=3 C TRUNCATION ERROR OF Z(2) IS OF ORDER H**4. C FOR REFERENCE, SEE C (1) F.B.HILDEBRAND, INTRODUCTION TO NUMERICAL ANALYSIS, C MCGRAW-HILL, NEW YORK/TORONTO/LONDON, 1956, PP.71-76. C (2) R.ZURMUEHL, PRAKTISCHE MATHEMATIK FUER INGENIEURE UND C PHYSIKER, SPRINGER, BERLIN/GOETTINGEN/HEIDELBERG, 1963, C PP.214-221. C C .................................................................. C SUBROUTINE QSF(H,Y,Z,NDIM) C C DIMENSION Y(*),Z(*) C HT=.3333333*H IF(NDIM-5)7,8,1 C C NDIM IS GREATER THAN 5. PREPARATIONS OF INTEGRATION LOOP 1 SUM1=Y(2)+Y(2) SUM1=SUM1+SUM1 SUM1=HT*(Y(1)+SUM1+Y(3)) AUX1=Y(4)+Y(4) AUX1=AUX1+AUX1 AUX1=SUM1+HT*(Y(3)+AUX1+Y(5)) AUX2=HT*(Y(1)+3.875*(Y(2)+Y(5))+2.625*(Y(3)+Y(4))+Y(6)) SUM2=Y(5)+Y(5) SUM2=SUM2+SUM2 SUM2=AUX2-HT*(Y(4)+SUM2+Y(6)) Z(1)=0. AUX=Y(3)+Y(3) AUX=AUX+AUX Z(2)=SUM2-HT*(Y(2)+AUX+Y(4)) Z(3)=SUM1 Z(4)=SUM2 IF(NDIM-6)5,5,2 C C INTEGRATION LOOP 2 DO 4 I=7,NDIM,2 SUM1=AUX1 SUM2=AUX2 AUX1=Y(I-1)+Y(I-1) AUX1=AUX1+AUX1 AUX1=SUM1+HT*(Y(I-2)+AUX1+Y(I)) Z(I-2)=SUM1 IF(I-NDIM)3,6,6 3 AUX2=Y(I)+Y(I) AUX2=AUX2+AUX2 AUX2=SUM2+HT*(Y(I-1)+AUX2+Y(I+1)) 4 Z(I-1)=SUM2 5 Z(NDIM-1)=AUX1 Z(NDIM)=AUX2 RETURN 6 Z(NDIM-1)=SUM2 Z(NDIM)=AUX1 RETURN C END OF INTEGRATION LOOP C 7 IF(NDIM-3)12,11,8 C C NDIM IS EQUAL TO 4 OR 5 8 SUM2=1.125*HT*(Y(1)+Y(2)+Y(2)+Y(2)+Y(3)+Y(3)+Y(3)+Y(4)) SUM1=Y(2)+Y(2) SUM1=SUM1+SUM1 SUM1=HT*(Y(1)+SUM1+Y(3)) Z(1)=0. AUX1=Y(3)+Y(3) AUX1=AUX1+AUX1 Z(2)=SUM2-HT*(Y(2)+AUX1+Y(4)) IF(NDIM-5)10,9,9 9 AUX1=Y(4)+Y(4) AUX1=AUX1+AUX1 Z(5)=SUM1+HT*(Y(3)+AUX1+Y(5)) 10 Z(3)=SUM1 Z(4)=SUM2 RETURN C C NDIM IS EQUAL TO 3 11 SUM1=HT*(1.25*Y(1)+Y(2)+Y(2)-.25*Y(3)) SUM2=Y(2)+Y(2) SUM2=SUM2+SUM2 Z(3)=HT*(Y(1)+SUM2+Y(3)) Z(1)=0. Z(2)=SUM1 12 RETURN END C This subroutine gives eastward (u), northward (v), and temperature (T) C at a given pressure height (zp) and local time (time) at 18 degrees N. C geographic latitude for a modified version of the Forbes [1982a,b] (parts I C and II) model....modified according to some average non-winter profiles C given by Harper [1981]. The main modifications are in phase of semidiurnal C and diurnal winds above 110 km (to earlier times), and multiplication of the C diurnal wind amplitudes by 0.8 . The following tables give wind values C from zp=-15.0 (approx 50 km) to zp=+5.0 (approx 400 km) in steps of zp=1.0. C Interpolation is done using spline under tension. subroutine tide18v(zp,mz,time,us,vs,ts) dimension zp(mz),us(mz),vs(mz),ts(mz) dimension ua(2,21),va(2,21),ta(2,21) dimension up(2,21),vp(2,21),tp(2,21) dimension store(2,6,21),phz(2),am(2) c the v's in the following data statements are southward data ( (store(1,j,i), j=1,6), i=1,10)/ & 5.,22., 5., 4., 2.,20., & 7.,16., 11.,22., 2.,18., & 10.,10., 18.,15., 2.,15.4, & 15.,4., 26.,10., 2.,12.3, & 25.,21., 40.,6., 2.,8., & 32.,16., 49.,2., 2.,2.7, & 38.,12., 54.,19., 2.,19., & 45.,7., 60.,13., 2.,11.7, & 50.,3., 70.,8., 2.,5., & 50.,23., 78.,3.8, 3.6,0./ data ( (store(1,j,i), j=1,6), i=11,21)/ & 38.,12., 72.,0., 3.8,15., & 30.,5., 37.,12., 4.,7., & 35.,23., 34.,4., 37.,0., & 42.,20., 35.,24., 75.,14.5, & 52.,17., 40.,20.7, 97.,14.5, & 58.,16.5, 43.,20., 120.,14.5, & 57.,17.2, 41.,21.0, 128.,14.5, & 55.,18.6, 40.,22.0, 136.,14.5, & 55.,19.1, 40.,23.5, 142.,14.5, & 55.,19.6, 40.,23.0, 150.,14.5, & 55.,19.6, 40.,23.0, 150.,14.5/ data ( (store(2,j,i), j=1,6), i=1,10)/ & 2.,9.5, 1.,2.7, 1.,2.7, & 1.7,9.0, 2.,.4, 1.,3., & 2.5,8.6, 3.,12., 1.,3.1, & 3.,8.5, 4.,11.7, 1.,3.2, & 5.,8.4, 7.5,11.5, 1.,3.3, & 5.,8.2, 9.,10.5, 1.,3.2, & 3.,8.0, 8.,8., 3.,3.1, & 4.,8., 11.,6., 3.5,2.6, & 9.,8.7, 25.,4.5, 2.5,1.5, & 17.,9.7, 47.,3., 3.,0./ data ( (store(2,j,i), j=1,6), i=11,21)/ & 22.,8.5, 54.,0.0, 12.5,11., & 37.,5.5, 54.,9.0, 25.,10.1, & 47.,3.0, 52.,6., 20.,9.2, & 50.,2.5, 47.,4.0, 17.,8.4, & 42.,0.7, 38.,2.5, 23.,7.4, & 30.,11.9, 35.,1.4, 34.,6.5, & 27.,10.8, 35.,0.5, 44.,6.2, & 26.,10.1, 35.,0.0, 49.,6.0, & 25.,9.8, 35.,11.5, 52.,6.0, & 25.,9.8, 35.,11.5, 55.,6.0, & 25.,9.8, 35.,11.5, 55.,6.0/ am(1)=0.8 am(2)=1.0 phz(1)=0. phz(2)=0. do k=1,2 do i=1,21 ua(k,i)=store(k,1,i)*am(k) up(k,i)=store(k,2,i)+phz(k) va(k,i)=store(k,3,i)*am(k) vp(k,i)=store(k,4,i)+phz(k) ta(k,i)=store(k,5,i) tp(k,i)=store(k,6,i) end do end do call interpv(ua,up,time,zp,mz,us) call interpv(va,vp,time,zp,mz,vs) call interpv(ta,tp,time,zp,mz,ts) c convert v to northward do i=1,mz vs(i)=-vs(i) end do return end subroutine interpv(fa,fp,time,zp,mz,fs) dimension fa(2,21),fp(2,21),ft(21),fs(mz),zz(21),zp(mz) dimension yp(21),temp(21) d=.2617992 s=.5235983 t=time/3600. do i=1,21 C arbitrary height scale (0 to 20 instead of -15 to 5), so we C can have positive argument in interpolation zz(i)=(i-1)*1.0 ft(i)=fa(1,i)*cos(d*(t-fp(1,i)))+fa(2,i)*cos(s*(t-fp(2,i))) end do sigma=2.0 call curv1(21,zz,ft,slp1,slpn,3,yp,temp,sigma,ierr) do i=1,mz C arbitrary height scale again zpr=zp(i)+15. fzpr=curv2(zpr,21,zz,ft,yp,sigma) fs(i)=fzpr end do return end function intrvl (t,x,n) c integer n real t,x(n) c c coded by alan kaylor cline c from fitpack -- january 26, 1987 c a curve and surface fitting package c a product of pleasant valley software c 8603 altus cove, austin, texas 78759, usa c c this function determines the index of the interval c (determined by a given increasing sequence) in which c a given value lies. c c on input-- c c t is the given value. c c x is a vector of strictly increasing values. c c and c c n is the length of x (n .ge. 2). c c on output-- c c intrvl returns an integer i such that c c i = 1 if e t .lt. x(2) , c i = n-1 if x(n-1) .le. t , c otherwise x(i) .le. t .le. x(i+1), c c none of the input parameters are altered. c c----------------------------------------------------------- c save i data i /1/ c tt = t c c check for illegal i c if (i .ge. n) i = n/2 c c check old interval and extremes c if (tt .lt. x(i)) then if (tt .le. x(2)) then i = 1 intrvl = 1 return else il = 2 ih = i end if else if (tt .le. x(i+1)) then intrvl = i return else if (tt .ge. x(n-1)) then i = n-1 intrvl = n-1 return else il = i+1 ih = n-1 end if c c binary search loop c 1 i = (il+ih)/2 if (tt .lt. x(i)) then ih = i else if (tt .gt. x(i+1)) then il = i+1 else intrvl = i return end if go to 1 end