Event Study Workflow

From OpenGGCM

Contents

Getting and using solar wind driver data

ACE Merged IMF and Solar Wind 64-second Averages: http://www.srl.caltech.edu/ACE/ASC/level2/index.html

1) Click on "MAG/SWEPAM Data"
2) Select a year
3) Click on "MAG and SWEPAM 64-second Averaged IMF and Solar Wind Ion Parameters"
4) Enable the box "check here to select all fields"
5) Properly fill out the time range boxes:
   Starting "YR/DOY":      Ending "YR/DOY":
   If you need the Day Of Year (DOY) for a particular date, you can use the perl utility:
   > ~/GGCM/trunk/util/perl/day2doy.pl
   Month > 8
   Day > 17
   Year > 2003
   HH:MM:SS:MS > 00:00:00:00

   User Input: 8/17/2003 0:0:0:0
   Computed: FP Day Of Year Number = 229
   Computed: FP Year = 2003.62739726027

   Alternatively, the built in UNIX utility cal can be used with the '-j' option:
   > cal 08 2003 -j
           August 2003
   Sun Mon Tue Wed Thu Fri Sat
                       213 214
   215 216 217 218 219 220 221
   222 223 224 225 226 227 228
   229 230 231 232 233 234 235
   236 237 238 239 240 241 242
   243

   On the ACE web page, continue by selecting 'text file download'
   Choose Data Format:      X-Y Plot      Text file download

6) Retrieve Data and save to your local disk.

Make a directory for the data file, because the next process will split the
downloaded file into many little files.

> mkdir event.2003.08.17
Move your downloaded data into this directory. The ACE convention is that it is
named 'ACE_MAGSWE_Data.txt'

Now copy this master data file into a working name, many of the OpenGGCM scripts
use a simple prefix to denote a particular satllite, e.g. 'ac' for ACE, and 'wi'
for WIND, etc.

> cp ACE_MAGSWE_Data.txt ac.txt

> ~/GGCM/trunk/util/perl/satfilter.pl --style=ggcm --sat=ACE ac.txt

A lot of stuff should fly by on the screen. When your are done:
~/event.2003.08.17> ls
ac.ACEepoch     ac.Bmag              ac.fp_doy     ac.pos_gse_y  ac.sec      ac.V_gsm_z
ac.Alpha_ratio  ac.B_rtn_n           ac.fp_year    ac.pos_gse_z  ac.Tp       ac.Vp
ac.B_gse_x      ac.B_rtn_r           ac.hr         ac.pos_gsm_x  ac.txt      ac.V_rtn_n
ac.B_gse_y      ac.B_rtn_t           ac.Lambda     ac.pos_gsm_y  ac.V_gse_x  ac.V_rtn_r
ac.B_gse_z      ac.day               ac.MAG_pts    ac.pos_gsm_z  ac.V_gse_y  ac.V_rtn_t
ac.B_gsm_x      ac.dBrms             ac.min        ac.pos_hs_x   ac.V_gse_z  ac.year
ac.B_gsm_y      ac.Delta             ac.Np         ac.pos_hs_y   ac.V_gsm_x
ac.B_gsm_z      ACE_MAGSWE_Data.txt  ac.pos_gse_x  ac.pos_hs_z   ac.V_gsm_y

The next preprocessing script requires that your data files have lowercase extensions.
In addition, these are the required files:
ac.np
ac.rr
ac.pp
ac.temp
ac.bxgse
ac.bygse
ac.bzgse
ac.btot
ac.vxgse
ac.vygse
ac.vzgse
ac.vtot
ac.vth
ac.xgse
ac.ygse
ac.zgse

In this example we are using ACE, so 'ac' is the prefix. If you use a different satellite,
you will have a different prefix, e.g. GEOTAIL uses 'ge'.

So let us make and edit a script:
> /bin/ls -1 | awk '{print "mv "$1 OFS $1}' > tmp.rename
  Note: If you have an alias for 'ls', it is a good idea to use /bin/ls because your
  alias may enter unexpected characters into your tmp file.
  After editing:
