cCc This is SWAMPRA in range dependent form ccccccccccccccccccccccc cccccc for Unix environment (adapted by graduate student Mallikarjuna Kilari) cccccc COMPILE with f95 -o outexe *.f cccccc swamp (shallow water acoustic mode program) ccccccccccccccccc cccccc Developed at NRL-SSC Numerical modeling Code 7181 ccccccccccccc cccccc Sponsered by NRL Ed Franchi & 6.1 and ONR 6.1 Jeff Simmins cccc cccccc Last modified Nov. 15, 1996 Integration routine speeded up ccc cccccc This code has up-dated routines MAKECU and MAKECD ccccccccccccc cccccc It also has an adaptive algorithm to determine optimal step ccc cccccc sizes dr. The only arbitrary parameter has DNRA has an cccccccc cccccc effect on step size dr. It is defaulted to 0.5 if you choose cc cccccc anything less then 0.5. It determined the maximim percentage cc cccccc in depth change dh which determines in part the step size dr. c cccccc A value of 1.0 means that dh is allowed to change 1% of the cccc cccccc water column depth- h- at that point. A value of 2.0 would be c cccccc 1/DNRA % or 1/2 % or 1/200. A value of 4 would mean 1/400. cccc cccccc This only assures that the % change will be no greater than ccc cccccc 1/DNRA %. It could even be less depending on the local rate ccc cccccc of change of bathymetry, the number of modes and the cccccccccc cccccc frequency. The adaptive step algotithm is a function of water c cccccc column depth, local rate of change of bathymetry,number of cccc cccccc modes, frequency and bottom hardness. These models work ccccccc cccccc better down-slope than up slope since flux loss to the bottom c cccccc is greater and less predictable going up slope. The addaptive c cccccc step size algorithm does not necessarily decrease dr with ccccc cccccc increasing frequency since the allowed step size dr is also ccc cccccc inversely related to the number of modes because dr is also ccc cccccc a function of the rate of change in mode number because ccccccc cccccc a larger number of modes tolerate greater changes in dr. DNRA c cccccc values of 1.0 are conservative,a value of 2.0 is very ccccccccc cccccc conservative in favour of precision and a value of 4.0 is cccc cccccc very likey overkill. Increasing DNRA by two increases the run c cccccc time by almost 2.0. Larger values of DNRA are preferred up cccc cccccc slope than down slope. All other paremeters in this code are cc cccccc physical or output parameters. Also the value delr determines c cccccc the minimum value of dr that will be tolerated. ccccccccccccccc cccccc The physical parameters are as follows: cccccccccccccccccccccc cccccc h- the initial water column hight from rrr1 to rstrt which is c cccccc assumed constant. cccccc hhh- the final water column height between rall to rend which cc cccccc is then assumed to be constant between rall and rend ccccc cccccc hs- is the source depth. cccccccccccccccccccccccccccccccccccccc cccccc hr- is the reciever depth. If hr is =0 then the reciever ccccc cccccc depths are set at a spacing of dh from the surface to a ccc cccccc place in the water column or the bottom. cccccccccccccccccc cccccc dh- when hr=0 then there are multiple reciever depths at cccccc cccccc spacings from dh from the water surface to the bottom. cccc cccccc dh should be smaller than the smallest value of the water c cccccc column height in along the bathemetry. It is ignored when cc cccccc hr is not set to 0.0. ccccccccccccccccccccccccccccccccccccc cccccc hhh- the water column height at rall. cccccc rstrt- the value after the first calulated value of r=rrr1 cccc cccccc at which the bathymetry changes. rstrt> rrr1. The cccc cccccc water column hight is h between rrr1 and rstrt ccccccccc cccccc rall- the final point in range for which the bathymetry cccccc cccccc changes. r varies from h at rstrt to hhh at rall ccccccc cccccc rall>rstrt. cccccccccccccccccccccccccccccccccccccccccccc cccccc rend- The final value of r for which the bathymetry is ccccccc cccccc constant. rend>rall. ccccccccccccccccccccccccccccccccccc cccccc freq- It is the frequency of the monochromatic point source cc cccccc rho0- Density in the water. Usually set to 1.0. cccccccccccccc cccccc rho1- Density is the bottom half-space rho1>rho0 usually. cccc cccccc c11- Constant sound speed in the bottom half-space. c11> cccc cccccc the minimum or mean sound speed in the water column. ccc cccccc alpha- the loss in the bottom half-space in dB/wave length. ccc cccccc y(i)- the value of the sound speed at depth z(i) ccccccccccccc cccccc z(i)- the depth measured from the top of the water column in c cccccc the half-space. These pair of values y(i) and z(i) ccccc cccccc must be input for the largest values of the water cccccc cccccc column depth which may be h or hhh in the present form.c cccccc Termination of reading the pair y(i) and z(i) is done cc cccccc via -1 -1. cccccc The following are control parameters cccccccccccccccccccccccccc cccccc nwhat- determined the order of the velocity profile in the cccc cccccc water column. When nwhat=0 then it reads y(i),z(i). cccc cccccc When nwhat>0 then it reads the pair as z(i),y(i) ccccccc cccccc rrr1- the first print out value of the transmission loss ccccc cccccc rrr1 may be anything between 0.0 to rstrt-delrs cccccccc cccccc delrs- this is the step side for print out. It is wise to ccccc cccccc choose delrs so that rrr1, rstrt, rall and rend are cccc cccccc integral multiples of delrs. ccccccccccccccccccccccccccc cccccc delrm- the minumum tolerated step size allowed for dr. cccccccc cccccc Although delr is set in the code along bathemetry you cc cccccc judge that delr shuold never go below some value delrm.c cccccc It is wise to choose delrm to be an integral fraction cc cccccc of delrs. cccccccccccccccccccccccccccccccccccccccccccccc ccccc Note it is easy to generalize to more diverse bathymetry. We ccc ccccc are doing that presently. Further a layered model is almost cccc ccccc completed and will be issued when tests are complete. This ccccc ccccc has been bench marked against the PE code RAM by Michael ccccccc ccccc Collins for several difficult up-slope and down-slope problems.c ccccc RAM is a magnificent and impressive code. There are many ccccccc ccccc precision parameters in RAM that determine the run time which cc ccccc include the number of pade' terms ( must be set > 2 and cccccccc ccccc probably > 6.0. The step size in depth dz. If the pressure ccccc ccccc release bottom is set not far below the ocean bottom then cccccc ccccc dz should be set at least 1/10 th of the wave number. For cccccc ccccc deeper bottoms and with greater impedance it ought to be even cc ccccc smaller. Typically dz should be set at about 75/freq to be ccccc ccccc safe or about 1/20 of a wave length. The reference value of cccc ccccc sound velocity in RAM may be chosen as the mean value in the ccc ccccc water column. Perhaps if the half-space sound speed is large ccc ccccc then the mean value of the mean value of the water column and cc ccccc that in the bottom would be a safe value. cccccccccccccccccccccc ccccc The step size dr for RAM depends on the rate of change of cccccc ccccc bathemetry as well as restrictions on the split step Pade'cccccc ccccc method used by Collins. One may vary dr to obtain stability cccc ccccc in TL but our experience is that much larger step sizes may be c ccccc used in RAM than in the PE code UMPE which places RAM at a ccccc ccccc much favourable advantage over UMPE. Further UMPE is a poor cccc ccccc shallow water model due to Gibbs phenomena for any reasonable cc ccccc bottom parameters and because variable density is assumed to cc ccccc commute with other physical parameters which is a bad cccccccccc ccccc approximation for any reasonable bottom. Thus, RAM is a much ccc ccccc better shallow water model, allows for much higher order ccccccc ccccc angles than UMPE and even with high precision parameters it is c ccccc more favourable in run time. ccccccccccccccccccccccccccccccccccc ccccc Since SWAMPR uses an adaptive step size determinator it is ccccc ccccc some what faster than RAM when we have run RAM using ccccccccccc ccccc describes above for RAM. For a sea mount problem it took SWAMP c ccccc 81 sec. to go over a range of 50 KM which includes going up a cc ccccc steep sea mount. For equivalent accuracy it took RAM over 1330 c ccccc sec. which is an advantage of a factor of 16. SWAMPR slows ccccc ccccc down in run time at a greater rate than RAM in it's present cccc ccccc form. For both codes they speed up appreciably as the cccccccccc ccccc rate of change of bathymetry decreases. RAM may incorporate cccc ccccc an adaptive step size algorithm for which case it will compete c ccccc better with SWAMPR. SWAMPR may output modal group velocites, ccc ccccc modal eigenvalues, model attenuation coefficients, etc. which cc ccccc helps in the analysis of data. SWAMPR and probably RAM have cccc ccccc potential to run faster subject to improved algorthms which are c ccccc ongoing. Also SWAMPR allows for expressing the individual cccccc ccccc eigenfunctions in a spherical representation required for cccccc ccccc scattering in a wave guide. cccccccccccccccccccccccccccccccccccc ccccc The next updates to SWAMPR are more versatile input parameters c ccccc which is a trivial matter and the addition of multiple ocean ccc ccccc bottom layering which is tedious but straitfoward and will cccc ccccc not affect run time greatly. A version which allows for the cccc ccccc inclusion of bottom layers with elastic boundary and/or viscous c ccccc parameters will be developed pending demand and funding. ccccccc ccccc SWAMP and SWAMPR are meant to be high frequency models and ccccc ccccc SWAMP has been run with over 1000 modes. The conversion of ccccc ccccc of number of modes to frequency is roughly N=0.0004*h*f. That cc ccccc means that f*h=2500*N or for 1000 modes in 100 meters of water c ccccc f=25 kHz. In 10 meters it is 250 kHz. That is not to say SWAMP ccc ccccc and SWAMPR cann't run for N>1000 but SWAMP ran well on a 1992 cccc ccccc vintage SGI with no problem. We have performed time domain ccccc ccccc calculations on a PC-Pentium 133 with 128 MB of memory and a ccc ccccc SCSI large disk with a mean number of modes of 100 for 1600 cccc ccccc frequencies in about 81 min. That should translate to about 8 cc ccccc minutes on an optimized SGI with an R8000 cpu and on the order c ccccc of seconds on a modern CRAY. SWAMPR is based on matrix and ccccc ccccc iterative methods and should be easy to vectorize or addopt to c ccccc a pipe line computer. cccccccccccccccccccccccccccccccccccccccccc ccccc It runs on a PC in Micro-soft power fortron, on an SGI and on cc ccccc a DEC alpha. It has been run on a CRAY YMP at a very fast cccccc ccccc rate as well as on an SGI 8 processor Chellenger at great speeds. cccccccccccccc SWAMPRA with new matrix method Generalized eigenvalue problen ccccc implicit real*8(a-h,o-z) parameter (np=1000,ny=350) common/blkmain1/h,hs,hr,hb,rr1,freq,rho0,rho1,c00,c11,alpha common/blkmain3/c2(np) common/blockncl/n2,b2,g2,atn2,atn2b common/blockma/n02,akz2,akr2,a2 common/blksico/sn2(np),cs2(np) character*20 DTG complex*16 c1,c2,sum1,ci real*8 freq0 integer npp,nos,Nf,iib,nnn dimension c1(np),akr1(np),akr2(np),akz1(np),akz2(np),sn1(np) dimension a1(np,np),a2(np,np),g1(np),g2(np),b1(np),b2(np),cs1(np) dimension y(ny),z(ny),ya(ny),za(ny),yaa(ny),zaa(ny),atn1b(np) dimension aslp(ny),atn1(np),atn2(np),atn2b(np) CALL DATIME(DTG) WRITE ( 6, '(1X,''START : '',A/)' ) DTG open(unit=24,file='SWAMPRA.OUT',form='formatted',status='unknown') open(unit=5,status='old',file='pulse.in') open(unit=4,status='replace',file='swampgreen.dat', &RECL=8,form='binary') read(5,*)freq0,npp,nos,Nf c freq1=400 c freqn=600 freq1=freq0-freq0/nos freqn=freq0+freq0/nos dfreq=freq0*npp/Nf iib=(-1) nnr1=0 do 3333 ii=0,Nf/2,1 freqr=real(ii)*dfreq if(freqr.lt.freq1) goto 3333 if(freqr.gt.freqn) goto 3333 open(unit=28,file='SWAMPRA.IN',form='formatted',status='old') iib=iib+1 print*,ii,freqr,dfreq rewind(28) read(28,*)nwhat,h,hhh,hs,hr,freq read(28,*)rho0,rho1,c11,alpha read(28,*)dh,rstrt,rall,rend read(28,*)rrr1,delrm,delrs,DNRA freq=freqr hb=h IF(DNRA.LT.0.050D0)DNRA=0.050D0 nnsm=int(delrs/delrm) if(nnsm.lt.1)nnsm=1 delrm=delrs/dfloat(nnsm) if(hb.lt.hhh)hb=hhh nnzz=int((hb+0.1*dh)/dh) numrr=0 ci=(0.0d0,1.0d0) do 100,i=1,ny if(nwhat.le.0)read(28,*)y(i),z(i) if(nwhat.gt.0)read(28,*)z(i),y(i) if(y(i).lt.0.00)go to 400 ya(i)=y(i) za(i)=z(i) 100 continue go to 500 400 m=i-1 500 continue if(hr.ne.0.0)dh=hr if(hr.ne.0.0)nnzz=1 number=i-1 do 750 ll=2,number aslp(ll-1)=(ya(ll)-ya(ll-1))/(za(ll)-za(ll-1)) 750 continue nelmer=number rr=rrr1 do 850 ll=1,number if(h.lt.za(ll))go to 950 yaa(ll)=ya(ll) zaa(ll)=za(ll) if(h.eq.za(ll))go to 877 850 continue go to 860 950 yaa(ll)=aslp(ll-1)*(h-za(ll-1))+ya(ll-1) zaa(ll)=h 860 numbera=ll-1 877 numbera=ll call nmcode(numbera,yaa,zaa) number=nelmer m=number do 1514 j=1,n2 sum1=(0.0d0,0.0d0) do 2114 k=1,n02 sum1=sum1+a2(j,k)*dsin(akz2(k)*hs) 2114 continue c2(j)=sum1/dsqrt(akr2(j)) 1514 continue call nmtll(iib,nnn,0,n2,nnzz,dh,rrr1,rstrt,delrs,rho0, ^rho1,h,b2,g2,atn2,atn2b,c2) if(h.eq.hhh)go to 2468 if(rall.eq.rstrt)go to 2468 h1=h angt=atan((h-hhh)/(rall-rstrt)) atann=dabs((h-hhh)/(rall-rstrt)) rr=rstrt n01=n02 n1=n2 do 171 i=1,n2 c1(i)=c2(i) b1(i)=b2(i) akr1(i)=akr2(i) g1(i)=g2(i) atn1(i)=atn2(i) atn1b(i)=atn2b(i) do 171 j=1,n02 a1(i,j)=a2(i,j) 171 continue do 172 i=1,n02 akz1(i)=akz2(i) sn1(i)=sn2(i) cs1(i)=cs2(i) 172 continue r0=0.0d0 rr2=rstrt lp=0 777 continue lp=lp+1 c12=c00/c11 DHARD=dsqrt(1.0d0-c12*c12)/c00 ana=2.0d0*freq*h*DHARD an=DNRA*100.0d0/ana if(an.lt.2.0d0)an=2.0d0 if(an.gt.20.0d0)an=20.0d0 delr=h/(ana*an*atann) n3=int(delr/delrs) if(n3.lt.1)n3=1 delr=n3*delrs if(delr.lt.delrm)delr=delrm rr=rr+delr if(rr.gt.rall)go to 1002 rr1=rr2 rr2=rr h2=(rall-rr)*tan(angt)+hhh h=h2 do 859 ll=1,number if(h.le.za(ll))go to 959 yaa(ll)=ya(ll) zaa(ll)=za(ll) 859 continue go to 869 959 yaa(ll)=aslp(ll-1)*(h-za(ll-1))+ya(ll-1) zaa(ll)=h 869 numbera=ll call nmcode(numbera,yaa,zaa) number=nelmer m=number nelmer=number 2220 continue if(h2.gt.h1)go to 1023 call makecu ^(n1,n01,n2,n02,r0,rr1,h1,h2,rho0,rho1,akr1,akr2,akz1,akz2, ^atn1,g1,g2,b1,b2,c1,c2,sn1,cs1,sn2,cs2,a1,a2) go to 1024 1023 call makecd ^(n1,n01,n2,n02,r0,rr1,h1,h2,rho0,rho1,akr1,akr2,akz1,akz2, ^atn1,g1,g2,b1,b2,c1,c2,sn1,cs1,sn2,cs2,a1,a2) 1024 continue call nmtll(iib,nnn,lp,n2,nnzz,dh,rr1,rr,delrs,rho0, ^rho1,h,b2,g2,atn2,atn2b,c2) 6701 continue do 405 i=1,n2 b1(i)=b2(i) g1(i)=g2(i) akr1(i)=akr2(i) c1(i)=c2(i) atn1(i)=atn2(i) atn1b(i)=atn2b(i) do 455 j=1,n02 a1(i,j)=a2(i,j) 455 continue 405 continue do 415 i=1,n02 akz1(i)=akz2(i) sn1(i)=sn2(i) cs1(i)=cs2(i) 415 continue n1=n2 n01=n02 h1=h2 r0=rr1 1000 go to 777 1001 continue go to 1003 1002 continue 1003 continue rdif=rall-rr1 h=hhh h2=h if(h1.eq.h2)go to 2021 do 1859 ll=1,number if(h.le.za(ll))go to 1959 yaa(ll)=ya(ll) zaa(ll)=za(ll) 1859 continue go to 1869 1959 yaa(ll)=aslp(ll-1)*(h-za(ll-1))+ya(ll-1) zaa(ll)=h 1869 numbera=ll call nmcode(numbera,yaa,zaa) if(h2.gt.h1)go to 2023 call makecu ^(n1,n01,n2,n02,r0,rall,h1,h2,rho0,rho1,akr1,akr2,akz1,akz2, ^atn1,g1,g2,b1,b2,c1,c2,sn1,cs1,sn2,cs2,a1,a2) go to 2024 2023 call makecd ^(n1,n01,n2,n02,r0,rall,h1,h2,rho0,rho1,akr1,akr2,akz1,akz2, ^atn1,g1,g2,b1,b2,c1,c2,sn1,cs1,sn2,cs2,a1,a2) 2024 continue 2021 continue number=nelmer m=number nelmer=number if(rr2.ge.rall)rr2=rall rend=rend+delrs call nmtll(iib,nnn,2,n2,nnzz,dh,rr2,rend,delrs,rho0, ^rho1,h,b2,g2,atn2,atn2b,c2) 5439 continue 2468 continue close(28) if(iib.lt.1) then write(*,*)'check parameters in"PULSE.IN" ' write(*,*)'nr=',nnn,'nz=',nnzz end if 3333 continue close(24) close(4) close(5) CALL DATIME(DTG) WRITE ( 6, '(1X,''STOP : '',A,/)' ) DTG end subroutine makecu(n1,n01,n2,n02,r0,ri0,h1,h2,rho0,rho1,akr1,akr2, ^akz1,akz2,atn1,g1,g2,b1,b2,c1,c2,sl1,cl1,sm2,cm2,a1,a2) implicit real*8(a-h,o-z) parameter (np=1000) complex*16 c1,c2,ci,sum3,cxy,c1j(np),yyy dimension c1(np),c2(np),akr1(np),akz1(np),b2(np), ^akz2(np),atn1(np),akr2(np) dimension a1(NP,NP),a2(NP,NP),g1(np),g2(np),b1(np) dimension sl2(np),sl1(np),cl2(np),cl1(np),sm2(np), ^cm2(np),f3k(np) dimension aa(NP,NP) eps=0.000000000010d0 rhs=dsqrt(rho0*rho1) dri=ri0-r0 ci=dcmplx(0.00,1.00d0) h21=h2-h1 do 5 l=1,n01 x2=akz1(l)*h2 SL2(l)=dsin(x2) 5 cl2(l)=dcos(x2) do 7 k=1,n2 x=h21*b2(k) if(x.gt.0.0d0)x=-x 7 f3k(k)=dexp(x) do 8 jj=1,n1 X5=AKR1(JJ)*DRI y5=-atn1(jj)*dri yyy=cdexp(ci*x5)*dexp(y5) yya=cdabs(yyy) c1j(jj)=c1(jj)*cdexp(ci*x5)*dexp(y5) 8 continue do 20 j=1,n1 do 20 k=1,n2 sum=0.0d0 sum5=0.0d0 do 30 l=1,n01 akcl1=akz1(l)*cl1(l) akcl2=akz1(l)*cl2(l) akz11=akz1(l)*akz1(l) do 40 m=1,n02 ad=(akz2(m)-akz1(l))*(akz1(l)+akz2(m)) if(dabs(ad).lt.eps)go to 700 akcm2=akz2(m)*cm2(m) f12=akcl2*sm2(m)-akcm2*sl2(l) f12=f12/ad go to 800 700 ah22=2.0d0*akz1(l)*h2 f12=0.50d0*(h2-dsin(ah22)) 800 continue sum=sum+a1(j,l)*a2(k,m)*f12 40 continue f3=f3k(k) f4=b2(k)*sl2(l)+akcl2 f5=b2(k)*sl1(l)+akcl1 f6=b2(k)*b2(k)+akz11 sum=sum*rho0+rhs*g2(k)*a1(j,l)*(f4-f3*f5)/f6 30 continue aa(j,k)=(sum+f3 ^*g1(j)*g2(k)*rho1/(b1(j)+b2(k))) 20 continue do 999 kk=1,n2 sum3=cmplx(0.0d0,0.0d0) sum4=cmplx(0.0d0,0.0d0) do 50 jj=1,n1 y7=dsqrt(akr1(jj)/akr2(kk)) y7=(y7+1.0d0/y7) sum3=sum3+c1j(jj)*aa(jj,kk)*y7 50 continue c2(kk)=sum3*0.50d0 999 continue return end subroutine makecd(n1,n01,n2,n02,r0,ri0,h1,h2,rho0,rho1,akr1,akr2, ^akz1,akz2,atn1,g1,g2,b1,b2,c1,c2,sl1,cl1,sm2,cm2,a1,a2) implicit real*8(a-h,o-z) parameter (np=1000) complex*16 c1,c2,ci,sum3,cxy,c1j(np) real*8 sum dimension c1(np),c2(np),akr1(np),akz1(np), ^akz2(np),atn1(np),akr2(np) dimension a1(NP,NP),a2(NP,NP),g1(np),g2(np),b1(np), ^b2(np) dimension aa(NP,NP) dimension sm2(np),sm1(np),cm2(np),cm1(np),sl1(np), ^cl1(np),f3j(np) eps=0.000000000010d0 rhs=dsqrt(rho0*rho1) dri=ri0-r0 ci=dcmplx(0.00,1.00d0) h12=h1-h2 do 5 m=1,n02 x1=akz2(m)*h1 sm1(m)=dsin(x1) 5 cm1(m)=dcos(x1) do 7 j=1,n1 x=h12*b1(j) if(x.gt.0.0d0)x=-x 7 f3j(j)=dexp(x) do 8 jj=1,n1 X5=AKR1(JJ)*DRI y5=-atn1(jj)*dri cxy=dcmplx(y5,x5) cxy=cdexp(cxy) c1j(jj)=c1(jj)*cxy 8 continue do 20 k=1,n2 do 20 j=1,n1 sum=0.0d0 do 40 m=1,n02 akcm1=akz2(m)*cm1(m) akcm2=akz2(m)*cm2(m) akam=akz2(M)*akz2(M) do 30 l=1,n01 x1=(akz1(l)-akz2(m)) x2=(akz2(m)+akz1(l)) ad=x1*x2 if(dabs(ad).lt.eps)go to 700 acl=akz1(l)*cl1(l) f12=akcm1*sl1(l)-acl*sm1(m) f12=f12/ad go to 800 700 ah22=2.0d0*akz1(l)*h2 f12=0.50d0*(h2-dsin(ah22)) 800 continue sum=sum+a1(j,l)*a2(k,m)*f12 30 continue f4=b1(j)*sm1(m)+akcm1 f5=b1(j)*sm2(m)+akcm2 f6=b1(j)*b1(j)+akam f3=f3j(j) sum=sum*rho0+RHS*g1(j)*a2(k,m)*(f4-f3*f5)/f6 40 continue aa(j,k)=(sum+f3 ^*g1(j)*g2(k)*rho1/(b1(j)+b2(k))) 20 continue do 999 kk=1,n2 sum3=cmplx(0.0d0,0.0d0) do 50 jj=1,n1 y7=dsqrt(akr1(jj)/akr2(kk)) y7=(y7+1.0d0/y7) sum3=sum3+c1j(jj)*aa(jj,kk)*y7 50 continue c2(kk)=sum3*0.50d0 999 continue return end subroutine nmcode(mum,ya,za) implicit real*8 (a-h,o-z) parameter (np=1000,npp=100*np+10,ny=350) common/blkmain1/h,hs,hr,hhh,r3,freq,rho0,rho1,c0,c1,alpha common/ammm/am,akk dimension thl(npp,2),theta(np),akr(np),akz(np) dimension za(ny),ya(ny) number=mum epss=1.d-13 sss=0.0d0 do 77 i=1,mum-1 77 sss=sss+(ya(i+1)+ya(i))*(za(i+1)-za(i)) cave=sss*0.50d0/(za(mum)-za(1)) c0=cave c01a=c0/c1 cc=c0 thcra=dasin(c01a) anmode=2.0d0*freq*h*dcos(thcra)/c0+0.50d0 if(anmode.ge.1.0d0)go to 987 if(anmode.lt.1.0)stop 987 continue 59 continue c01=c0/c1 eps2=0.0050d0 eps=0.000010d0 two=2.0d0 one=1.0d0 pi=4.0d0*datan(1.0d0) rad=pi/180.0d0 ang=180.0d0/pi pi2=two*pi w=pi2*freq continue ak=w/c0 akk=ak akr1=ak*c0/c1 r01=rho0/rho1 c01=c0/c1 c10=c1/c0 thcr=dasin(c01) anmode=2.0d0*freq*h*dcos(thcr)/c0+0.50d0 nmode=anmode nptry=100*nmode+10 c102=c10*c10 th0=dasin(c01) delx1=(pi*0.50d0-th0)/dfloat(nptry) x=(pi+delx1)*0.50d0 do 125 j=1,nptry x=x-delx1 thl(j,1)=x sn=dsin(x) sn2=sn*sn cs=dcos(x) wxyz=c102*sn2-one b1=dsqrt(wxyz) y=datan(c01*r01*b1/cs) thl(j,2)=y 125 continue note=0 imin=1 do 1001 m=1,nmode 99 continue am=dfloat(m) y0=dcos(thl(1,1)) y=h*ak*y0-(am-0.50d0)*pi z=thl(1,2) amin=(z-y)*(z-y) do 225 i=imin+1,npp y0=dcos(thl(i,1)) y=h*ak*y0-(am-0.50d0)*pi z=thl(i,2) zy=(z-y)*(z-y) if(zy.le.amin)imin=i if(zy.le.amin)amin=zy if(zy.le.eps2)go to 56 225 continue 55 try=thl(imin,1)*ang print*,'did not terminate with root will try',imin,try go to 35 56 imin=i 35 continue x3=thl(imin,1) x2=thl(imin-1,1) x1=thl(imin-2,1) z3=g(x3) z2=g(x2) z1=g(x1) do 2 n=1,4000 if(dabs(z3).lt.eps)go to 25 x=f(x1,x2,x3,z1,z2,z3) x1=x2 x2=x3 x3=x z1=z2 z2=z3 z3=g(x) 2 continue 26 print*,' no root at ka=',x,' after',n,' trys','fx=',x3 am=dfloat(m) y0=dcos(thl(1,1)) y=h*ak*y0-(am-0.50d0)*pi z=thl(1,2) amin=(z-y)*(z-y) go to 25 3 format(//9x,'n',24x,'x',26x,'f(x)',//) 4 format(5x,i5,5x,d36.28,5x,e10.3) 25 continue theta(m)=pi*0.50d0-x akr1=ak*dsin(x) akz1=ak*dcos(x) akz(m)=akz1 akr(m)=akr1 1001 continue 1000 continue call Makeq(nmode,number,YA,ZA,akz,akr) 1234 continue return end function f(x1,x2,x3,z1,z2,z3) implicit real*8 (a-h,o-z) x31=x3-x1 x21=x2-x1 z12=z1-z2 z13=z1-z3 z21=z2-z1 z31=z3-z1 zc12=z1*x1-z2*x2 zc13=z1*x1-z3*x3 c=(x21*zc13-x31*zc12)/(x31*z12-x21*z13) b=-(z12*c+zc12)/x21 a=z1*(c+x1)-b*x1 f=-a/b return end function g(x) implicit real*8 (a-h,o-z) common/blkmain1/h,hs,hr,hhh,r3,freq,rho0,rho1,c0,c1,alpha common/ammm/am,ak pi=4.0*datan(1.0d0) one=1.0d0 c01=c0/c1 r01=rho0/rho1 c102=1.0/(c01*c01) cs=dcos(x) ss=dsin(x) z=h*ak*cs-(am-0.50d0)*pi sn2=ss*ss b1=dsqrt(c102*sn2-one) y=datan(c01*r01*b1/cs) g=(z-y)*(z-y) return end subroutine Makeq(nm,number,y1,z1,akz,akr) implicit real*8 (a-h,o-z) parameter (np=1000,ny=350) common/blkmain1/h,hs,hr,hhh,r3,freq,rho0,rho1,c0,c1,alpha common/blksico/si1(np),co1(np) dimension y(ny),akr(np),alm(np,np), ^ y1(ny),q(np,np),vm(np),z1(ny) dimension ab(ny),ama(ny),akz(np) dimension si(np,ny),co(np,ny),ALM2(np,np) two = 2.0d0 pi=4.0*datan(1.0d0) ak=2*pi*freq/c0 m=number cc=c0 crk2=rho0*ak*ak*0.50d0 do 40 i=1,m y(i)=2.0d0*(y1(i)-c0)/c0 40 continue ak1=ak*c0/c1 do 161 k=1,number-1 ama(k)=(y(k+1)-y(k))/(z1(k+1)-z1(k)) ab(k)=y(k)-ama(k)*z1(k) 161 continue do 200 k=1,nm axx=akz(k) do 300 n=1,number xxx=axx*z1(n) si(k,n)=dsin(xxx) 300 co(k,n)=dcos(xxx) 200 continue rh2=rho0*0.50d0 ak2=ak1*ak1 cont=rh2*rho0/rho1 do 787 i=1,nm si1(i)=si(i,number) co1(i)=co(i,number) akr1=akr(i) bm=dsqrt(akr1*akr1-ak2) akz1=akz(i) sn2=si(i,number)*si(i,number) cosi=si(i,number)*co(i,number) alm(i,i)=rh2*(h-cosi/akz1) alm2(i,i)=alm(i,i) vm(i)=cont*sn2/bm+alm(i,i) alm(i,i)=alm(i,i)/vm(i) vm(i)=1.0d0/dsqrt(vm(i)) 787 continue do 50 i=1,nm aki=akz(i) a1=aki*co(i,number) do 50 j=i+1,nm akj=akz(j) b1=akj*co(j,number) abij=(aki-akj)*(aki+akj) alm1=(b1*si(i,number)-a1*si(j,number))/abij alm2(i,j)=alm1 alm2(j,i)=alm1 alm(i,j)=alm1*vm(i)*vm(j) alm(j,i)=alm(i,j) 50 continue 742 continue 741 continue do 750 i=1,nm aki=akz(i) do 760 j=i+1,nm sum=0.0d0 akj=akz(j) call integ(number,i,j,aki,akj,ama,ab,z1,si,co,sum2) 785 q(j,i)=crk2*sum2*(vm(i)*vm(j)) 760 continue 750 continue do 777 i=1,nm aki=akz(i) call integ2(number,i,aki,ama,ab,z1,si,co,sum2) q(i,i)=crk2*sum2*vm(i)*vm(i) 777 continue do 778 i=1,nm do 779 j=i+1,nm q(i,j)=q(j,i) 779 continue 778 continue call MAKEA(nm,ak,vm,akz,akr,alm,alm2,q) return end subroutine integ(number,I,J,aig,ajg,am,ab,z,si,co,sum) implicit real*8 (a-h,o-z) parameter (np=1000,numb=350) dimension z(number),am(number),ab(number) dimension si(np,numb),co(np,numb) sum=0.0d0 dij=aig-ajg dd=1.0d0/dij sij=aig+ajg ss=1.0d0/sij z1=z(1) AA=SI(I,1)*CO(J,1) BB=SI(J,1)*CO(I,1) sns1=(AA+BB) snd1=(AA-BB) AA=SI(I,1)*SI(J,1) BB=CO(J,1)*CO(I,1) cns1=(BB-AA) cnd1=(AA+BB) sd1=(cnd1*dd+z1*snd1)*dd ss1=(cns1*ss+z1*sns1)*ss sds1=snd1*dd sss1=sns1*ss do 10 k=1,number-1 z2=z(k+1) AA=SI(I,k+1)*CO(J,k+1) BB=SI(J,k+1)*CO(I,k+1) sns2=AA+BB snd2=AA-BB AA=SI(I,k+1)*SI(J,k+1) BB=CO(J,k+1)*CO(I,k+1) cns2=BB-AA cnd2=AA+BB sd2=(cnd2*dd+z2*snd2)*dd ss2=(cns2*ss+z2*sns2)*ss 21 sx=am(k)*(sd2-sd1-ss2+ss1) sds2=snd2*dd sss2=sns2*ss sc=ab(k)*(sds2-sds1-sss2+sss1) sum=sum+sc+sx sd1=sd2 ss1=ss2 sds1=sds2 sss1=sss2 10 continue return end subroutine integ2(number,I,aig,am,ab,z,si,co,sum) implicit real*8 (a-h,o-z) parameter (np=1000,numb=350) dimension z(number),am(number),ab(number) dimension si(np,numb),co(np,numb) sum=0.0d0 sij=aig+aig ss=1.0d0/sij ss2=ss*ss z1=z(1) sns1=2.0d0*SI(I,1)*CO(I,1) cns1=2.0d0*CO(I,1)*CO(I,1)-1.0d0 sd1=z1*(z1*0.50d0-sns1*ss)-cns1*ss2 sds1=z1-sns1*ss do 10 k=1,number-1 z2=z(k+1) sns2=2.0d0*SI(I,k+1)*CO(I,k+1) cns2=2.0d0*CO(I,k+1)*CO(I,k+1)-1.0d0 sd2=z2*(z2*0.50d0-sns2*ss)-cns2*ss2 sx=am(k)*(sd2-sd1) sds2=z2-sns2*ss sc=ab(k)*(sds2-sds1) sum=sum+sc+sx sd1=sd2 sds1=sds2 10 continue return end subroutine nmcal(nm,alm) parameter (np=1000) implicit real*8 (a-h,o-z) common/blockncl/nmm,b,g1,atn,atnb common/blockma/nm2,akzz,akrr,a common/blkmain1/h,hs,hr,hhh,r3,freq,rho0,rho1,c0,c1,alpha common/blksico/sn(np),cs(np) dimension a(np,np),akzz(np),b(np),alm(np,np) dimension akrr(np),g1(np),atn(np),atnb(np),avm1(np) nmm=nm nm2=nm pi=4.0d0*Datan(1.0d0) ak=2.0*pi*freq/c0 ak1=ak*c0/c1 cdvid=1.0d0/(c1*c1*20.0*dlog10(dexp(1.0d0))) ckonb=pi*alpha*freq*freq*cdvid ckon=rho1*ckonb cons=(rho0/rho1) do 511 j=1,nm akr1=akrr(j) if(akr1*akr1.le.ak1*ak1)go to 252 b(j)=dsqrt(akr1*akr1-ak1*ak1) 511 continue go to 253 252 CONTINUE nmm=j-1 253 continue 136 continue do 516 j=1,nmm sum2=0.0d0 sum3=0.0d0 do 1116 k=1,nm sum3=sum3+a(j,k)*sn(k) 1116 continue g1(j)=cons*sum3 516 continue rho15=0.50d0*rho1 do 78 i=1,nmm sum0=0.0d0 bm=b(i) do 221 jk=1,nm akzj=akzz(jk) do 221 kk=1,nm s2=0.0d0 akzk=akzz(kk) 331 sum0=sum0+alm(jk,kk)*a(i,jk)*a(i,kk) 221 continue avm1(i)=(rho0*sum0+rho15*g1(i)*g1(i)/bm) avm1(i)=1.0d0/dsqrt(avm1(i)) 78 continue do 1235 k=1,nmm g1(k)=g1(k)*avm1(k) cont=g1(k)*g1(k) atn(k)=ckon*cont/(b(k)*akrr(k)) atnb(k)=ckonb/b(k) 1235 continue do 5321 j=1,nmm do 5321 i=1,nm 5321 a(j,i)=a(j,i)*avm1(j) return end subroutine nmtll(iib,nnn,lp,nmm,nnzz ^ ,dh,rr0,rrf,delrs,rho0,rho1, ^hb,b2,g2,atn,atnb,coef) parameter (np=1000,numz=501,numbr1=51,numbr=20001) implicit real*8 (a-h,o-z) complex*16 sum,ci,cnorm,acmx,coef(np),rrrc character rfl*12 common/blockma/nm,akzz,akrr,a complex*8 pres(numz,numbr) real*4 rl(numbr),zl(numz),tl(numz,numbr),ar,ai,rrr,ttt dimension b2(np),g2(np),a(np,np),atnb(np) dimension r(numbr),akrr(np),akzz(np) ^,atn(np),zr(numz),y(numz,np) ^,psin(np,numz) ci=(0.0,1.0d0) nfln=lp+1 pi=4.0d0*Datan(1.0d0) mmm=nnzz rr00=rr0 if(lp.eq.0)rr00=0.0d0 nnn=0 rri=rr0 999 if(rri.ge.rrf)go to 331 400 nnn=nnn+1 r(nnn)=rri rl(nnn)=sngl(rri) rri=rri+delrs go to 999 331 continue if(rr0.eq.rrf)nnn=1 if(rr0.eq.rrf)r(1)=rr0 if(rr0.eq.rrf)rl(nnn)=sngl(rrf) do 32 i=1,mmm zr(i)=i*dh zl(i)=sngl(zr(i)) 32 continue mmb=mmm do 3001 i=1,mmm if(zr(i).gt.hb)go to 3002 3001 continue go to 3003 3002 mmb=i-1 3003 continue do 2001 k=1,nm do 2002 l=1,mmb akz1=akzz(k)*zr(l) psin(k,l)=dsin(akz1) 2002 continue 2001 continue do 513 j=1,nmm do 513 l=1,mmb sum2=0.0d0 do 1113 k=1,nm sum2=sum2+a(j,k)*psin(k,l) sum2=rho0*sum2 1113 continue 513 y(l,j)=sum2 if(mmb.ge.mmm)go to 1515 do 1514 j=1,nmm sum2=0.0d0 do 1513 l=mmb+1,mmm ccc=(b2(j))*(zr(l)-hb) y(l,j)=rho1*g2(j)*dexp(-ccc) 1513 continue 1514 continue 1515 CNORM=dsqrt(1.0d0*pi) do 40 i=1,nnn do 45 l=1,mmb sum=(0.0,0.0) do 50 j=1,nmm akr=akrr(j) akr1=akr*r(i)+1.0d0 akr=akrr(j) akr1=akr*r(i)+1.0d0 rrs=-atn(j)*(r(i)-rr00) rrrc=ci*(r(i)-rr00)*akr sum=sum+coef(j)*y(l,j)*cdexp(rrrc)*dexp(rrs)/dsqrt(akr1) 50 continue ar=sngl(dreal(cnorm*sum)) ai=sngl(dimag(cnorm*sum)) ttt=10*alog10(ar*ar+ai*ai) rrr=rl(i) tl(l,i)=ttt pres(l,i)=cmplx(ar,ai) 45 continue if(mmb.ge.mmm)go to 1257 do 945 l=mmb+1,mmm sum=(0.0,0.0) do 950 j=1,nmm akr=akrr(j) akr1=akr*r(i)+1.0d0 rrb=-atnb(j)*r(i) rrs=-atn(j)*(r(i)-rr00) rrrc=ci*(r(i)-rr00)*akr sum=sum+coef(j)*y(l,j)*cdexp(rrrc)*dexp(rrs)/dsqrt(akr1) 950 continue ar=sngl(dreal(cnorm*sum)) ai=sngl(dimag(cnorm*sum)) ttt=10*alog10(ar*ar+ai*ai) tl(l,i)=ttt pres(l,i)=cmplx(ar,ai) 945 continue 1257 continue 40 continue do 439 ir=1,nnn do 439 iz=1,mmm if(iib.lt.1) then write(24,*)rl(ir),tl(iz,ir) end if write(UNIT=4)pres(iz,ir) 439 continue if(nnn.lt.100000)go to 2007 nmun=0 write(rfl,1111)'YLM',nfln 1111 format(a,I3.2,a) open(unit=9,status='unknown', ^form='formatted',access='sequential',file=rfl) nmun=nmun+1 do 2005 l=1,mmm write(9,*)l do 2005 m=1,nmm write(9,*)y(l,m) 2005 continue close(9) 2007 continue return end subroutine eigqr(m,lambda,t) parameter(np=1000) implicit real*8 (a-h,o-z) real*8 t(np,np),lambda(np),e(np) integer la(np) call tred2(t,np,m,lambda,e) call tqli(lambda,e,np,m,T) j=0 do 100 i=1,m amax=0.0d0 j=j+1 do 200 k=j,m if(lambda(k).lt.amax)go to 200 amax=lambda(k) la(i)=k 200 continue x=lambda(i) lambda(i)=amax lambda(la(i))=x do 400 n=1,m x=t(n,i) t(n,i)=t(n,la(i)) 400 t(n,la(i))=x 100 continue 300 continue return end subroutine tqli(d,e,np,n,z) implicit real*8 (a-h,o-z) dimension d(np),e(np),z(np,np) if (n.gt.1) then do 11 i=2,n e(i-1)=e(i) 11 continue e(n)=0. do 15 l=1,n iter=0 1 do 12 m=l,n-1 dd=dabs(d(m))+dabs(d(m+1)) if (dabs(e(m))+dd.eq.dd) go to 2 12 continue m=n 2 if(m.ne.l)then iter=iter+1 g=(d(l+1)-d(l))/(2.0d0*e(l)) r=dsqrt(g*g+1.0d0) g=d(m)-d(l)+e(l)/(g+sign(r,g)) s=1.0d0 c=1.0d0 p=0.0d0 do 14 i=m-1,l,-1 f=s*e(i) b=c*e(i) if(dabs(f).ge.dabs(g))then c=g/f r=sqrt(c*c+1.0d0) e(i+1)=f*r s=1.0d0/r c=c*s else s=f/g r=sqrt(s*s+1.0d0) e(i+1)=g*r c=1.0d0/r s=s*c endif g=d(i+1)-p r=(d(i)-g)*s+2.0d0*c*b p=s*r d(i+1)=g+p g=c*r-b do 13 k=1,n f=z(k,i+1) z(k,i+1)=s*z(k,i)+c*f z(k,i)=c*z(k,i)-s*f 13 continue 14 continue d(l)=d(l)-p e(l)=g e(m)=0.0d0 go to 1 endif 15 continue endif return end SUBROUTINE tred2(a,np,n,d,e) IMPLICIT REAL*8 (a-h,o-z) DIMENSION a(np,np),d(np),e(np) do 18 i=n,2,-1 l=i-1 h=0. scale=0.0d0 if(l.gt.1)then do 11 k=1,l scale=scale+dabs(a(i,k)) 11 continue if(scale.eq.0.0d0)then e(i)=a(i,l) else do 12 k=1,l a(i,k)=a(i,k)/scale h=h+a(i,k)**2 12 continue f=a(i,l) g=-sign(dsqrt(h),f) e(i)=scale*g h=h-f*g hd=1.0d0/h a(i,l)=f-g f=0.0d0 do 15 j=1,l a(j,i)=a(i,j)/h g=0.0d0 do 13 k=1,j g=g+a(j,k)*a(i,k) 13 continue do 14 k=j+1,l g=g+a(k,j)*a(i,k) 14 continue e(j)=g/h f=f+e(j)*a(i,j) 15 continue hh=f/(h+h) do 17 j=1,l f=a(i,j) g=e(j)-hh*f e(j)=g do 16 k=1,j a(j,k)=a(j,k)-f*e(k)-g*a(i,k) 16 continue 17 continue endif else e(i)=a(i,l) endif d(i)=h 18 continue d(1)=0.0d0 e(1)=0.0d0 do 24 i=1,n l=i-1 if(d(i).ne.0.0d0)then do 22 j=1,l g=0. do 19 k=1,l g=g+a(i,k)*a(k,j) 19 continue do 21 k=1,l a(k,j)=a(k,j)-g*a(k,i) 21 continue 22 continue endif d(i)=a(i,i) a(i,i)=1. do 23 j=1,l a(i,j)=0.0d0 a(j,i)=0.0d0 23 continue 24 continue return END SUBROUTINE chol(a,n,np) IMPLICIT REAL*8 (a-h,o-z) DIMENSION a(np,np) a(1,1)=dsqrt(a(1,1)) aa=1.0d0/a(1,1) do 10 j=2,n a(1,j)=a(1,j)*aa 10 continue do 20 i=2,n sum=a(i,i) do 30 k=1,i-1 30 sum=sum-a(k,i)*a(k,i) a(i,i)=dsqrt(sum) aa=1.0d0/a(i,i) do 40 j=i+1,n sum=a(i,j) do 50 k=1,i-1 sum=sum-a(k,i)*a(k,j) 50 continue a(i,j)=sum*aa 40 continue 20 continue return END subroutine MAKEA(n,ak,vm,akz,akr,ALM,alm2,q) implicit real*8 (a-h,o-z) parameter (np=1000) common/blockma/nm,akzz,akrr,a dimension am(np,np),vm(np),q(np,np),ALM(NP,NP), ^ a(np,np),eig(np),akr(np),akz(np),akzz(np),akrr(np),ALM2(NP,NP) nm=n call chol(alm,n,np) call GETY(q,alm,n,np) call GETC(q,alm,n,np) do 179 i=1,n ak1=akr(i) ak2=ak1*ak1 a(i,i)=AK2-q(i,i) do 79 j=i+1,n a(i,j)=-q(i,j) 79 a(j,i)=a(i,j) 179 continue call eigqr(n,eig,a) do 123 i=1,n do 125 j=1,n 125 am(i,j)=a(j,i) 123 continue call GETX(am,alm,n,np) do 75 ir =1,n d=am(ir,ir) dd=1.0d0/d do 70 ic=1,n a(ir,ic)=Am(ir,ic)*vm(ic)*dd 70 continue 75 CONTINUE do 2000 i=1,n akzz(i)=akz(i) akrr(i)=dsqrt(eig(i)) 2000 continue call nmcal(n,alm2) return end SUBROUTINE GETY(a,al,n,np) IMPLICIT REAL*8 (a-h,o-z) DIMENSION a(np,np),al(np,np) aa=1.0d0/al(1,1) do 10 jc=1,n a(1,jc)=a(1,jc)*aa 10 continue do 20 ir=2,n aa=1.0d0/al(ir,ir) do 30 jc=1,ir sum=a(ir,jc) do 40 k=1,ir-1 40 sum=sum-al(ir,k)*a(k,jc) 30 a(ir,jc)=sum*aa 20 continue do 50 ir=1,n do 50 jc=1,ir 50 a(jc,ir)=a(ir,jc) return end SUBROUTINE GETX(a,al,n,np) IMPLICIT REAL*8 (a-h,o-z) DIMENSION a(np,np),al(np,np) aa=1.0d0/al(n,n) do 10 ir=1,n a(ir,n)=a(ir,n)*aa 10 continue do 20 jc=n-1,1 aa=1.0d0/al(jc,jc) do 30 ir=1,n sum=a(ir,jc) do 40 k=jc+1,n 40 sum=sum-al(ir,k)*a(k,jc) 30 a(ir,jc)=sum*aa 20 continue return end SUBROUTINE GETC(a,al,n,np) IMPLICIT REAL*8 (a-h,o-z) DIMENSION a(np,np),al(np,np) aa=1.0d0/al(n,n) do 10 i=1,n a(i,n)=a(i,n)*aa 10 continue do 20 j=n-1,1,-1 aa=1.0d0/al(j,j) do 30 i=1,j sum=a(i,j) do 40 k=j+1,n 40 sum=sum-a(i,k)*al(j,k) 30 a(i,j)=sum*aa 20 continue do 50 j=n,1,-1 do 50 i=1,j 50 a(j,i)=a(i,j) return end subroutine datime(dtg) character dtg*(*) cc CHANGING THE NAME OF VARIABLE DATE TO ADATE SINCE DATE IS A FUNCTION character adate*11, months(12)*3, time*11 integer*2 day,hour,minute,second,i100th,month,year data months / 'Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', &'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec' / cc ADDED BY UTSAV/KILARI FOR SOLARIS integer date_time(8) character*10 b(3) character newdate*11 character newtime*8 integer newmonth cc call date (newdate) cc call gettim (hour,minute,second,i100th) cc write(date,'(I2,1X,A3,1X,I4)') day,months(month),year call date_and_time(b(1),b(2),b(3),date_time) newmonth = date_time(2) write(newdate,'(I2,1X,A3,1X,I4)') date_time(3),months(newmonth),date_time(1) cc added by UTSAV/KILARI for debugging cc write (*,'(I2,1X,A3,1X,I2)') newdate cc write(newtime,'(2(I2,1H:),I2)') date_time(5),date_time(6),date_time(7) cc write(time,'(2(I2,1H:),I2)') hour,minute,second cc write(dtg,'(A11,1X,A8)') date,time write(dtg,'(A11,1X,A8)') newdate,newtime return end cc subroutine datime(dtg) cc character dtg*(*) cc character date*11, months(12)*3, time*11 cc integer*2 day,hour,minute,second,i100th,month,year cc data months / 'Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', cc &'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec' / cc call getdat (year,month,day) cc call gettim (hour,minute,second,i100th) cc write(date,'(I2,1X,A3,1X,I4)') day,months(month),year cc write(time,'(2(I2,1H:),I2)') hour,minute,second cc write(dtg,'(A11,1X,A8)') date,time cc return cc end