C Program residence C Alan Burns 2007 parameter (imx=73,jmx=36,ijtt=20,nmx=5,imxint=imx*nmx) parameter (jmxint=jmx*nmx,nut=210,imx2=(imx-1)*nmx+1) parameter (jmx2=(jmx-1)*nmx+1,jmx3=jmx2+2,mut=121) parameter (ntj=4) C Declare variables C electron density real ne(imx,jmx,nut) real ne2(imx,jmx) real ed(imx,jmx) C ion velocity real uiz(imx,jmx,nut) real viz(imx,jmx,nut) real ui(imx,jmx) real vi(imx,jmx) C Frictional heating real htjout(imx,jmx,nut) real htjouz1(imx,jmx,nut) real htjoum(imx,jmx,nut) real htjouti(imx,jmx) real htjouz1i(imx,jmx) real htjoumi(imx,jmx) C latitude line of number fluxes real alinend(jmx3) C do we need to interpret time real workdata(imxint,jmxint,ijtt) real holdwork(imx2,jmx3) real uint(imx2,jmx3) real vint(imx2,jmx3) real edint(imx2,jmx3) real edout(imx2,jmx3,mut) real uiout(imx2,jmx3,mut) real viout(imx2,jmx3,mut) common /velsin/uiout,viout real htjoutint(imx2,jmx3) real htjouz1int(imx2,jmx3) real htjoumint(imx2,jmx3) real htjoutout(imx2,jmx3,mut) real htjouz1out(imx2,jmx3,mut) real htjoumout(imx2,jmx3,mut) common /htin/htjoutout,htjouz1out,htjoumout real alat(mut),alon(mut) real alata(mut),alona(mut) real xlt(mut) real xltg(mut) common /trajout1/alon,alat,alona,alata,xlt common /trajout2/xltg real htjoua(mut) real htjoub(mut) real htjouc(mut) real alat1(ntj,mut),alon1(ntj,mut) real alata1(ntj,mut),alona1(ntj,mut) real xlt1(ntj,mut) real outbound(24) real inbound(24) real outbound2(24,mut) real inbound2(24,mut) real outboundmag(24,2) real inboundmag(24,2) real alonn(40,24),alatn(40,24) real time(40,24) common /timout/time,alonn,alatn real hintf(40,24) real hintp(40,24) real hintm(40,24) common /htout/hintf,hintp,hintm C Need location of magnetic poles in geographic coordinates C COORDS OF S (LAT,LON) AND N OR CENTERED (LAT,LON) DIPOLE real mag(4) DATA MAG/-74.5,127.,79.,-70./ real work(4) C Set up common blocks pi = atan(1.) C put lat lon magnetic north pole into variables C Change this to radians DTOR=ATAN(1.)/45. rlatmp = mag(3)*ATAN(1.)/45. rlonmp = mag(4)*dtor C Read in characteristic energy, flux(or just electron density) and ion C drifts in geographic coordinates C Keep 20 minutes of information, use geographic coordinates, so corotation in C Read in electron densities open(11,file='input/chren.dat') read(11,5000)ne 5000 format(10E11.4) close(11) C Read in zonal ion drifts in geographic coordinates open(11,file='input/uig.dat') read(11,5000)uiz close(11) C Read in meridional ion drifts in geographic coordinates open(11,file='input/vig.dat') read(11,5000)viz close(11) C Read in frictional heating data 1. Total height integrated per unit C volume open(11,file='input/htintjoul.dat') read(11,5000)htjout close(11) C write(*,*)htjout C Read in frictional heating data 2. Above z=-1 height integrated per unit C volume open(11,file='input/htintjoul2.dat') read(11,5000)htjouz1 close(11) C Read in frictional heating data 3. Total height integrated per unit C mass open(11,file='input/htintjoulm.dat') read(11,5000)htjoum close(11) write(*,*)'finished reading' C Use linear interpolation to fill in grid (use nxn, where n=nmx, needs to be C set in subroutine as well.) jut=102 do i=1,imx,1 do j=1,jmx,1 ne2(i,j)=ne(i,j,jut) enddo enddo ipol=1 call interp(ne2,holdwork,ipol) C write(13,*)holdwork C Loop through local times (magnetic or geographic? do i=1,24,1 ii = 12 ix = (i-1)*3*5+1 write(*,*)'ix',ix C ix=10 write(*,*)'ix',ix do j = 1,jmx3,1 C Try this without the reverse alinend(j) = holdwork(ix,j) C if(holdwork(ix,j) .eq. 0.)write(*,*)'zero',i,j enddo C call gtm(glat,glon,rlatmp,rlonmp,rlatm,rlonm,dip,dec,work,1) write(*,*)'finished gtm' C Find the position of the auroral oval for every hour of local time in C geomagnetic coordinates or geographic coordinates call searchoval(alinend,ifrst,jfrst,ilast,jlast) outbound(i) = jfrst - 87.5 inbound(i) = jlast - 87.5 if(inbound(i) .gt. 85.)inbound(i)=85. C magnetic stuff that was not used eventually gnlon = (i-1)*15.*dtor gnlat = outbound(i)*dtor C write(*,*)'gnlon gnlat',gnlon,gnlat call gtm(gnlat,gnlon,rlatmp,rlonmp,rlatm,rlonm,dip,dec,work,1) outboundmag(i,1) = rlatm/dtor outboundmag(i,2) = rlonm/dtor gnlat = inbound(i)*dtor call gtm(gnlat,gnlon,rlatmp,rlonmp,rlatm,rlonm,dip,dec,work,1) inboundmag(i,1) = rlatm/dtor inboundmag(i,2) = rlonm/dtor C write(*,*)outboundmag(i,1),outboundmag(i,2) C Run back trajectories until parcel outside oval find residence time C- perhaps 20 minutes + is C better do this in a subroutine C Use interpolation to find characteristic energy + number flux along tajectory C Joule heating Te too? enddo C Write out boundaries in m agnetic coordinates open(11,file='outbound.dat') write(11,5000)outbound close(11) open(11,file='inbound.dat') write(11,5000)inbound close(11) C interpolate ui and vi over 2 hours of time in one minute intervals do kut = 1,7,1 do i=1,imx,1 do j=1,jmx,1 ui(i,j)=uiz(i,j,jut-kut+1) enddo enddo do i=1,imx,1 do j=1,jmx,1 vi(i,j)=viz(i,j,jut-kut+1) enddo enddo if(kut .eq. 1)then open(12,file='uiout.dat') write(12,5000)ui close(12) open(12,file='viout.dat') write(12,5000)vi close(12) endif ipol=-1 call interp(ui,uint,ipol) call interp(vi,vint,ipol) do i=1,imx2,1 do j=1,jmx3,1 uiout(i,j,(kut-1)*20+1) = uint(i,j) viout(i,j,(kut-1)*20+1) = vint(i,j) enddo enddo enddo C interpolate htjout etc over 2 hours of time in one minute intervals C + spatially do kut = 1,7,1 do i=1,imx,1 do j=1,jmx,1 htjouti(i,j)=htjout(i,j,jut-kut+1) enddo enddo ipol=-1 call interp(htjouti,htjoutint,ipol) do i=1,imx2,1 do j=1,jmx3,1 htjoutout(i,j,(kut-1)*20+1) = htjoutint(i,j) enddo enddo enddo do kut = 1,7,1 do i=1,imx,1 do j=1,jmx,1 htjouz1i(i,j)=htjouz1(i,j,jut-kut+1) enddo enddo ipol=-1 call interp(htjouz1i,htjouz1int,ipol) do i=1,imx2,1 do j=1,jmx3,1 htjouz1out(i,j,(kut-1)*20+1) = htjouz1int(i,j) enddo enddo enddo do kut = 1,7,1 do i=1,imx,1 do j=1,jmx,1 htjoumi(i,j)=htjoum(i,j,jut-kut+1) enddo enddo ipol=-1 call interp(htjoumi,htjoumint,ipol) do i=1,imx2,1 do j=1,jmx3,1 htjoumout(i,j,(kut-1)*20+1) = htjoumint(i,j) enddo enddo enddo C open(18,file='outbound.dat') C write(18,5000)outbound C close(18) C open(18,file='inbound.dat') C write(18,5000)inbound C close(18) C interpolate uint and vint over 2 hours of time in one minute intervals call intertime(uiout) call intertime(viout) C interpolate htjoutint over 2 hours of time in one minute intervals call intertime(htjoutout) call intertime(htjouz1out) call intertime(htjoumout) C open(13,file='uiout.dat') C write(13,5000)uiout C close(13) C open(13,file='viout.dat') C write(13,5000)viout c close(13) C write(*,*)htjoutout utstrt = 1200. alons=75. C this is 10 versus 10 given that UT is 12 if(alons .gt. 180.)alons=alons-360. alats = inbound(17) -15 call traj(alons,alats,utstrt) do jut=1,121,1 alat1(1,jut) = alat(jut) alon1(1,jut) = alon(jut) xlt1(1,jut) = xlt(jut) enddo C Calculate frictional heating along the trajectory C alata and alona are in geographic cordinates so can be used C to get i and j for height integrated Joule heating. do jut=1,121,1 i = alona(jut) + 181. j = alata(jut) + 88.5 htjoua(jut) = htjoutout(i,j,jut) htjoub(jut) = htjouz1out(i,j,jut) htjouc(jut) = htjoumout(i,j,jut) enddo open(18,file='traj1lat.dat') write(18,5000)alata close(18) open(18,file='traj1lon.dat') write(18,5000)alona close(18) open(18,file='traj1xlt.dat') write(18,5000)xltg close(18) open(18,file='traj1htja.dat') write(18,5000)htjoua close(18) open(18,file='traj1htjb.dat') write(18,5000)htjoub close(18) open(18,file='traj1htjc.dat') write(18,5000)htjouc close(18) write(*,*)alata(1),alata(121) alons=-3. C this is 10 versus 10 given that UT is 12 if(alons .gt. 180.)alons=alons-360. alats = inbound(12) -5. call traj(alons,alats,utstrt) do jut=1,121,1 alat1(1,jut) = alat(jut) alon1(1,jut) = alon(jut) xlt1(1,jut) = xlt(jut) enddo do jut=1,121,1 i = alona(jut) + 181. j = alata(jut) + 88.5 htjoua(jut) = htjoutout(i,j,jut) htjoub(jut) = htjouz1out(i,j,jut) htjouc(jut) = htjoumout(i,j,jut) enddo open(18,file='traj2lat.dat') write(18,5000)alata close(18) open(18,file='traj2lon.dat') write(18,5000)alona close(18) open(18,file='traj2xlt.dat') write(18,5000)xltg close(18) open(18,file='traj2htja.dat') write(18,5000)htjoua close(18) open(18,file='traj2htjb.dat') write(18,5000)htjoub close(18) open(18,file='traj2htjc.dat') write(18,5000)htjouc close(18) alons=-90. C this is 10 versus 10 given that UT is 12 if(alons .gt. 180.)alons=alons-360. alats = inbound(6) -15. call traj(alons,alats,utstrt) do jut=1,121,1 alat1(1,jut) = alat(jut) alon1(1,jut) = alon(jut) xlt1(1,jut) = xlt(jut) enddo do jut=1,121,1 i = alona(jut) + 181. j = alata(jut) + 88.5 htjoua(jut) = htjoutout(i,j,jut) htjoub(jut) = htjouz1out(i,j,jut) htjouc(jut) = htjoumout(i,j,jut) enddo open(18,file='traj3lat.dat') write(18,5000)alata close(18) open(18,file='traj3lon.dat') write(18,5000)alona close(18) open(18,file='traj3xlt.dat') write(18,5000)xltg close(18) open(18,file='traj3htja.dat') write(18,5000)htjoua close(18) open(18,file='traj3htjb.dat') write(18,5000)htjoub close(18) open(18,file='traj3htjc.dat') write(18,5000)htjouc close(18) alons=-180. C this is 10 versus 10 given that UT is 12 if(alons .gt. 180.)alons=alons-360. alats = inbound(1) -5. call traj(alons,alats,utstrt) do jut=1,121,1 alat1(1,jut) = alat(jut) alon1(1,jut) = alon(jut) xlt1(1,jut) = xlt(jut) enddo do jut=1,121,1 i = alona(jut) + 181. j = alata(jut) + 88.5 htjoua(jut) = htjoutout(i,j,jut) htjoub(jut) = htjouz1out(i,j,jut) htjouc(jut) = htjoumout(i,j,jut) enddo open(18,file='traj4lat.dat') write(18,5000)alata close(18) open(18,file='traj4lon.dat') write(18,5000)alona close(18) open(18,file='traj4xlt.dat') write(18,5000)xltg close(18) open(18,file='traj4htja.dat') write(18,5000)htjoua close(18) open(18,file='traj4htjb.dat') write(18,5000)htjoub close(18) open(18,file='traj4htjc.dat') write(18,5000)htjouc close(18) do kut = 1,7,1 do i=1,imx,1 do j=1,jmx,1 ed(i,j)=ne(i,j,jut-kut+1) enddo enddo ipol=1 call interp(ed,edint,ipol) do i=1,imx2,1 do j=1,jmx3,1 edout(i,j,(kut-1)*20+1) = edint(i,j) enddo enddo enddo call intertime(edout) do jut = 1,121,1 do i=1,imx2,1 do j=1,jmx3,1 holdwork(i,j) = edout(i,j,jut) enddo enddo do i=1,24,1 ix = (i-1)*3*5+1 do j = 1,jmx3,1 alinend(j) = holdwork(ix,j) enddo C Find the position of the auroral oval for every hour of local time in C geomagnetic coordinates or geographic coordinates call searchoval(alinend,ifrst,jfrst,ilast,jlast) outbound2(i,jut) = jfrst - 87.5 inbound2(i,jut) = jlast - 87.5 if(inbound2(i,jut) .gt. 85.)inbound2(i,jut)=85. enddo enddo write(*,*)inbound(12),inbound2(12,1), 1outbound(12),outbound2(12,1) C work out time until trajectory is out of oval call reside(inbound,outbound) C C output residence times and positions of end points C output trajectory positions and values along each trajectory C end open(12,file='alatn.dat') write(12,5000)alatn close(12) open(12,file='alonn.dat') write(12,5000)alonn close(12) open(12,file='time.dat') write(12,5000)time close(12) open(12,file='hintf.dat') write(12,5000)hintf close(12) open(12,file='hintp.dat') write(12,5000)hintp close(12) open(12,file='hintm.dat') write(12,5000)hintm close(12) stop end C C ************************************************************************ C ************************************************************************ C ********** Subroutine interp ********** C ************************************************************************ C ************************************************************************ C C Routine to interpolate in space subroutine interp(input,work,ipol) parameter (imx=73,jmx=36,ijtt=20,nmx=5,imxint=imx*nmx) parameter (jmxint=jmx*nmx,nut=210,imx2=(imx-1)*nmx+1) parameter (jmx2=(jmx-1)*nmx+1,jmx3=jmx2+2) real input(imx,jmx),work(imx2,jmx3) real hold(imxint,jmxint) do j = 1,jmx,1 do i = 1,imx,1 work((i-1)*nmx+1,(j-1)*nmx+1) = input(i,j) C if(input(i,j) .eq. 0.)write(*,*)'input 0',i,j enddo enddo C write(*,*)'jmx2',jmx2,imx2 C write(*,*)input(imx,jmx),work((imx-1)*nmx+1,(jmx-1)*nmx+1) do j = 1,jmx,1 do i = 1,imx-1,1 do n = 1,nmx,1 ii = (i-1)*nmx + n wt1 = (float(nmx)-n+1)/float(nmx) wt2 = (n-1)/float(nmx) work(ii,(j-1)*nmx+1) = 1 work((i-1)*nmx+1,(j-1)*nmx+1)*wt1+ 2 work(i*nmx+1,(j-1)*nmx+1)*wt2 C if(work((i)*nmx+1,(j-1)*nmx+1) .eq. 0.) C 1 write(*,*)'work is zero',i,j C write(*,*)wt1,wt2,ii,i,j,hold(ii,j),input(i,j) enddo enddo enddo C write(14,*)hold do i = 1,imx2,1 do j = 1,jmx-1,1 do n = 1,nmx,1 jj = (j-1)*nmx + n wt1 = (float(nmx)-n+1)/float(nmx) wt2 = (n-1)/float(nmx) work(i,jj) = 1 work(i,(j-1)*nmx+1)*wt1+ 2 work(i,j*nmx+1)*wt2 C write(*,*)wt1,wt2,hold(i,j),work(i,jj) enddo enddo enddo c fill in northern points up to 89.5 C Note that I have not yet done this for the southern hemisphere C Note also that I am hardwiring the resolution at present. AGB do i = 1,imx2,1 do j = jmx2+1,jmx2+2,1 i2 = i + imx2/2 if(i2 .gt. imx2)i2 = i2 - imx2 C write(*,*)'i2 i',i2,i wt1 = (float(nmx)-j-jmx2)/float(nmx) wt2 = (j-jmx2)/float(nmx) work(i,j) = wt1*work(i,jmx2) + ipol*wt2*work(i2,jmx2) enddo enddo C write(*,*)'out of interp' C write(*,*)work(10,(jmx-1)*nmx+1) return end C C ************************************************************************ C ************************************************************************ C ********** Subroutine intertime ********** C ************************************************************************ C ************************************************************************ C subroutine intertime(vout) parameter (imx=73,jmx=36,ijtt=20,nmx=5) parameter (imx2=(imx-1)*nmx+1,mut=121) parameter (jmx2=(jmx-1)*nmx+1,jmx3=jmx2+2) real vout(imx2,jmx3,mut) do kut=1,6,1 do lut = 1,19,1 do i = 1,imx2,1 do j = 1,jmx3,1 vout(i,j,lut+1+(kut-1)*20) = vout(i,j,(kut-1)*20+1)* 1 (20.-float(lut))/20. + vout(i,j,kut*20+1)*float(lut)/20. enddo enddo enddo enddo return end C C ************************************************************************ C ************************************************************************ C ********** Subroutine searchoval ********** C ************************************************************************ C ************************************************************************ C subroutine searchoval(alinend,ifrst,jfrst,ilast,jlast) parameter (imx=73,jmx=36,ijtt=20,nmx=5,imxint=imx*nmx) parameter (jmx2=(jmx-1)*nmx+1,jmx3=jmx2+2) real alinend(jmx3) C Find the maximum precipitation along a line of latitude ig =0 ival =90 ign=0 ign2=0 ign3=0 pmax = 0. do i = jmx2-90,jmx3 if(alinend(i) .gt. pmax)pmax=alinend(i) enddo do i = jmx2-90,jmx3 if(ign .ne. 1)then if(alinend(i) .gt. 0.5)then ign=1 jfrst= i endif endif if(ign .eq. 1)then if(alinend(i+1) .lt. alinend(i))then ign2=1 endif endif if((ign2 .eq. 1) .and. (ign3 .ne. 1))then if((((alinend(i+1) .gt. alinend(i))) .and. 1 (i .gt. (jfrst+10))) 1 .or. (alinend(i) .lt. pmax/3.))then jlast=i ign3=1 endif endif enddo write(*,*)jfrst,jlast C Search over all latitudes from the pole for end of auroral oval C Assume maxval of number flux / div is a reasonable boundary return end C C ************************************************************************ C ************************************************************************ C ********** Subroutine reside ********** C ************************************************************************ C ************************************************************************ C subroutine reside(inbound,outbound) parameter (imx=73,jmx=36,ijtt=20,nmx=5,imxint=imx*nmx) parameter (jmxint=jmx*nmx,nut=210,imx2=(imx-1)*nmx+1) parameter (jmx2=(jmx-1)*nmx+1,jmx3=jmx2+2,mut=121) C real inbound(24,mut),outbound(24,mut) real inbound(24),outbound(24) real alonn(40,24),alatn(40,24) real time(40,24) real uiout(imx2,jmx3,mut) real viout(imx2,jmx3,mut) common /velsin/uiout,viout common /timout/time,alonn,alatn real htjoutout(imx2,jmx3,mut) real htjouz1out(imx2,jmx3,mut) real htjoumout(imx2,jmx3,mut) common /htin/htjoutout,htjouz1out,htjoumout real hintf(40,24) real hintp(40,24) real hintm(40,24) common /htout/hintf,hintp,hintm real alat(mut),alon(mut) real alata(mut),alona(mut) real xlt(mut) common /trajout1/alon,alat,alona,alata,xlt utstrt = 1200. C Set up arrayof start points C do iwid = 1,20,1 do lt=1,24,1 do iwid = 1,40,1 alonn(iwid,lt) = (lt-1)*15. -180. alatn(iwid,lt) = inbound(lt)-iwid+1 time(iwid,lt) = 0. c if((inbound(lt,1)-iwid+1) .gt. outbound(lt,1))then if((inbound(lt)-iwid+1) .gt. outbound(lt))then call traj(alonn(iwid,lt),alatn(iwid,lt),utstrt) C Find out how long trajectory is in oval istop = 0 joulct1 = 0. joulct2 = 0. joulct3 = 0. do jut = 1,121,1 ilt = ifix(alona(jut)/15.)+13 if(ilt .le. 0)ilt=ilt+24 if(ilt .gt. 24)ilt=ilt-24 iii = (alona(jut) + 180.) + 1 jjj = (alata(jut) + 87.5) + 1 if(istop .eq. 0)then C write(*,*)htjoutout joulct1 = joulct1 + htjoutout(iii,jjj,jut) joulct2 = joulct2 + htjouz1out(iii,jjj,jut) joulct3 = joulct3 + htjoumout(iii,jjj,jut) endif if(istop .eq. 1)then hintf(iwid,lt) = joulct1 hintp(iwid,lt) = joulct2 hintm(iwid,lt) = joulct3 endif C if((alata(jut) .gt. inbound(ilt,jut)) .and. (istop .eq. 0))then if((alata(jut) .gt. inbound(ilt)) .and. (istop .eq. 0))then if(jut .eq. 1)then write(*,*)alata(jut),inbound(ilt),ilt,lt,alatn(iwid,lt), 1 alonn(iwid,lt),alona(jut),'inbound' endif time(iwid,lt) = jut-1 istop = 1 endif C if((alata(jut) .lt. outbound(ilt,jut)) .and. (istop .eq. 0))then if((alata(jut) .lt. outbound(ilt)) .and. (istop .eq. 0))then if(jut .eq. 1)then write(*,*)alata(jut),outbound(ilt),ilt,lt,alatn(iwid,lt), 1 alonn(iwid,lt),alona(jut),'outbound' endif time(iwid,lt) = jut-1 istop = 1 endif if((jut .eq. 121) .and. (istop .eq. 0))then time(iwid,lt) = jut istop = 1 endif enddo C write(*,*)lt,time(iwid,lt),alatn(iwid,lt),alonn(iwid,lt) C write(*,*)inbound(lt),outbound(lt),iwid C write(*,*)lt,ilt,iwid,istop,time(iwid,lt) endif enddo C write(*,*)lt,time(iwid,lt),alatn(iwid,lt),alonn(iwid,lt)yy enddo return end