diff -crN gcpm_v22/Makefile gcpm_v22_r1/Makefile *** gcpm_v22/Makefile Thu Jan 1 09:00:00 1970 --- gcpm_v22_r1/Makefile Sat Feb 13 14:13:01 2010 *************** *** 0 **** --- 1,39 ---- + CC = gcc + F77 = gfortran + + CFLAGS = -O + FFLAGS = -O -fd-lines-as-comments -ffixed-line-length-none + + LIBRARY = libgcpm.a + + OBJS = bulge.o \ + ne_inner_ps_trough.o \ + switchon.o \ + gcpm_v22.o \ + ne_iri_cap.o \ + iri_ps_bridge.o \ + ne_iri_ps_trough.o \ + iri_ps_eq_bridge.o \ + ne_iri_ps_trough_eq.o \ + iri_sm.o \ + pp_profile.o \ + tec.o + + .SUFFIXES: .c .for + + .c.o: + $(CC) $(CFLAGS) -c $< + .for.o: + $(F77) $(FFLAGS) -c $< + + all: $(LIBRARY) + + $(LIBRARY): $(OBJS) + ar cru $(LIBRARY) $(OBJS) + ranlib $(LIBRARY) + + clean:; rm -f $(OBJS) *~ + + veryclean:; rm -f $(OBJS) $(LIBRARY) *~ + + library: $(LIBRARY) diff -crN gcpm_v22/gcpm_v22.for gcpm_v22_r1/gcpm_v22.for *** gcpm_v22/gcpm_v22.for Sat Feb 13 13:19:37 2010 --- gcpm_v22_r1/gcpm_v22.for Sat Feb 13 14:29:33 2010 *************** *** 165,171 **** c data oldmlt/-1.0/,oldkp/-1.0/ ! common /irioutput/ rz12,f107 c c altrans = the half width in L-shell of over which the transition c takes place between the trough and polar cap models. --- 165,171 ---- c data oldmlt/-1.0/,oldkp/-1.0/ ! common /irioutput/ rz12,f107,ne,nh,nhe,no c c altrans = the half width in L-shell of over which the transition c takes place between the trough and polar cap models. *************** *** 179,185 **** c invarient latitude of the polarward edge of the auroral zone. c This location is determined as a function of MLT and Kp from c the array PN. ! if(oldmlt.ne.amlt .or. oldkp.ne.akp) then oldmlt=amlt oldkp=akp bmlt=amlt*3.0+1.0 --- 179,185 ---- c invarient latitude of the polarward edge of the auroral zone. c This location is determined as a function of MLT and Kp from c the array PN. ! CCC if(oldmlt.ne.amlt .or. oldkp.ne.akp) then oldmlt=amlt oldkp=akp bmlt=amlt*3.0+1.0 *************** *** 210,216 **** tranlow=alcrit-altrans tranhigh=alcrit+altrans d type *,'auroral zone:',alcrit,altrans,tranlow,tranhigh ! endif c We need to obtain the L-shell of the location given, while limiting c the maximum L-shell that will be used. Higher latitudes and L-shells --- 210,216 ---- tranlow=alcrit-altrans tranhigh=alcrit+altrans d type *,'auroral zone:',alcrit,altrans,tranlow,tranhigh ! CCC endif c We need to obtain the L-shell of the location given, while limiting c the maximum L-shell that will be used. Higher latitudes and L-shells diff -crN gcpm_v22/iri_ps_bridge.for gcpm_v22_r1/iri_ps_bridge.for *** gcpm_v22/iri_ps_bridge.for Sat Feb 13 13:19:37 2010 --- gcpm_v22_r1/iri_ps_bridge.for Sat Feb 13 13:38:11 2010 *************** *** 38,52 **** c alatr=acos(cosrl) !Trevor Garner found error assuming north only, now pass lat along=amod((amlt+12.0),24.0)*amltrad cosrl=amin1(sqrt(r/al),1.0) ! alatrl=cos(cosrl)*ahemisphere call iri_sm(alatrl,along,r,itime,outf,oarr) r2=oarr(2)/re+1.0 cosrl=amin1(sqrt(r2/al),1.0) ! alatrl=cos(cosrl)*ahemisphere call iri_sm(alatrl,along,r2,itime,outf,oarr) r2=oarr(2)/re+1.0 cosrl=amin1(sqrt(r2/al),1.0) ! alatrl=cos(cosrl)*ahemisphere call iri_sm(alatrl,along,r2,itime,outf,oarr) c approximate the F2 peak along the L-shell=al --- 38,52 ---- c alatr=acos(cosrl) !Trevor Garner found error assuming north only, now pass lat along=amod((amlt+12.0),24.0)*amltrad cosrl=amin1(sqrt(r/al),1.0) ! alatrl=acos(cosrl)*ahemisphere call iri_sm(alatrl,along,r,itime,outf,oarr) r2=oarr(2)/re+1.0 cosrl=amin1(sqrt(r2/al),1.0) ! alatrl=acos(cosrl)*ahemisphere call iri_sm(alatrl,along,r2,itime,outf,oarr) r2=oarr(2)/re+1.0 cosrl=amin1(sqrt(r2/al),1.0) ! alatrl=acos(cosrl)*ahemisphere call iri_sm(alatrl,along,r2,itime,outf,oarr) c approximate the F2 peak along the L-shell=al *************** *** 85,91 **** r2=ah2/re+1.0 c get the density at the maximum slope height cosrl=amin1(sqrt(ro/al),1.0) ! alatrl=cos(cosrl)*ahemisphere call iri_sm(alatrl,along,ro,itime,outf,oarr) antransh=outf(1,1) --- 85,91 ---- r2=ah2/re+1.0 c get the density at the maximum slope height cosrl=amin1(sqrt(ro/al),1.0) ! alatrl=acos(cosrl)*ahemisphere call iri_sm(alatrl,along,ro,itime,outf,oarr) antransh=outf(1,1) *************** *** 95,105 **** c determine the ionospheric densities above and below to support initial c calculation of the power law function. cosrl=amin1(sqrt(r1/al),1.0) ! alatrl=cos(cosrl)*ahemisphere call iri_sm(alatrl,along,r1,itime,outf,oarr) an1=outf(1,1) cosrl=amin1(sqrt(r2/al),1.0) ! alatrl=cos(cosrl)*ahemisphere call iri_sm(alatrl,along,r2,itime,outf,oarr) an2=outf(1,1) --- 95,105 ---- c determine the ionospheric densities above and below to support initial c calculation of the power law function. cosrl=amin1(sqrt(r1/al),1.0) ! alatrl=acos(cosrl)*ahemisphere call iri_sm(alatrl,along,r1,itime,outf,oarr) an1=outf(1,1) cosrl=amin1(sqrt(r2/al),1.0) ! alatrl=acos(cosrl)*ahemisphere call iri_sm(alatrl,along,r2,itime,outf,oarr) an2=outf(1,1) diff -crN gcpm_v22/iri_sm.for gcpm_v22_r1/iri_sm.for *** gcpm_v22/iri_sm.for Sat Feb 13 13:19:37 2010 --- gcpm_v22_r1/iri_sm.for Sat Feb 13 14:30:22 2010 *************** *** 11,16 **** --- 11,17 ---- real pos_sm(3),pos_geo(3),blatr,blatd,blong,rtemp integer*4 ddd,jmag,itime(2),yyyy,i logical jf(30) + real buf(7) common /irioutput/ rz12,f107,ne,nh,nhe,no *************** *** 44,51 **** jmag=0 d type *,'iri called with:' d type *,blatd,blongd,yyyy,-ddd,dhour,aheight ! call iri_sub(jf,jmag,blatd,blongd,yyyy,-ddd,dhour, ! & aheight,aheight,delh,outf,oarr) d type *,'density returned=',outf(1,1),oarr(2) outf(1,1)=amax1(0.0,outf(1,1)) --- 45,68 ---- jmag=0 d type *,'iri called with:' d type *,blatd,blongd,yyyy,-ddd,dhour,aheight ! ! call check_iridb(blatd,blongd,yyyy,-ddd,dhour,aheight, ! & buf,istatus) ! if (istatus.eq.1) then ! outf(1,1) = buf(1) ! oarr(33) = buf(2) ! oarr(41) = buf(3) ! outf(6,1) = buf(4) ! outf(7,1) = buf(5) ! outf(5,1) = buf(6) ! oarr(2) = buf(7) ! else ! call iri_sub(jf,jmag,blatd,blongd,yyyy,-ddd,dhour, ! & aheight,aheight,delh,outf,oarr) ! call add_iridb(blatd,blongd,yyyy,-ddd,dhour,aheight, ! & outf,oarr) ! endif ! d type *,'density returned=',outf(1,1),oarr(2) outf(1,1)=amax1(0.0,outf(1,1)) *************** *** 62,64 **** --- 79,148 ---- return end + subroutine check_iridb(blatd,blongd,yyyy,ddd,dhour,aheight, + & buf,istatus) + implicit none + + real buf(7),blatd,blongd,dhour,aheight + integer*4 yyyy,ddd,istatus,i,j,ip + real db(2000,11) + integer idb(2000,2) + data ip/0/ + common /iridb/ db,idb,ip + + istatus = 0 + do i=1,2000 + if (blatd.eq.db(i,1).and. + & blongd.eq.db(i,2).and. + & dhour.eq.db(i,3).and. + & aheight.eq.db(i,4).and. + & yyyy.eq.idb(i,1).and. + & ddd.eq.idb(i,2)) then + buf(1) = db(i,5) + buf(2) = db(i,6) + buf(3) = db(i,7) + buf(4) = db(i,8) + buf(5) = db(i,9) + buf(6) = db(i,10) + buf(7) = db(i,11) + istatus = 1 + endif + enddo + + return + end + + subroutine add_iridb(blatd,blongd,yyyy,ddd,dhour,aheight, + & outf,oarr) + implicit none + + real outf(20,100),oarr(50) + real blatd,blongd,dhour,aheight + integer*4 yyyy,ddd,ip + real db(2000,11) + integer idb(2000,2) + common /iridb/ db,idb,ip + + ip = ip+1 + if (ip.gt.2000) then + ip = 1 + endif + + db(ip,1) = blatd + db(ip,2) = blongd + db(ip,3) = dhour + db(ip,4) = aheight + + idb(ip,1) = yyyy + idb(ip,2) = ddd + + db(ip,5) = outf(1,1) + db(ip,6) = oarr(33) + db(ip,7) = oarr(41) + db(ip,8) = outf(6,1) + db(ip,9) = outf(7,1) + db(ip,10)= outf(5,1) + db(ip,11)= oarr(2) + + return + end diff -crN gcpm_v22/ne_inner_ps_trough.for gcpm_v22_r1/ne_inner_ps_trough.for *** gcpm_v22/ne_inner_ps_trough.for Sat Feb 13 13:19:37 2010 --- gcpm_v22_r1/ne_inner_ps_trough.for Sat Feb 13 14:31:22 2010 *************** *** 27,32 **** --- 27,33 ---- integer*4 icount data a6old/0.0/,a7old/0.0/,a8old/0.0/ + common /geo/ geosync_trough d type *,'ne_eq_trough caled with',al,amlt,akp c *************** *** 128,139 **** b1=a7 c return c This routine checks to make sure that a8 is not beyond the point where c the inner plasmasphere and trough density models would cross were there c no plasmapause. If a8 is beyond that point, then the crossing point is c substituted for a8. ! entry check_crossing(a8,am1,b1,amlt,akp) d type *,'initial crossing=',a8,am1,b1,amlt,akp c Determine where the inner plasmasphere plus plasmapause profile --- 129,142 ---- b1=a7 c return + end c This routine checks to make sure that a8 is not beyond the point where c the inner plasmasphere and trough density models would cross were there c no plasmapause. If a8 is beyond that point, then the crossing point is c substituted for a8. ! function check_crossing(a8,am1,b1,amlt,akp) ! common /geo/ geosync_trough d type *,'initial crossing=',a8,am1,b1,amlt,akp c Determine where the inner plasmasphere plus plasmapause profile *************** *** 157,165 **** d type *,'crossing:',zl,stepl,diff icount=icount+1 if (icount.gt.100) then ! type *,'check_crossing is loop-bound:',am1,b1,zl,amlt,akp,a8, & geosync_trough,stepl ! type *,'STOPPING***********' stop endif enddo --- 160,168 ---- d type *,'crossing:',zl,stepl,diff icount=icount+1 if (icount.gt.100) then ! write(*,*) 'check_crossing is loop-bound:',am1,b1,zl,amlt,akp,a8, & geosync_trough,stepl ! write(*,*) 'STOPPING***********' stop endif enddo diff -crN gcpm_v22/ne_iri_ps_trough.for gcpm_v22_r1/ne_iri_ps_trough.for *** gcpm_v22/ne_iri_ps_trough.for Sat Feb 13 13:19:37 2010 --- gcpm_v22_r1/ne_iri_ps_trough.for Sat Feb 13 13:55:37 2010 *************** *** 34,42 **** d type *,'back from ne_iri_ps_trough_eq=',eq_iri_ps_trough d type *,'bridge? ',amlt,amlt_o,akp,akp_o,al,al_o d type *,' ? ',itime(1),itime1_o,itime(2),itime2_o ! if(amlt.ne.amlt_o .or. akp.ne.akp_o .or. ! & itime(1).ne.itime1_o .or. itime(2).ne.itime2_o .or. ! & abs(al-al_o).gt.1.0e-5) then d type *,'Evaluate field aligned power law function' d type *,'showit:',(amlt-amlt_o),(akp-akp_o),(al-al_o) d type *,(itime(1)-itime1_o),(itime(2)-itime2_o) --- 34,42 ---- d type *,'back from ne_iri_ps_trough_eq=',eq_iri_ps_trough d type *,'bridge? ',amlt,amlt_o,akp,akp_o,al,al_o d type *,' ? ',itime(1),itime1_o,itime(2),itime2_o ! CCC if(amlt.ne.amlt_o .or. akp.ne.akp_o .or. ! CCC & itime(1).ne.itime1_o .or. itime(2).ne.itime2_o .or. ! CCC & abs(al-al_o).gt.1.0e-5) then d type *,'Evaluate field aligned power law function' d type *,'showit:',(amlt-amlt_o),(akp-akp_o),(al-al_o) d type *,(itime(1)-itime1_o),(itime(2)-itime2_o) *************** *** 60,66 **** itime2_o=itime(2) d type *,'ne_iri_ps_trough transh=',transh,alpha,dno,co,istat d type *,' >',r,al,alatr,amlt,itime,eq_iri_ps_trough ! end if c compute density as given by the bridge function eq_bridge=dno*aheight**(-alpha) + co --- 60,66 ---- itime2_o=itime(2) d type *,'ne_iri_ps_trough transh=',transh,alpha,dno,co,istat d type *,' >',r,al,alatr,amlt,itime,eq_iri_ps_trough ! CCC end if c compute density as given by the bridge function eq_bridge=dno*aheight**(-alpha) + co diff -crN gcpm_v22/ne_iri_ps_trough_eq.for gcpm_v22_r1/ne_iri_ps_trough_eq.for *** gcpm_v22/ne_iri_ps_trough_eq.for Sat Feb 13 13:19:37 2010 --- gcpm_v22_r1/ne_iri_ps_trough_eq.for Sat Feb 13 13:53:20 2010 *************** *** 38,45 **** ps_inner=ne_inner_ps(r,amlt,am1,b1)*1.0e6 d type *,'from ne_inner_ps:',ps_inner,r,amlt,am1,b1 ! if(amlt.ne.amlt_o .or. akp.ne.akp_o .or. ! & itime(1).ne.itime1_o .or. itime(2).ne.itime2_o) then c c determine the height power law fit parameters that connect the c topside ionosphere to the plasmasphere --- 38,45 ---- ps_inner=ne_inner_ps(r,amlt,am1,b1)*1.0e6 d type *,'from ne_inner_ps:',ps_inner,r,amlt,am1,b1 ! C if(amlt.ne.amlt_o .or. akp.ne.akp_o .or. ! C & itime(1).ne.itime1_o .or. itime(2).ne.itime2_o) then c c determine the height power law fit parameters that connect the c topside ionosphere to the plasmasphere *************** *** 52,58 **** akp_o=akp itime1_o=itime(1) itime2_o=itime(2) ! end if d type *,'values 1:',a8,am1,b1 ps_bridge=ano*aheight**(-alpha) --- 52,58 ---- akp_o=akp itime1_o=itime(1) itime2_o=itime(2) ! C end if d type *,'values 1:',a8,am1,b1 ps_bridge=ano*aheight**(-alpha) diff -crN gcpm_v22/pp_profile.for gcpm_v22_r1/pp_profile.for *** gcpm_v22/pp_profile.for Sat Feb 13 13:19:37 2010 --- gcpm_v22_r1/pp_profile.for Sat Feb 13 13:52:04 2010 *************** *** 24,31 **** d type *,'into pp_profile:',al,amlt,akp,a8 c c Allow for mlt rotation of the buldge with Kpmax ! if((akp.ne.akp_old) .or. (amlt.ne.amlt_old)) ! & call bulge(amlt,akp,a8,a9,centroid) akp_old=akp amlt_old=amlt c --- 24,32 ---- d type *,'into pp_profile:',al,amlt,akp,a8 c c Allow for mlt rotation of the buldge with Kpmax ! CCC if((akp.ne.akp_old) .or. (amlt.ne.amlt_old)) ! CCC & call bulge(amlt,akp,a8,a9,centroid) ! call bulge(amlt,akp,a8,a9,centroid) akp_old=akp amlt_old=amlt c diff -crN gcpm_v22/tec.c gcpm_v22_r1/tec.c *** gcpm_v22/tec.c Thu Jan 1 09:00:00 1970 --- gcpm_v22_r1/tec.c Sat Feb 13 13:05:24 2010 *************** *** 0 **** --- 1,58 ---- + #include + + #define RE (6370.0) + #define NREG (5) + + void atec(int itime[2], float amlt, float alatr, float akp, float *upto500, float *upto1000, float *upto20000){ + + float ec; + float rr,r; + float outn[4]; + float rstart[NREG]={ 0.0, 500.0, 1000.0, 2000.0, 10000.0 }; + float rend [NREG]={ 500.0, 1000.0, 2000.0, 10000.0, 20000.0 }; + float rstep [NREG]={ 10.0, 20.0, 50.0, 100.0, 1000.0 }; + int i,ist,ied; + + *upto500 = 0.0; + *upto1000 = 0.0; + *upto20000 = 0.0; + for(i=0; i