> cat tmp.rename
mv ac.B_gse_x ac.bxgse
mv ac.B_gse_y ac.bygse
mv ac.B_gse_z ac.bzgse
mv ac.Bmag ac.btot
mv ac.Np ac.np
mv ac.pos_gse_x ac.xgse
mv ac.pos_gse_y ac.ygse
mv ac.pos_gse_z ac.zgse
mv ac.Tp ac.temp
mv ac.V_gse_x ac.vxgse
mv ac.V_gse_y ac.vygse
mv ac.V_gse_z ac.vzgse
mv ac.Vp ac.vtot

Execute this script:
> sh tmp.rename

There are probably more clever shell tricks than just editing a file, but it is quick.

The following is a simple fortran code used by Wenhui Li to make the solar wind data files. 
You must edit it for your own use. When download the ACE data, select year,day,hour,min,sec,
H+ density, H+ temp, H+ speed,Vx_gse,Vy_gse, Vz_gse,Bx_gse,By_gse,Bz_gse,gse_x,gse_y,gse_z.
--------------------------------------------------------------------------------------
c  makeMinvarFiles.f

      program makeMinvarFiles

      integer npoint, year(100000), mon(100000), day(100000)
      integer RDSTAT, hh(100000), mm(100000)
      integer mon1, day1, hhsimf
      real sec(100000), Np(100000), Tp(100000), Vp(100000)
      real Vx(100000), Vy(100000), Vz(100000), Bx(100000)
      real By(100000), Bz(100000), x(100000), y(100000), z(100000)
      character dummy*160, infile*80, outf(30)*80

      mon1 = 8
      day1 = 31
      infile = '/home/wenhuil/eventInput/2003Aug31/ACE_MAGSWE_Data.txt'
      OPEN(UNIT=10, FILE=infile, STATUS='OLD')
      RDSTAT = 0
      DO i=1, 50
         read(10, *) dummy
         write(0,*) dummy
      ENDDO

      i = 0
      DO WHILE (RDSTAT .EQ. 0)
         i = i + 1
         read(10, *, IOSTAT=RDSTAT) year(i),day(i),hh(i),mm(i),
     *        sec(i),Np(i),Tp(i),Vp(i),Vx(i),Vy(i),Vz(i),Bx(i),By(i),
     *        Bz(i),x(i),y(i),z(i)
