CDIR$ NOLIST PROGRAM GWAVE2D CDIR$ BOUNDS SAVE c c This program is used to illustrate the linear wave produced by a c gaussian wave form in time. Note that the horizontal is periodic c with fixed wavelength scale PARAMETER (MX=21,ITMX=21,MZ=96,KMX=96,MT=50,NV=2,KMXP=96) PARAMETER (ITX=12) c PARAMETER (ITX=288) PARAMETER (IMXP2=MX+1) C include "blank.h" include "solars.h" include "sodium.h" COMMON/TRAVLH/RK(50),ALP(20),RB(20),RKM(60) COMMON/NOZTRF/XNNOZ(KMX) COMMON/EDDTRF/EDY(KMX) COMMON/WINDTRF/UN(KMX),VN(KMX),WN(KMX),UI(KMX),VI(KMX),WI(KMX) 1,WGAM(KMX),DWGAM(KMX),UIT,VUT COMMON/INDEXX/ITIME,IAUR,ISOLPRO COMMON/PERTRF/TNO(KMX),TNP(KMX),FDIFP(KMX) COMMON/IONTRFF/XNPI(KMX),XNNI(KMX),XNEEE(KMX) C DIMENSION ZPP(MZ),WOUT(7,MX,MZ),WOUTT(5,MX,MZ) c DIMENSION ZPP(MZ),X(MX),WOUT(7,MX,MZ),WOUTT(5,MX,MZ) c DIMENSION ZPP(MZ),X(MX),WOUT(7,MX,MZ),WOUTT(5,MX,MZ), c 1SHTO(MZ) C common/plt4/ AWGAM(ITX,KMX),AVGAM(ITX,KMX),AXNASP(ITX,KMX), 1 AXNAS(ITX,KMX),TIMES(ITX),ZHTT(ITX,KMX) dimension fplt(itx,kmx,4),vp(4) equivalence(fplt,awgam) character*8 flab(4) character*56 toplab data flab /'WGAM ','VGAM ','XNASP ','XNAS '/ c COMMON/TIMETRF/TIMEE common/iyrday/IYD COMMON/OXOXOX2/O3OR(KMX,-1:IMXP2),XNOX(KMX), 1 OOXR(KMX,-1:IMXP2),FNOX(KMX,-1:IMXP2),XNOX1, 1 PS2B(-1:IMXP2) COMMON /PARA/ CP,RG,GA,H,G,TEMPO,DTO,TO,WO, < IZM,XKX,XKZ,OMEGO,DR,PI,ALPHA,ALTO,DZP,ZPO < ,DX,XO,IXM C COMMON/NOZNOZ2/SNO2NO(KMX,-1:IMXP2),RNONOZ(KMX,-1:IMXP2), 1 FNNOZ(KMX,-1:IMXP2) DATA DX/400.0e+5/,XO/0.0/,DZP/0.2/ > ,ZPO/-14.0/ c common/ehist/shto(mz),eox(mx,mz),eo(mx,mz),eo3(mx,mz),ehox(mx,mz), + eh(mx,mz),eho2(mx,mz),eoh(mx,mz),enoz(mx,mz),eno(mx,mz), + eno2(mx,mz),en4s(mx,mz),enas(mx,mz),enasp(mx,mz),ene(mx,mz), + etion(mx,mz),eo2(mx,mz),en2(mx,mz),eo1d(mx,mz),ete(mx,mz), + eti(mx,mz),etn(mx,mz),eiop(mx,mz),en2d(mx,mz),enh2o(mx,mz), + enh2(mx,mz),ench4(mx,mz),enh2o2(mx,mz),enco(mx,mz),enco2(mx,mz) + ,enax(mx,mz),enao(mx,mz),enao2(mx,mz),enaoh(mx,mz) c c Arrays in /pcontour/ are for contouring (see plt.f): C C HTINT(NHTINT) = heights at which to interpolate C (set to 66,68,...,198,200 by BF, 5/90) C FINT(MX,NHTINT) = interpolated field for contouring C PARAMETER(NHTINT=68,HT1=66.,DHT=2.) common/pcontour/ ETNN(MX,MZ),UWIND(MX,MZ),VWIND(MX,MZ), + WWIND(MX,MZ),ZHT(MX,MZ), + PKE(MX,MZ),pui(mx,mz),pvi(mx,mz),pwi(mx,mz),pwgam(mx,mz), + POX(MX,MZ),PNE(MX,MZ),PO(MX,MZ),PO3(MX,MZ),PHOX(MX,MZ), + PH(MX,MZ),POH(MX,MZ),PHO2(MX,MZ),PNOZ(MX,MZ),PNO(MX,MZ), + PNO2(MX,MZ),PN4S(MX,MZ),PLNAS(MX,MZ),PLNASP(MX,MZ), + eun(mx,mz),ewn(mx,mz),eke(mx,mz),UWINDZ(MX,MZ), + EDDT(MX,MZ),TNPER(MX,MZ),HTINT(NHTINT),FINT(MX,NHTINT),x(mx) dimension GINT(ITX,NHTINT) COMMON/EDDYTRE/FDIFO(KMX) COMMON/WIND/UU(KMX),WW(KMX) C COMMON/FIELDS2/FT(KMX,-1:IMXP2),FNO2(KMX,-1:IMXP2), 1 FNO(KMX,-1:IMXP2),FNN2(KMX,-1:IMXP2),FNAR(KMX,-1:IMXP2), 2 FNHE(KMX,-1:IMXP2),FNNO(KMX,-1:IMXP2),FNN2D(KMX,-1:IMXP2), 3 FNN4S(KMX,-1:IMXP2),FRJ(KMX,-1:IMXP2),FW(KMX,-1:IMXP2), 4 FNHOX(KMX,-1:IMXP2),PPN2D(KMX,-1:IMXP2),FNE(KMX,-1:IMXP2), 5 PPN4S(KMX,-1:IMXP2),FK4O2P(KMX,-1:IMXP2), 6 FK5O2P(KMX,-1:IMXP2),FF107,FTE(KMX,-1:IMXP2), 7 FDIFK(KMX,-1:IMXP2),FU(KMX,-1:IMXP2),FNAS(KMX,-1:IMXP2), 8 FNAO(KMX,-1:IMXP2),FNAO2(KMX,-1:IMXP2),FNAOH(KMX,-1:IMXP2) C COMMON/FIELDT2/FNH2O(KMX,-1:IMXP2),FNH2O2(KMX,-1:IMXP2), 1 FNO1D(KMX,-1:IMXP2),FNH2(KMX,-1:IMXP2),FNOH(KMX,-1:IMXP2), 2 FNHO2(KMX,-1:IMXP2),FNO3(KMX,-1:IMXP2),FNH(KMX,-1:IMXP2), 3 FNCH4(KMX,-1:IMXP2),FNCO(KMX,-1:IMXP2),FNCO2(KMX,-1:IMXP2), 4 FNNO2(KMX,-1:IMXP2) C COMMON/SODIUMP/PNA(KMXP,-1:IMXP2),XLNA(KMXP,-1:IMXP2), 1XLBNA(-1:IMXP2),PNAI(KMXP,-1:IMXP2),XLNAI(KMXP,-1:IMXP2), 2XLBNAI(-1:IMXP2),RNAO(KMXP,-1:IMXP2),RNAO2(KMXP,-1:IMXP2), 3RNAOH(KMXP,-1:IMXP2) C CALL OPNGKS INUMXY=0 CHARSZ = 0.025 CHSZNEG = -CHARSZ EXPSZ = 0.019 DO 50 I=1,NHTINT 50 HTINT(I) = HT1-DHT+FLOAT(I)*DHT C MXM=MX MZM=MZ F107=70. F107A=70. IYD=76080 DAY=80. STL=12. GLAT=18.3 GLONG=0. RE=6371.E+5 BOLTZ=1.38E-16 ITIME=0 ISOLPRO=0 UT=0. c c If IW < 0 restart will read from abs(iw), otherwise c it will write to iw c IW=-90 c IW=-5 CALL RESTART(UT,IW) DO 15 IX=1,MX X(IX)=XO+(IX-1)*DX 15 CONTINUE DO 25 IZ=1,MZ ZP(IZ)=(IZ-1)*DZP+ZPO ZPP(IZ)=ZP(IZ) 25 CONTINUE DO 88 J=1,MZ EDDYC=2.0E-6 EDDYMB=8.5E-7 ZPMB=-9.0 EDDYLB=8.0E-8 EZB=-14. TURZP=-6.5 FDIFO(J)=EDDYC*EXP((TURZP-ZP(J))) IF(ZP(J).LE.TURZP) FDIFO(J)=EDDYMB*EXP(ALOG(EDDYC/EDDYMB)*(ZP(J)- 1ZPMB)/(TURZP-ZPMB)) IF(ZP(J).LE.ZPMB) FDIFO(J)=EDDYLB*EXP(ALOG(EDDYMB/EDDYLB)*(ZP(J)- 1ZP(1))/(ZPMB-EZB)) XNAPCO2(J)=1.E-20 XNAPH2O(J)=1.E-20 XNAOP(J)=1.E-20 XNAPN2(J)=1.E-20 88 CONTINUE DTIME=300. C**************************** DO 500 IT=1,ITX C**************************** c write(6,"('500 loop: it=',i3,' itx=',i3)") it,itx TIME=DTIME+(IT-1)*DTIME TIMEE=TIME/3600. TIMES(IT)=TIMEE C ******************************************************* C C RAY YOU ONLY REQIURE KNOWLEDGE OF GRID C THEN A SINGLE CALL TO SUBROUTINE TIDE C C CALL TIDES(TIME,X,MXM,ZPP,MZM,TNO,SHTO,WOUTT) ISEASAV=0 IUTAV=0 GMLAT=29.7 GMLONG=7.8 UTT=TIMEE-4. IF(UTT.LT.0) UTT=UTT+24. IF(UTT.GT.24.) UTT=UTT-24. DAYP=DAY C CALL EFIELD(GMLAT,GMLONG,DAYP,UTT,ISEASAV,IUTAV,POT,VU,VE) C UIT=VE*1.E+2 C VUT=VU*1.E+2 UIT=0. VUT=0. C C C ******************************************************* DO 444 I=1,MX DO 444 J=1,MZ EDDT(I,J)=0. C UWIND(I,J)=WOUTT(1,I,J)*1.E+2 C VWIND(I,J)=WOUTT(2,I,J)*1.E+2 C WWIND(I,J)=WOUTT(3,I,J)*SHTO(J) C TNPER(I,J)=WOUTT(4,I,J) C ZHT(I,J)=ZPHT(J)+WOUTT(5,I,J)*1.E-5 UWIND(I,J)=0. VWIND(I,J)=0. WWIND(I,J)=0. TNPER(I,J)=0. ZHT(I,J)=ZPHT(J) EUN(I,J)=VWIND(I,J) EWN(I,J)=WWIND(I,J) ETN(I,J)=TNO(J)+TNPER(I,J) EKE(I,J)=FDIFO(J)+EDDT(I,J) FNO2(J,I)=EO2(I,J) FNO(J,I)=EO(I,J) FNO3(J,I)=EO3(I,J) FNOX(J,I)=EOX(I,J) FNN2(J,I)=EN2(I,J) FNNOZ(J,I)=ENOZ(I,J) FNNO(J,I)=ENO(I,J) FNNO2(J,I)=ENO2(I,J) FNN2D(J,I)=EN2D(I,J) FNN4S(J,I)=EN4S(I,J) FNCH4(J,I)=ENCH4(I,J) FNH2(J,I)=ENH2(I,J) FNCO2(J,I)=ENCO2(I,J) FNCO(J,I)=ENCO(I,J) FNHOX(J,I)=EHOX(I,J) FNH2O(J,I)=ENH2O(I,J) FNH(J,I)=EH(I,J) FNHO2(J,I)=EHO2(I,J) FNOH(J,I)=EOH(I,J) FNH2O2(J,I)=ENH2O2(I,J) FNAS(J,I)=ENAS(I,J) FNHE(J,I)=ENASP(I,J) FNAR(J,I)=ENAX(I,J) FNAO(J,I)=ENAO(I,J) FNAO2(J,I)=ENAO2(I,J) FNAOH(J,I)=ENAOH(I,J) 444 CONTINUE DO 445 I=1,MX DO 446 J=1,MZ UN(J)=UWIND(I,J) VN(J)=VWIND(I,J) WN(J)=WWIND(I,J) EDY(J)=EKE(I,J) TNP(J)=TNPER(I,J) FDIFP(J)=EDDT(I,J) XNOX(J)=EOX(I,J) XNO(J)=EO(I,J) XNO3(J)=EO3(I,J) XNHOX(J)=EHOX(I,J) XNH(J)=EH(I,J) XNOH(J)=EOH(I,J) XNHO2(J)=EHO2(I,J) XNNOZ(J)=ENOZ(I,J) XNNO(J)=ENO(I,J) XNNO2(J)=ENO2(I,J) XN4S(J)=EN4S(I,J) XNAS(J)=ENAS(I,J) XNASP(J)=ENASP(I,J) XNAX(J)=ENAX(I,J) XNE(J)=ENE(I,J) XNPI(J)=ETION(I,J) XNO2(J)=EO2(I,J) XNN2(J)=EN2(I,J) XNO1D(J)=EO1D(I,J) RHO(J)=(16.*XNO(J)+32.*XNO2(J)+28.*XNN2(J))*1.66E-24 TN(J)=ETN(I,J) TE(J)=ETE(I,J) TI(J)=ETI(I,J) XIOP(J)=EIOP(I,J) XN2D(J)=EN2D(I,J) XNH2O(J)=ENH2O(I,J) XNH2(J)=ENH2(I,J) XNCH4(J)=ENCH4(I,J) XNH2O2(J)=ENH2O2(I,J) XNCO(J)=ENCO(I,J) XNCO2(J)=ENCO2(I,J) UU(J)=VWIND(I,J) WW(J)=WWIND(I,J)/SHTO(J) 446 CONTINUE CALL COLUMN(UT) CALL SOLHEAT(UT) CALL DENMOD(UT,DTIME) C CALL NADIFSV(UT,DTIME) CALL COMPSEU(DTIME,I,IT) DO 447 J=1,MZ ENE(I,J)=XNE(J) ETION(I,J)=XNPI(J) EO1D(I,J)=XNO1D(J) ETE(I,J)=TE(J) ETI(I,J)=TI(J) ETN(I,J)=TN(J) EIOP(I,J)=XIOP(J) EN2D(I,J)=XN2D(J) PUI(I,J)=UI(J) PVI(I,J)=VI(J) PWI(I,J)=WI(J) PWGAM(I,J)=WGAM(J) 447 CONTINUE 445 CONTINUE DS=0.2 SB=-14. CALL FACE2(DTIME,DS,SB,DX) DO 450 I=1,MX DO 450 J=1,MZ EO2(I,J)=FNO2(J,I) EO(I,J)=FNO(J,I) EO3(I,J)=FNO3(J,I) EOX(I,J)=FNOX(J,I) EN2(I,J)=FNN2(J,I) ENOZ(I,J)=FNNOZ(J,I) ENO(I,J)=FNNO(J,I) ENO2(I,J)=FNNO2(J,I) EN4S(I,J)=FNN4S(J,I) ENCH4(I,J)=FNCH4(J,I) ENH2(I,J)=FNH2(J,I) ENCO2(I,J)=FNCO2(J,I) ENCO(I,J)=FNCO(J,I) EHOX(I,J)=FNHOX(J,I) ENH2O(I,J)=FNH2O(J,I) EH(I,J)=FNH(J,I) EHO2(I,J)=FNHO2(J,I) EOH(I,J)=FNOH(J,I) ENH2O2(I,J)=FNH2O2(J,I) ENASP(I,J)=FNHE(J,I) ENAX(I,J)=FNAR(J,I) FNAS(J,I)=FNAR(J,I)/(1.+RNAO(J,I)+RNAO2(J,I)+RNAOH(J,I)) ENAS(I,J)=FNAS(J,I) FNAO(J,I)=FNAS(J,I)*RNAO(J,I) ENAO(I,J)=FNAO(J,I) FNAO2(J,I)=FNAS(J,I)*RNAO2(J,I) ENAO2(I,J)=FNAO2(J,I) FNAOH(J,I)=FNAS(J,I)*RNAOH(J,I) ENAOH(I,J)=FNAOH(J,I) 450 CONTINUE DO 548 J=1,MZ AVGAM(IT,J)=PVI(10,J) AWGAM(IT,J)=PWGAM(10,J) AXNASP(IT,J)=ALOG10(AMAX1(ENASP(10,J),1.E-20)) AXNAS(IT,J)=ALOG10(AMAX1(ENAS(10,J),1.E-20)) ZHTT(IT,J)=ZHT(10,J) 999 FORMAT(1X,I3,2X,5E12.3) 548 CONTINUE C************************ C IF(MOD(IT,12).EQ.0) GO TO 449 c IF(MOD(IT,144).EQ.0) GO TO 449 IF(MOD(IT,4).EQ.0) GO TO 449 GO TO 500 449 CONTINUE C************************ DO 448 I=1,MX DO 448 J=1,MZ PKE(I,J)=ALOG10(EKE(I,J)) POX(I,J)=FNOX(J,I) PO(I,J)=FNO(J,I) PO3(I,J)=ALOG10(FNO3(J,I)) PHOX(I,J)=ALOG10(FNHOX(J,I)) PH(I,J)=ALOG10(FNH(J,I)) POH(I,J)=ALOG10(FNOH(J,I)) PHO2(I,J)=ALOG10(FNHO2(J,I)) PNOZ(I,J)=ALOG10(FNNOZ(J,I)) PNO(I,J)=ALOG10(FNNO(J,I)) PNO2(I,J)=ALOG10(FNNO2(J,I)) C PNE(I,J)=ALOG10(AMAX1(ENE(I,J),1.E-20)) PNE(I,J)=ALOG10(AMAX1(ETION(I,J),1.E-20)) PN4S(I,J)=ALOG10(FNN4S(J,I)) PLNAS(I,J)=ALOG10(AMAX1(FNAS(J,I),1.E-20)) PLNASP(I,J)=ALOG10(AMAX1(FNHE(J,I),1.E-20)) 448 CONTINUE c c Contour at current time: call plt 500 continue c c Contour with time on x-axis, zp on y-axis: call cpset(times,itx,zp,kmx,vp,xmid) finc = -6. write(6,"('gwave2d contour 4 fields with zp on y-axis:', + ' time=',f6.2,' f107=',f5.1,' day=',i5)") timee,f107,iyd do ip=1,4 call cpcnrc(fplt(1,1,ip),itx,itx,kmx,0.,0.,finc,1,-1,-1634B) call labrect(times,itx,zp,kmx,'TIME','ZPRES',0.) call clearstr(toplab) write(toplab,"('FIELD=',a,' TIME=',f5.2,' F107=',f5.1, + ' DAY=',i5)") flab(ip),timee,f107,iyd call wrlab(toplab(1:lenstr(toplab)),xmid,vp(4)+.05,.018) call frame enddo c c Contour with time on x-axis, ht on y-axis: call cpset(times,itx,htint,nhtint,vp,xmid) write(6,"('gwave2d contour 4 fields with ht on y-axis:', + ' time=',f6.2,' f107=',f5.1,' day=',i5)") timee,f107,iyd do ip=1,4 call twodint(fplt(1,1,ip),zhtt,itx,kmx,gint,htint,nhtint,0,0) call cpcnrc(fplt(1,1,ip),itx,itx,kmx,0.,0.,finc,1,-1,-1634B) call labrect(times,itx,htint,nhtint,'TIME','HT (KM)',0.) call clearstr(toplab) write(toplab,"('FIELD=',a,' TIME=',f5.2,' F107=',f5.1, + ' DAY=',i5)") flab(ip),timee,f107,iyd call wrlab(toplab(1:lenstr(toplab)),xmid,vp(4)+.05,.018) call frame enddo c c If IW < 0 restart will read from abs(iw), otherwise c it will write to iw c c IW=1 IW=91 CALL RESTART(UT,IW) CALL CLSGKS STOP END subroutine twodint(fin,fht,nx,ny,fout,hts,nhts,ilog,iprnt) c c Interpolate 2d field fin(nx,ny) to hts(nhts) heights. c Interpolated output field is returned in fout(nx,nhts) c (All heights are in km) c c On input: c fin(nx,ny) = input field to be interpolated c (pressure points are in second dimension) c fht(nx,ny) = height field corresponding to fin c (pressure points are in second dimension) c nx,ny = dimensions of fin and fht (will be 21,96 c for use with Robles twogw model) c hts(nhts) = array of target heights at which fin is to be found c (hts(1) < min(fht) and hts(nhts) < max(fht)) c nhts = number of target heights c ilog = log interpolation flag: if ilog = 1, then do a log c interpolation; if ilog.ne.1, do linear interpolation c iprnt = print flag: if iprnt=1 then print out various stuff, c otherwise twodint works silently unless an error is c detected c c On output: c fout(nx,nhts) = interpolated output field c (interpolation height points in second dimension) c (input fields are unchanged) c dimension fin(nx,ny),fht(nx,ny),fout(nx,nhts),hts(nhts) c c Check input: c if (iprnt.eq.1) write(6,"(' twodint: on entry: nx ny=',2i4, + ' nhts=',i4,' ilog=',i2,' hts=',/(6f13.3))") + nx,ny,nhts,ilog,(hts(i),i=1,nhts) if (nx.le.0.or.ny.le.0) then write(6,"(' >>> twodint error: bad nx ny=',2i6)") nx,ny stop 'twodint' endif if (nhts.le.0) then write(6,"(' >>> twodint error: bad nhts=',i6)") nhts stop 'twodint' endif c c Find min/max of fht and check against hts array: c fhtmax = -1.e36 fhtmin = +1.e36 do 100 ix = 1,nx do 100 iy = 1,ny fhtmax = amax1(fhtmax,fht(ix,iy)) fhtmin = amin1(fhtmin,fht(ix,iy)) 100 continue if (iprnt.eq.1) write(6,"(' twodint: fhtmin=',f12.3,' fhtmax=', + f12.3,' hts(1)=',f12.3,' hts(nhts)=',f12.3)") + fhtmin,fhtmax,hts(1),hts(nhts) if (hts(1).le.fhtmin) then write(6,"(' >>> twodint error: hts(1) below min height: ', + 'hts(1)=',f12.3,' fhtmin=',f12.3)") hts(1),fhtmin stop 'twodint' endif if (hts(nhts).ge.fhtmax) then write(6,"(' >>> twodint error: hts(nhts) above max height: ', + 'hts(nhts)=',f12.2,' fhtmax=',f12.2)") hts(nhts),fhtmax stop 'twodint' endif c c X loop: c do 200 ix = 1,nx c c Height loop: c do 225 ih = 1,nhts c c hts array should be monotonically increasing: c if (ih.gt.1.and.hts(ih).lt.hts(ih-1)) then write(6,"(' >>> twodint error: hts not monotonically ', + 'increasing: ih=',i3,' hts(ih)=',f12.2,' hts(ih-1)=', + f12.2)") ih,hts(ih),hts(ih-1) stop 'twodint' endif c c Find indices in fht that bracket desired height: c (fht(ix,ialt1) below hts(ih), and fht(ix,ialt2) above) c do 250 iy = 1,ny-1 if (hts(ih).ge.fht(ix,iy).and.hts(ih).le.fht(ix,iy+1)) then ialt1 = iy ialt2 = iy+1 goto 255 endif 250 continue write(6,"(' >>> twodint error: could not find index ', + 'in fht to hts: ix=',i3,' ih=',i3,' hts(ih)=',f12.2)") + ix,ih,hts(ih) stop 'twodint' 255 continue if (iprnt.eq.1.and.mod(ix,7).eq.0) + write(6,"(' twodint: ix=',i3,' ih=',i3, + ' hts(ih)=',f12.3,' ialt1=',i3,' fht(ix,ialt1)=',f12.3, + ' ialt2=',i3,' fht(ix,ialt2)=',f12.3)") + ix,ih,hts(ih),ialt1,fht(ix,ialt1),ialt2,fht(ix,ialt2) c c Do log interpolation: c if (ilog.eq.1) then rlogarg = fin(ix,ialt2) / fin(ix,ialt1) if (rlogarg.le.0.) then write(6,"(' >>> twodint error: rlogarg <= 0: rlogarg=', + e12.4,' ix=',i3,' ih=',i3,' hts(ih)=',f12.3,' ialt1,2=', + 2i3,/' fht(ix,ialt1)=',f12.3,' fht(ix,ialt2)=',f12.3, + /' fin(ix,ialt1)=',e12.4,' fin(ix,ialt2)=',e12.4)") + rlogarg,ix,ih,hts(ih),ialt1,ialt2,fht(ix,ialt1), + fht(ix,ialt2),fin(ix,ialt1),fin(ix,ialt2) stop 'rlogarg' endif exparg = (alog(fin(ix,ialt2) / fin(ix,ialt1)) / + (fht(ix,ialt2) - fht(ix,ialt1))) * + (hts(ih) - fht(ix,ialt1)) fout(ix,ih) = fin(ix,ialt1) * exp(exparg) if (iprnt.eq.1.and.mod(ix,7).eq.0) + write(6,"(' twodint log interp: ix=',i3, + ' ih=',i3,' hts(ih)=',f12.3,' fout(ix,ih)=',e12.4)") + ix,ih,hts(ih),fout(ix,ih) else c c Do linear interpolation: c f1 = (hts(ih)-fht(ix,ialt1)) / (fht(ix,ialt2)-fht(ix,ialt1)) fout(ix,ih) = f1*fin(ix,ialt2) + (1.-f1)*fin(ix,ialt1) if (iprnt.eq.1.and.mod(ix,7).eq.0) + write(6,"(' twodint linear interp: ix=',i3, + ' ih=',i3,' hts(ih)=',f12.3,' fout(ix,ih)=',e12.4)") + ix,ih,hts(ih),fout(ix,ih) endif c c End height loop: c 225 continue c c End x loop: 200 continue return end