program mix use MixReader ! Run: gmake; then run: driver real, allocatable :: gridX(:,:),gridY(:,:) real, allocatable :: potN(:,:),pots(:,:) real, allocatable :: facN(:,:),facs(:,:) real, allocatable :: pedN(:,:),peds(:,:) real, allocatable :: hallN(:,:),halls(:,:) real, allocatable :: engN(:,:),engs(:,:) real, allocatable :: fluxN(:,:),fluxs(:,:) real, allocatable :: nwN(:,:),nws(:,:) integer, parameter :: ndata=365*12*24, iyr=2006 ! integer, parameter :: ndata=365*12*24, iyr=2007 ! integer, parameter :: ndata=366*12*24, iyr=2008 character*31 mixdir06 ! for Aug 2013 output from Dec 13-16, 2006 ! character*53 mixdir06 ! for TIEGCM 1.94.2 version to 0 UT day 350 ! character*47 mixdir06 ! for TIEGCM 1.94.1 version to ~5 UT day 349 character*38 filtim06 ! character*85 filenam06 ! for 1.94.1 version character*91 filenam06 ! for 1.94.2 version character*20 datehrmn real :: radius(27,181),theta(181) integer nji(2),njlat,nimlt,j,i,nmph,dmin,n23 real :: mlatdeg(27),colatrad(27),mltn(181),mlts(181),mlt(181),colatdeg(27) real :: pii,dtr,htr real :: latnx,latnn,latsx,latsn,latnamx,latnpmn,latsamx,latspmn real :: mltnx,mltnn,mltsx,mltsn,mltnamx,mltnpmn,mltsamx,mltspmn real :: theta0(2),offa(2),dskofa(2),phid(2),phin(2),rrad(2),offc(2),dskofc(2),model_ctpoten(2),byimf real :: date(ndata),bx(ndata),by(ndata),bz(ndata),swvel(ndata),swden(ndata) real :: epotSH(181,27),epotNH(181,27),efluxSH(181,27),efluxNH(181,27) real :: offcen(2),dskofcen(2),radcen(2) integer :: inn,inx,jnn,jnx,isn,isx,jsn,jsx real :: mltday,mltnit ! 09/12 bae: eflux=flux*fac_p2e*alfa fac_p2e=2*1.602e-9 for alfa, so half for eng=2*alfa real, parameter:: fac_p2e = 1.602e-9 ! convert from particle to energy flux real :: Re,areatot,dmlon,hpi_sh_mix,hpi_nh_mix real :: da(181),mlat2pol(27) real :: mlat_sh_max,mlat_nh_max,hpi_sh_max,hpi_nh_max pii = 4.*atan(1.) dtr = pii/180. rtd = 1./dtr htr = pii/12. ! nh0 = 0 usually and nm1 = 1 usually except for ifwr=1 for 1 time only nh0 = 0 nm1 = 1 ! nmph = number of minutes per hour (30x2min WHI, 12x5min AGU Dec06) ! dmin = number of minutes between CMIT files ! ifwr = 1,0 if do,not write out epot for plotting with plottxt.m ifwr = 0 ! Set avmincrad for CMIT = 11.1 deg from fit of CP to crad_used ! Default avmincrad = 10 deg from min crit(1)=15 deg (crit(1)=crad+5deg) avmincrad = 11.1 ! mixdir06= "/hao/aim4/schmitt/data/LTR-2_1_4/Agu2006/CMITs/" ! mixdir06= "/hao/aim4/schmitt/data/LTR-2_1_5b/CMIT_Agu2006/CMITs/" mixdir06= "/hao/aim3/emery/cmit_dec06_hdf/" dmin = 5 ida0 = 13 imon = 12 nd0 = 347 nd1 = 347 nd2 = 350 nmph = 60/dmin ! Change values of nh0 and nm1 to get 1 time only for ifwr=1 starting at the begin time istep = 0 nd = 0; do ndnum=nd1,nd2 nd = nd + 1 ndmo = ida0 + nd1 - nd0 + nd - 1 ! Would get segmentation fault if nhr=5 and nm=3,12 in 0530, but find nm=3,5 or nm=6,12?? n23 = 23 nh0 = 0 do nhr=nh0,n23 ! Dec 2006 AGU storm is every 5 min do nm=nm1,nmph nmin = (nm-1)*dmin istep = istep + 1 write (filtim06,"('AGU_CMITs_mix_2006-12-',i2.2,'T',i2.2,'-',i2.2,'-00Z.hdf')") ndmo,nhr,nmin write (datehrmn,"('2006 ',i3.3,' 12 ',i2.2,1x,i2.2,1x,i2.2)") ndnum,ndmo,nhr,nmin ! write (filenam06,"(a47,a38)") mixdir06,filtim06 ! for 1.94.1 version ! write (filenam06,"(a53,a38)") mixdir06,filtim06 write (filenam06,"(a31,a38)") mixdir06,filtim06 ! write (*,*) "nd filenam06 =",nd,filenam06 call openMixFile(filenam06) call readVariable("Grid X", gridX) ! write(*,*) "Grid X dims", shape(gridX) call readVariable("Grid Y", gridY) ! write(*,*) "Grid Y dims", shape(gridY) call readVariable("Potential North [V]",potN) ! write(*,*) "Potential North [V] dims", shape(potN) call readVariable("Potential South [V]",potS) ! write(*,*) "Potential South [V] dims", shape(potS) ! call readVariable("FAC North [A/m^2]",facN) ! write(*,*) "FAC North [A/m^2] dims", shape(facN) ! call readVariable("FAC South [A/m^2]",facS) ! write(*,*) "FAC South [A/m^2] dims", shape(facS) ! call readVariable("Pedersen conductance North [S]",pedN) ! write(*,*) "Pedersen conductance North [S] dims", shape(pedN) ! call readVariable("Pedersen conductance South [S]",pedS) ! write(*,*) "Pedersen conductance South [S] dims", shape(pedS) ! call readVariable("Hall conductance North [S]",hallN) ! write(*,*) "Hall conductance North [S] dims", shape(hallN) ! call readVariable("Hall conductance South [S]",hallS) ! write(*,*) "Hall conductance South [S] dims", shape(hallS) call readVariable("Average energy North [keV]",engN) ! write(*,*) "Average energy North [keV] dims", shape(engN) call readVariable("Average energy South [keV]",engS) ! write(*,*) "Average energy South [keV] dims", shape(engS) call readVariable("Number flux North [1/cm^2 s]",fluxN) ! write(*,*) "Number flux North [1/cm^2 s] dims", shape(fluxN) call readVariable("Number flux South [1/cm^2 s]",fluxS) ! write(*,*) "Number flux South [1/cm^2 s] dims", shape(fluxS) ! call readVariable("Neutral wind speed North [m/s]",nwN) ! write(*,*) "Neutral wind speed North [m/s] dims", shape(nwN) ! call readVariable("Neutral wind speed South [m/s]",nwS) ! write(*,*) "Neutral wind speed South [m/s] dims", shape(nwS) call closeMixFile() !!! Do something with the data here!!! !!! Do something with the data here!!! !!! Do something with the data here!!! nji = shape(gridX) njlat = nji(1) nimlt = nji(2) ! write(*,*) "njlat nimlt =",njlat,nimlt ! Calculate rho and theta polar coordinates ! Calculate radius(27,181) 181 columns all the same, rows 1-27 from 0-.6947(0 to 39.8deg co-lat) do j=1,njlat do i=1,nimlt radius(j,i) = sqrt(gridX(j,i)**2 + gridY(j,i)**2) enddo enddo ! write(*,*) "radius1,27 =",radius(1,1),radius(njlat,1) do j=1,njlat colatrad(j) = radius(j,1) mlatdeg(j) = 90. - colatrad(j)/dtr ! colatdeg only used for printout where colat changed from 0-40 (radius) to 40-0deg ! mlat2pol is used with efluxS(N)H and epot colatdeg(j) = radius(28-j,1)/dtr mlat2pol(j) = 90. - radius(28-j,1)/dtr enddo if (nhr==0 .and. nm==1) write(*,*) "mlat2pol =",mlat2pol ! theta is 0 to +2pi, where 0=noon counterclockwise to later times for NH, ! but is clockwise to earlier times for SH ! Avoid atan2 for pole, since is a singularity - should be same for all j but j=1 j=njlat do i=1,nimlt theta(i) = atan2(gridY(j,i),gridX(j,i)) if (theta(i) < 0.) theta(i) = theta(i) + 2.*pii enddo ! write(*,*) theta(1),theta(nimlt) do i=1,nimlt mltn(i) = 12.+theta(i)/htr if (mltn(i)>=24.) mltn(i) = mltn(i)-24. ! 09/12 corrected mlts time going backwards from noon mlts(i) = 12.-theta(i)/htr; if (mlts(i)<0.) mlts(i) = mlts(i)+24. ! Want mlt=0-24 (theta always >=0 so goes 0-24) mlt(i) = theta(i)/htr enddo dmlt = (nimlt-1)/24. deltamlt = 1./dmlt ! if (nhr==0 .and. nm==1) write(*,*) "nimlt dmlt deltamlt MLT =", nimlt,dmlt,deltamlt,mlt ! Change epotn,s from 0-40 colat to 40-0 colat; add 180 deg to epotn so starts 0 MLT; ! change epotn so goes from 0-24 not noon-0-noon backwards in time ! Switch SH to later times is i(noon_later)=182-i(noon_earlier) ! Add 12 hours is +90 or -90 ! Also calculate efluxSH,NH from engS,N and fluxS,N using fac_p2e=1.602e-9 do i=1,181 do j=1,27 if (i .lt. 91) then epotSH(i,j) = potS(28-j,92-i) efluxSH(i,j) = engS(28-j,92-i)*fac_p2e*fluxS(28-j,92-i) epotNH(i,j) = potN(28-j,i+90) efluxNH(i,j) = engN(28-j,i+90)*fac_p2e*fluxN(28-j,i+90) else epotSH(i,j) = potS(28-j,272-i) efluxSH(i,j) = engS(28-j,272-i)*fac_p2e*fluxS(28-j,272-i) epotNH(i,j) = potN(28-j,i-90) efluxNH(i,j) = engN(28-j,i-90)*fac_p2e*fluxN(28-j,i-90) endif enddo enddo ! for i,j ! Calculate HP from efluxS,NH ! Re = radius of the Earth in m at the level of the ionosphere (use 6370e+3 in MIX) Re = 6371.e+3 areatot = 0. ! The MIX grid is 2 deg in maglon (irreg in maglat) dmlon = 360./180. do j=1,26 da(j) = Re*cos(mlat2pol(j)*dtr)*dmlon*dtr*(Re*dtr*(mlat2pol(j+1)-mlat2pol(j))) areatot = areatot + da(j)*180. enddo hpi_sh_mix = 0. hpi_nh_mix = 0. do i=1,180 hpi_sh_max = 0. hpi_nh_max = 0. do j=1,26 hpi_sh_mix = hpi_sh_mix + da(j)*efluxSH(i,j) hpi_nh_mix = hpi_nh_mix + da(j)*efluxNH(i,j) if (efluxSH(i,j) > hpi_sh_max) then hpi_sh_max = efluxSH(i,j) mlat_sh_max = mlat2pol(j) endif if (efluxNH(i,j) > hpi_nh_max) then hpi_nh_max = efluxNH(i,j) mlat_nh_max = mlat2pol(j) endif enddo ! do j=1,26 ! TEMP ! write (6,"(1x,'i mlt S,N_maxmlat S,N_eflux=',i3,3f7.2,2e12.4)") i,mlt(i),mlat_sh_max,mlat_nh_max,hpi_sh_max,hpi_nh_max enddo ! do i=1,180 ! Convert from mW/m2 * m2(da) * 1.e-12GW/mW to GW hpi_sh_mix = hpi_sh_mix*1.e-12 hpi_nh_mix = hpi_nh_mix*1.e-12 write (*,*) "HP SH,NH (GW) =",hpi_sh_mix,hpi_nh_mix ! TEMP write (61,"(i5,2f7.2,' MIX istep HPS HPN')") istep,hpi_sh_mix,hpi_nh_mix i02mlt = (nimlt-1)/2 + 2*dmlt i08mlt = (nimlt-1)/2 + 8*dmlt i14mlt = (14-12)*dmlt i20mlt = (21-12)*dmlt k02mlt = 2*dmlt + 1 k08mlt = 8*dmlt + 1 k14mlt = 14*dmlt + 1 k20mlt = 21*dmlt + 1 ! Find max/min Epot ! Find max Epot and position of max for 2-8 MLT ! Find min Epot and position of min for 14-20 MLT potnx = 0. potnn = 0. potnamx = 0. potnpmn = 0. potsx = 0. potsn = 0. potsamx = 0. potspmn = 0. do j=1,njlat do i=1,nimlt if (potN(j,i)>potnx*1.e+3) then potnx = potN(j,i)*1.e-3 latnx = mlatdeg(j) mltnx = mltn(i) if (abs(latnx)>89.99) mltnx = 6. inx = i jnx = j endif if (potN(j,i)89.99) mltnn = 18. inn = i jnn = j endif if (potS(j,i)>potsx*1.e+3) then potsx = potS(j,i)*1.e-3 latsx = mlatdeg(j) mltsx = mlts(i) if (abs(latsx)>89.99) mltsx = 6. isx = i jsx = j endif if (potS(j,i)89.99) mltsn = 18. isn = i jsn = j endif if (potN(j,i)>potnamx*1.e+3 .and. i>i02mlt .and. ii14mlt .and. ipotsamx*1.e+3 .and. i>k02mlt .and. ik14mlt .and. i