c         write(0,*) year(i),day(i),hh(i),mm(i),sec(i),Np(i),Tp(i)
      ENDDO
      npoint = i-1
      CLOSE(UNIT=10)
                                                                
      outf(1) = '/home/wenhuil/eventInput/2003Aug31/ac.np'
      outf(2) = '/home/wenhuil/eventInput/2003Aug31/ac.temp'
      outf(3) = '/home/wenhuil/eventInput/2003Aug31/ac.vtot'
      outf(4) = '/home/wenhuil/eventInput/2003Aug31/ac.vxgse'
      outf(5) = '/home/wenhuil/eventInput/2003Aug31/ac.vygse'
      outf(6) = '/home/wenhuil/eventInput/2003Aug31/ac.vzgse'
      outf(7) = '/home/wenhuil/eventInput/2003Aug31/ac.bxgse'
      outf(8) = '/home/wenhuil/eventInput/2003Aug31/ac.bygse'
      outf(9) = '/home/wenhuil/eventInput/2003Aug31/ac.bzgse'
      outf(10) = '/home/wenhuil/eventInput/2003Aug31/ac.xgse'
      outf(11) = '/home/wenhuil/eventInput/2003Aug31/ac.ygse'
      outf(12) = '/home/wenhuil/eventInput/2003Aug31/ac.zgse'

      OPEN(UNIT=12, FILE=outf(1), STATUS='NEW')
      do i=1, npoint
         if (Np(i) .GE. 0.0) then
            write(12, *)
     *             year(i),mon1,day1,hh(i),mm(i),sec(i),Np(i)
         endif
      enddo
      CLOSE(UNIT=12)
      OPEN(UNIT=12, FILE=outf(2), STATUS='NEW')
      do i=1, npoint
         if (Tp(i) .GE. 0.0) then
             write(12, *)
     *             year(i),mon1,day1,hh(i),mm(i),sec(i),Tp(i)
         endif
      enddo
      CLOSE(UNIT=12)

      OPEN(UNIT=12, FILE=outf(3), STATUS='NEW')
      do i=1, npoint
         if (Vp(i) .GE. 0.0) then
             write(12, *)
     *             year(i),mon1,day1,hh(i),mm(i),sec(i),Vp(i)
         endif
      enddo
      CLOSE(UNIT=12)

      OPEN(UNIT=12, FILE=outf(4), STATUS='NEW')
      do i=1, npoint
         if (Vx(i) .GE. -2000.0) then
             write(12, *)
     *             year(i),mon1,day1,hh(i),mm(i),sec(i),Vx(i)
         endif
      enddo
      CLOSE(UNIT=12)

      OPEN(UNIT=12, FILE=outf(5), STATUS='NEW')
      do i=1, npoint
         if (Vy(i) .GE. -1000.0) then
             write(12, *)
     *             year(i),mon1,day1,hh(i),mm(i),sec(i),Vy(i)
         endif
      enddo
      CLOSE(UNIT=12)

      OPEN(UNIT=12, FILE=outf(6), STATUS='NEW')
      do i=1, npoint
         if (Vz(i) .GE. -1000.0) then
             write(12, *)
     *             year(i),mon1,day1,hh(i),mm(i),sec(i),Vz(i)
         endif
      enddo
      CLOSE(UNIT=12)

      OPEN(UNIT=12, FILE=outf(7), STATUS='NEW')
      do i=1, npoint
         if (Bx(i) .GE. -200.0) then
             write(12, *)
     *             year(i),mon1,day1,hh(i),mm(i),sec(i),Bx(i)
         endif
      enddo
      CLOSE(UNIT=12)

      OPEN(UNIT=12, FILE=outf(8), STATUS='NEW')
      do i=1, npoint
         if (By(i) .GE. -200.0) then
             write(12, *)
     *             year(i),mon1,day1,hh(i),mm(i),sec(i),By(i)
         endif
      enddo
      CLOSE(UNIT=12)

      OPEN(UNIT=12, FILE=outf(9), STATUS='NEW')
      do i=1, npoint
         if (Bz(i) .GE. -200.0) then
             write(12, *)
     *             year(i),mon1,day1,hh(i),mm(i),sec(i),Bz(i)
         endif
      enddo
      CLOSE(UNIT=12)

      OPEN(UNIT=12, FILE=outf(10), STATUS='NEW')
      do i=1, npoint
         if (x(i) .GE. -1000.0) then
             write(12, *)
     *             year(i),mon1,day1,hh(i),mm(i),sec(i),x(i)/6378.0
         endif
      enddo
      CLOSE(UNIT=12)

      OPEN(UNIT=12, FILE=outf(11), STATUS='NEW')
      do i=1, npoint
         write(12, *)
     *        year(i),mon1,day1,hh(i),mm(i),sec(i),y(i)/6378.0
      enddo
      CLOSE(UNIT=12)

      OPEN(UNIT=12, FILE=outf(12), STATUS='NEW')
      do i=1, npoint
         write(12, *)
     *        year(i),mon1,day1,hh(i),mm(i),sec(i),z(i)/6378.0
      enddo
      CLOSE(UNIT=12)


      end

--------------------------------------------------------------------





You will now use a csh script. Here is the call chain of mak.prep csh script:
/usr/local/public/bin/mak.prep
   --> /usr/local/public/bin/amak.derived
      --> tser_arith, a binary file generated by
          /usr/local/public/src/tser_arith.com
          which depends on compile, and fppn
All of the call chain must exist in order for this to run correctly.

Now run mak.prep in our example directory:

~/event.2003.08.17> /usr/local/public/bin/mak.prep
Tool path is /usr/local/public/bin
Spacecraft prefix (ac, cl, g0, g8, ge, l0, l4, l7, l9, ma, po, sa, wi) ?
ac
WARNING: There will be no checks on your Timecode entries.
T1 (YEAR:MO:DD:HR:MM:0.0) = ?
2003:08:17:00:00:0.0
T2 (YEAR:MO:DD:HR:MM:0.0) = ?
2003:08:19:23:59:0.0
Fillgap, in seconds ?
30
Spacecraft= ac T1=2003:08:17:00:00:0.0 T2=2003:08:19:23:59:0.0 Fillgap=30

You should now have a bunch of files all with the same 30 second cadence.

Leave this folder, but keep track of it's exact path. This path and the folder
name will be assigned to the IDATA variable in your 'runme' file. When the
runme is executed, it will look in this folder to construct the solar wind
input for the simulation.

Getting and using satellite orbital data

Navigate your web browser to: http://sscweb.gsfc.nasa.gov/cgi-bin/sscweb/Locator.cgi

   Select Satellites and the start-stop range
   Click the "Ouput Settings" and select GSE coordinates in the XYZ  column.
   Click "Output Units/Filtering"
   Select yy/mm/dd and hh:mm:ss
   Enter "8" for the number of decimal places.
   Click "Submit query and wait for output"

Copy the ASCII output from your screen to a file, like  sscweb_yymmdd.txt

Edit this file by deleting from the top until "GROUP I," about 30  lines.
Then delete all the junk text at the bottom of the file.

Copy the sscweb_yymmdd.txt file into the run folder with the name  sscweb.txt

Make sure you have a valid runme in the same folder as the sscweb.txt
data file and execute the following command:
~/GGCM/branches/OpenGGCM_v3.1/run-template/mak.orbits

This perl script will generate a file called orbits.txt that will be
read by OpenGGCM during the run. Check to make sure it is not empty  and
that it looks like valid satllite data. If the file is empty, your  runme
may have an incorrect time range, or you downloaded the wrong range  from
the SSC website.

After the event run has completed, you must postprocess the  $RUN.timehistNNNNN
files that will be located in the data directory. This should be  done on a
storage node since the sorting process is memory intensive.
~/GGCM/branches/OpenGGCM_v3.1/run-template/mak.timehistplot --run  $RUN --sort

Edit your 'runme'

## ..... automatic SW file generation (comment section out if not used)
SWMON      ac           # solar wind monitor sc (wi/ac/i8/ge/fi), use 'fi' if none
SWFILE     minvar       # file with SW/IMF, can be auto/minvar
SETIMFBX   none         #  set Bx to fixed value, or "none"
INPAVE     30           #  averaging interval for input files, seconds
IDATA      /The/Full/Pathname/event.2003.08.17

Postprocessing

Read more about this software here and about how to get Jimmy's Software

Plotting the standard plots.

  * Batch driver for amak.ts_fix01, amak.ioseries2, amak.pl_ioseries2, amak.ts_fix01

> mak.std_*

Do the Eparallel etc stuff (Optional)

  * Batch driver for ftr4c & iocl4
  Does the field line tracing in 3df files to get various quantities
  This is *really* time consuming since it queues a job for each timestep.

> mak.3d-pp-batch01

Compute the Open-Closed boundaries from iof files

   produces io_opcl variable ionosphere field
   io_opcl=-1 for closed field, io_opcl=1 for
  * Batch driver for iocl3

> ~/bin/mak.ioopcl

Extract Open-Closed boundaries as function of MLT

  * Batch driver for iocl3-mlt

> ~/bin/mak.OpenClosed

process_iof.sh