(nonexistent sample file)

Initial Positions File

The initial positions file contains a separate line for each desired trajectory. Each line contains three spatial indices (U, V, W), a time index, and a fifth ignorable yet required paramter.

A sample initial_positions.txt file is:

      29.500  13.000   5.500  1.000    1.0
      29.500  13.000   6.500  1.000    1.0
      29.500  13.000   7.500  1.000    1.0

Tracmass

There is a PDF manual at the site, i.e.

http://doos.misu.su.se/tracmass/tracmass_documentation.pdf

which isn't terribly helpful in figuring out the variables in the two required input files. There is a PDF presentation included in the source code distribution which also isn't terribly helpful.

There is a much more helpful documentation source within the source code itself, i.e.

/raid2/GOMEX/lagrange/tracmass/tracmass/tracsvn/src/comment.f95

which contains much more practical information than either of the PDF documents.

The Quick and Dirty Method for Compilation

Note: These instructions are for a 64-bit Linux platform using the gfortran compiler. This will be different for other combinations.

Obtain the Tracmass source code package at:

http://doos.misu.su.se/tracmass/tracmass.zip

Unzip the bundle, e.g.

unzip tracmass.zip

to create a subdirectory named tracmass. Go to the tracsvn subdirectory, e.g.

cd tracsvn

Make the following changes - which may or may not be needed as explained below - in the Makefile and various Fortran source code files.

• Find an appropriate subdirectory under transcvn/projects, or create your own and figure out how to modify everything correctly for your case using the essentially undocumented package. I'll attempt to explain this below this list.

• Change the PROJECT variable to the name of the transcvn/projects subdirectory you have selected.

• Go to the tracsvn/src directory and edit the files getfile.f95, sw_stat.f95, sw_seck.f95, sw_dens0.f95 and sw_pres.f95. Replace all occurrences of ALLOCATABLE with POINTER. I had to do it using gfortran/gcc version 4.1.2 and the reason why is explained at:

  http://www.tek-tips.com/viewthread.cfm?qid=1618900

  Your mileage may vary. If your Fortran compiler doesn't complain about the ALLOCATABLE statements, then ignore this. A brief summary of why this may be necessary is:

  In fact, it really depends on your compiler. Using the ALLOCATABLE attribute for a dummy argument is a FORTRAN-2003 feature. That feature is a common extension of FORTRAN-95 compilers because it was the object of a Technical Recommendation (TR15581 in 1998). So most of recent FORTRAN-95 compilers implement it. But if you have a strict FORTRAN-95 compiler, then you are right : only the POINTER attribute is allowed for a dummy argument.

• If you're using a 64-bit rather than a 32-bit machine, you need to change the Makefile line:

  FF_FLAGS = -c -x f95-cpp-input -fconvert=big-endian -gdwarf-2 -m32

  to:

  FF_FLAGS = -c -x f95-cpp-input -fconvert=big-endian -gdwarf-2 -m64

• Change the LIB_DIR and INC_DIR lines to reflect the location of your NetCDF library, e.g. if your NetCDF library is installed at /opt/netcdf4 then change:

  LIB_DIR = -L/sw/lib -L/sw/lib/netcdf-gfortran/lib
  INC_DIR = -I/sw/include -I/sw/lib/netcdf-gfortran/include
  

  to

  LIB_DIR = -L/opt/netcdf4/lib
  INC_DIR = -I/opt/netcdf4/include

• Create a projects/gomex subdirectory and copy all of the tes subdirectory files into that. Change the Makefile pointer to gomex, and compile.

• Start modifying the gomex directory Makefile. First, comment out fl16 = -Dvarbottombox and try to compile. You get an error about dztb that can be corrected by commenting out the occurrences of dztb in the readfield.f95 file. Enter make clean and then make to successfully compile.

• Next comment out fl17 = -Dfreesurface in the gomex directory Makefile. Enter make clean and then make to successfully compile.

• Copy tes_grid.in to gomex_grid.in and make the following changes for the GOMEX project:

    GCMname     = 'GOMEX Project GCMs',
    ! GCMsource   = 'http://www.ccsm.ucar.edu/models/ccsm3.0/pop/',
    gridName    = 'GOMEX Standard Grid',
    ! gridSource  = 'http://www.whoi.edu/sbl/liteSite.do?litesiteid=23412&articleId=35629',
    gridDesc    = '0.05 deg. grid in GOM from -98 to -79 W and 18 to 31 N'/
  ...
    IMT      = 401,
    JMT      = 261,
    KM       = 1,
  ...
    NTRACMAX=1000/
  ...
    ngcm          = 24,             ! (5*24) hours between GCM datasets
    iter          = 5,              ! iteration between two gcm data sets
    intmax        = 365,            ! maximum length of RCO fields
    fieldsPerFile = 365/
  ...
    yearmin        = 0,
    yearmax        = 0,
  ...
    baseYear       = 2010/
  

  Enter make clean and then make to successfully compile.

• Copy tes_run.in to gomex_run.in and make the following changes in the INITRUNSEED section.

    nff      = 3,
    isec     = 3,
    idir     = 0,
    nqua     = 1, 
    partQuant= 1,    !(particles/gridcell or m3s-1/particle or m3/particle)
  

  Change:

    SeedType = 1,
  

  to

    SeedType = 2,
    varSeedFile = 2,
    seedFile = "/raid2/GOMEX/lagrange/drifters/drifters_seed.txt"/
  

  and then comment out the remainder of the INITRUNSEED section, i.e.

    !ist1     =   1,
    !ist2     =  -1,
    !jst1     =   1,
    !jst2     =  -1,
    !kst1     =   1,
    !kst2     =  -1/
  

  Enter make clean and then make to successfully compile.

• Delete the entire INITRUNTEMPSALT, INITRUNDIFFUSION and INITRUNSEDIMENT sections. Change the following line:

  timax =365./ ! maximum time length of a trajectory in days

  Enter make clean and then make to successfully compile.

• Remove all entries in the INITRUNGRID section and replace with:

    subGrid      = 0,
    subGridID    = 1/
  

  Enter make clean and then make to successfully compile.

• Change the INITRUNTIME, INITRUNDATE and INITRUNWRITE sections to:

  &INITRUNTIME
    intmin       = 1,
    intspin      = 1,
    intrun       = 365,
    intstep      = 1/
  &INITRUNDATE
    startHour         = 0,
    startDay          = 1,
    startMon          = 1,
    startYear         = 2010/
  &INITRUNWRITE
    ncoor   = 1,      
    kriva    = 1,
    inDataDir    = '/raid2/GOMEX/STEP1/1_IASNFS_Ko/tracmass/in',
    outDataDir   = '/raid2/GOMEX/STEP1/1_IASNFS_Ko/tracmass/out',
    outDataFile     = 'gomex_1_IASNFS'/
  

  Enter make clean and then make to successfully compile.

• Copy the projects/ncarCSM/readfield.f95 file to the gomex subdirectory. Edit it in the following ways.
  • Remove all the #ifdef tempsalt sections.
  • Set the value of the variable dstamp to H_iasroms_201001011200-surf.nc.
  • Modify the following statements as follows:

      uvel = get3DfieldNC(trim(fieldFile), 'eastward_sea_water_velocity') / 100
      vvel = get3DfieldNC(trim(fieldFile), 'northward_sea_water_velocity') / 100
    ...
      do k=1,km
         kk=km+1-k
         uflux(:,:,k,2) = uvel(:,:,kk) * dy * dz(k)
         vflux(:,:,k,2) = vvel(:,:,kk) * dx * dz(k)
      end do
    

• Edit the src/coord.f95 file and add:

  #ifdef gomex
  dx = 0.05d0
  dy = 0.05d0
  stlon1=0.
  stlat1=0.
  rmin=20.d0
  rmax=30.d0
  smin= 0.d0
  smax=40.d0
  tmin= 0.d0
  tmax=30.d0
  #endif
  
  #ifdef gomex
  dz = 100.0d0
  do i=1,imt
   do j=1,jmt
    dxdy(i,j)=dx*deg*dy*deg
   enddo
  enddo
  #endif
  

• The 2-D array kmt - in addition to specifying the number of vertical levels - also serves as a land/ocean point discriminator. This construct in loop.f95:

         case ('landError')
            if(kmt(ib,jb).eq.0) then
               if (verboseMess == 1) then
                  print *,'====================================='
                  print *,'Warning: Trajectory on land'
  

  shows us this. Thus we must define kmt in coord.f95 thusly:

  kmt=KM ! all grid cells are active in an AGCM

NCO Tricks Needed

To extract the first time step from H_iasroms_201001011200-surf.nc:

ncks -H -d time,0 H_iasroms_201001011200-surf.nc fred.nc

To remove the now redundant time dimension from fred.nc:

ncwa -a time fred.nc bob.nc

To convert the negative depths in model 5 into positive depths:

ncap2 -s "zlev=zlev*-1" MITGOM_HCAST_2010_v2.2-mod.nc MITGOM_HCAST_2010_v2.2-mod.nc

Note: The program convert.py was first written thinking that there was a binary problem and that reading and then writing it using Numpy would solve that. During that process, I noticed that the depths were defined as negative.

Pre-Processing Options

Many of the features of Tracmass are established via pre-processing options. The option in boldface is the one we will use in each category.

Project Options

These are options that select a specific project with a specific type of grid.

 -Dgomex       GOMEX project
 -Dorca        NEMO with one of the following ORCA grids
 -Dorca2       NEMO with 2    degree ORCA2 grid 
 -Dorca1       NEMO with 1    degree ORCA1 grid
 -Dorca05      NEMO with 0.5  degree ORCA05 grid
 -Dorca025     NEMO with 0.25 degree ORCA025 grid
 -Dorca12      NEMO with 1/12 degree ORCA12 grid
 -Docc         OCCAM 66 levles 1/4 deree resolution
 -Drco         RCO model 2 nm
 -Dfors        Forsmark model
 -Dsimp        Simpevarp model
 -Dtes         Test basin with analytical time depenent velocity fields
 -Dtun         ROMS GCM for the Med 
 -Datm         The AGCM (IFS) from ECMWF, which is part of EC-Earth

Time Options

Options for the temporal nature of the velocity fields to be used.

 -Dstat         Stationary velocity fields (does not work at the moment)
 -Dtimestep     Time changing velocity fields but using stationary solution to the differential
                equations with time step changing velocities
 -Dtimeanalyt   Time changing velocity fields with time analytical solutions to the differential
                equations based on Vries and Döös (2001). This scheme is not as robust as
                -Dtimestep.

Output Options

Options to choose what fields are desired as output.

 -Dstreamxy      Calculates the Lagrangian barotropic stream function.
 -Dstreamv       Calculates the Lagrangian vertical stream function as a function of model
                 level
 -Dstreamr       —————————————– " ——————————- temperature, salinity/
                 humidity, density/pressure.
    (The the latter two options have to be used with:)
 -Dtempsalt      Calculates the temperature, salinity and density for the ocean
                 and temperature, humidity and pressure for the atmosphere
 -Ddensity       Calculates only the density along the trajectory.
 -Dtracer        Stores trajectory particle positions as a simulated tracer

Non-Passive Advection Options

 -Dsediment       Sediment code developed for RCO
 -Dtwodim         Two-dimensional trajectories, which do not change depth. The vertical velocity
                  is simply set to zero.
 -Dturb           Adds parameterised sub-grid turbulent velocities u’ and v’
 -Ddiffusion      Adds a diffusive random position to the trajectory

It is presently a mystery as to whether any of these will work with our GOMEX grid/fields.

Implementing a New Grid

Page 43 of the Tracmass PDF document gives sketchy details on how to adapt the code to a grid that doesn't presently exist among the Tracmass choices. Start by making a new subdirectory in the projects subdirectory, i.e.

/raid2/GOMEX/lagrange/tracmass/tracmass/tracsvn/projects

which will in this case be gomex, i.e.

mkdir /raid2/GOMEX/lagrange/tracmass/tracmass/tracsvn/projects/gomex

Then copy the contents of one of the existing projects subdirectories into the new gomex subdirectory. Since we're using NetCDF files, we'll choose a subdirectory containing a model that also uses NetCDF files. These include eccoSOSE, gomoos, hycom, jplSCB, multcol, ncarCSM and orc. The ncarCSM example seems the simplest pattern to follow so we will copy that via:

cd /raid2/GOMEX/lagrange/tracmass/tracmass/tracsvn/projects
cp -R ncarCSM/* gomex
cd gomex
mv ncarCSM_run.in gomex_run.in
mv ncarCSM_grid.in gomex_grid.in

Note: The hycom directory file looks to be fairly useful, too.

The gomex directory now contains the following files:

gomex_grid.in
gomex_run.in
Makefile
readfield.f95
setupgrid.f95

The file readfield.f95 (itself called by loop.f95 in the src subdirectory) does all the heavy lifting concerning the input files. It calls subroutines dealing with NetCDF files that are buried over in the file:

/raid2/GOMEX/lagrange/tracmass/tracmass/tracsvn/src/getfile.f95

The getfile.f95 file contains generic calls designed to read 1-, 2- and 3-D fields, i.e.

get1DfieldNC, get2DfieldNC, get3DfieldNC

In the example readfield.f95 file we copied from the ncarCSM subdirectory these calls are accessed via the lines:

uvel = get3DfieldNC(trim(fieldFile), 'UVEL') / 100
vvel = get3DfieldNC(trim(fieldFile), 'VVEL') / 100

Input Files

The Tracmass code has two required input files, test_grid.in and test_run.in.

Grid File

The grid file defines the parameters of the input grid.

The pertinent variables for our situation are:

• IMT - 401 - grid points in x or longitudinal direction
• JMT - 261 - grid points in y or latitudinal direction
• KM - 1 - grid points in vertical direction
• LBT - 1 - stream function separation dimension
• NEND - 2 - number of end sections defined (NEND = LBT + 1)
• NTRACMAX - ? - no. of drifters times time steps times ?
• ngcm - 24 - hours between datasets
• iter - ? (5) - iterations between datasets
• intmax - 365 - maximum length of RCO fields
• fieldsPerFile - 365 - number of time fields per file
• yearmin - 0
• yearmax - 0
• baseSec - 0
• baseMin - 0 - reference basetime for velocity field
• baseHour - 0
• baseDay - 1
• baseMon - 1
• baseYear - 2010
• arcscale - ? - original arc meters -> km

An example of the input grid file located at:

/raid2/GOMEX/lagrange/tracmass/tracmass/tracsvn/projects/ncarCSM/ncarCSM_grid.in

is:

&INITGRIDVER
  ! === Used to check if this file has the correct syntax
  gridVerNum  = 2/
  !
  !
 &INITGRIDDESC
  !
  GCMname     = 'Community Climate System Model CCSM/POP',
  GCMsource   = 'http://www.ccsm.ucar.edu/models/ccsm3.0/pop/',
  gridName    = 'Scott Doney CSMBGC 3deg runs',
  gridSource  = 'http://www.whoi.edu/sbl/liteSite.do?litesiteid=23412&articleId=35629',
  gridDesc    = '3.5 deg run for the global Ocean.'/
  !
  !
$INITGRIDGRID
  !
  IMT      = 116,
  JMT      = 100,
  KM       = 25,
  LBT      = 1,
  NEND     = 6/                   ! NEND = LBT +1
  !
  !
$INITGRIDNTRAC
  !
  NTRACMAX=7000000/               ! 7*1000*1000 
  !
  !
$INITGRIDTIME
  !
  ngcm          = 24,             ! (5*24) hours between GCM datasets
  iter          = 5,              ! iteration between two gcm data sets
  intmax        = 365,            ! maximum length of RCO fields
  fieldsPerFile = 31/
  !
  !
$INITGRIDDATE
  !
  yearmin        = 0,
  yearmax        = 0,
  ! === Reference basetime for the velocity field.
  baseSec        = 0,
  baseMin        = 0,
  baseHour       = 0,
  baseDay        = 1,
  baseMon        = 1,
  baseYear       = 2004/
  !
  !
$INITGRIDARC
   arcscale      = 0.001/         ! orig arc meters -> km
!  arcscale      = 0.00001/       ! orig arc meters -> 100 km  
                                  ! (occ66 || ifs)

Run File

The run file contains values needed for various parameters for the trajectory run.

The pertinent variables for our purposes are:

• subGrid - 0 - 0 for using the full available grid
• intmin - 1 - starting value for initial dataset
• intspin - 1 - trajectory release period (in timesteps)
• intrun - 365 - number of timesteps to be used for run
• intstep - 1 - positive forward, negative backward
• startHour - 0 - starting hour for run
• startDay - 1 - starting day for run
• startMon - 1 - starting month for run
• startYear - 2010 - starting year for run
• ncoor - 1 - output in lon/lat coordinates
• kriva - 1 - write at time intervals of datasets
• inDataDir - ? - directory containing input fields
• outDataDir - ? - directory for saving output fields
• outDataFile - ? - name of present trajectory run for file naming
• nff - 3 - follow trajectories both ways
• isec - 3 - seed particles in horizontal plane
• idir - 0 - follow both directions
• nqua - 1 - constant number of particles in all boxes
• partQuant - 1 - particles per grid cell
• SeedType - 2 - use a list to define which cells to seed (1 = seed an area defined by ist, jst, kst)

seedType

The seedType variable offers three methods for seeding particles:

 1 = Seed an area defined by ist, jst, and kst.
 2 = Use a list to define which cells to seed.
 3 = Use a 2-D mask file.

Since we want to compare the trajectories induced by our various GOMEX model horizontal velocity fields, we will want to use option 2 and create a list of which grid cells we want to seed based on the data in the drifting-buoys-*.nc files in the /raid2/GOMEX/lagrange/drifters directory.

All of the examples in the projects directory use option 1, so we must dig into the source code file src/init_seed.f95 to figure out how to use option 2. The seedlist section of that file is listed below and tells us that we need to specify a variable varSeedFile as 2 to seed particles from a given listfile, and a variable seedFile giving the name of the the given listfile. We further find that the listfile has the format (6I6), with the six values being:

1. i - i or x or longitudinal coordinate
2. j - j or y or latitudinal coordiate
3. k - k or vertical coordinate, 1 = surface
4. idir - direction to follow, 0 = positive and negative
5. isec - particle seeding direction, 3 = horizontal plane
6. -1

Thus an example seedfile - e.g. drifters_seed.txt - will look like:

    27    42     1     0     3    -1
    85   123     1     0     3    -1
...

The $INITRUNSEED section of the gomex_run.in file will then be:

$INITRUNSEED
  nff      = 3,
  isec     = 3,
  idir     = 0,
  nqua     = 1,
  partQuant= 1,    !(particles/gridcell or m3s-1/particle or m3/particle)
  seedType = 2,
  varSeedFile = 2,
  seedFile = "/raid2/GOMEX/lagrange/drifters/drifters_seed.txt"/

The task here is to write a program to convert the initial lon/lat coordinates in each of the drifting-buoys-*.nc files into their corresponding i and j values in the standard GOMEX grid.

Here is the program written to do this. It takes as input all the drifting-buoys-*.nc files in sequence, and outputs two files:

• drifters_seed.txt - a file formatted for use with the TRACMASS program
• test_drifters.txt - a file containing a human readable summary about the contents of drifting-buoys-*.nc

The contents of the test_drifters.txt file are:

File Listing - test_drifters.txt

                    GOMEX       closest             
  1st      1st      i    j    GOMEX grid          days since 1950-01-01                 calendar date
  lon      lat     lon  lat    lon    lat       start           finish               start               finish
-79.899   20.516   362   50   -79.9   20.5   22141.7826389   22157.9520833   2010-08-15 18:47:00   2010-08-31 22:51:00
-86.43   20.559   231   51   -86.45   20.55   22078.975   22157.9548611   2010-06-13 23:24:00   2010-08-31 22:55:00
-83.795   20.555   284   51   -83.8   20.55   22101.6020833   22120.3861111   2010-07-06 14:27:00   2010-07-25 09:16:00
-79.745   29.22   365   224   -79.75   29.2   22085.3666667   22086.7048611   2010-06-20 08:48:00   2010-06-21 16:55:00
-83.796   24.872   284   137   -83.8   24.85   22101.5270833   22111.0104167   2010-07-06 12:39:00   2010-07-16 00:15:00
-86.833   24.737   223   135   -86.85   24.75   22105.5013889   22157.9993056   2010-07-10 12:02:00   2010-08-31 23:59:00
-88.111   27.756   198   195   -88.1   27.75   22112.0756944   22157.95625   2010-07-17 01:49:00   2010-08-31 22:57:00
-87.653   26.515   207   170   -87.65   26.5   22111.5965278   22141.8576389   2010-07-16 14:19:00   2010-08-15 20:35:00
-86.43   24.105   231   122   -86.45   24.1   22103.625   22157.9548611   2010-07-08 15:00:00   2010-08-31 22:55:00
-85.278   25.331   254   147   -85.3   25.35   22106.3604167   22157.9993056   2010-07-11 08:39:00   2010-08-31 23:59:00
-83.016   25.059   300   141   -83.0   25.05   22101.1534722   22157.9527778   2010-07-06 03:41:00   2010-08-31 22:52:00
-86.846   24.756   223   135   -86.85   24.75   22105.4993056   22157.9541667   2010-07-10 11:59:00   2010-08-31 22:54:00
-87.58   26.514   208   170   -87.6   26.5   22111.625   22142.3756944   2010-07-16 15:00:00   2010-08-16 09:01:00
-85.964   25.856   241   157   -85.95   25.85   22104.84375   22157.9340278   2010-07-09 20:15:00   2010-08-31 22:25:00
-85.964   25.857   241   157   -85.95   25.85   22104.8506944   22157.9548611   2010-07-09 20:25:00   2010-08-31 22:55:00
-84.696   25.526   266   151   -84.7   25.55   22106.75   22157.9513889   2010-07-11 18:00:00   2010-08-31 22:50:00
-83.035   25.013   299   140   -83.05   25.0   22101.2625   22108.9194444   2010-07-06 06:18:00   2010-07-13 22:04:00
-84.186   26.229   276   165   -84.2   26.25   22074.75   22157.9527778   2010-06-09 18:00:00   2010-08-31 22:52:00
-84.688   25.528   266   151   -84.7   25.55   22106.7583333   22157.9534722   2010-07-11 18:12:00   2010-08-31 22:53:00
-85.312   24.492   254   130   -85.3   24.5   22102.5875   22116.0798611   2010-07-07 14:06:00   2010-07-21 01:55:00
-85.754   24.88   245   138   -85.75   24.9   22104.1125   22157.9541667   2010-07-09 02:42:00   2010-08-31 22:54:00
-88.0   28.001   200   200   -88.0   28.0   22054.2013889   22063.2430556   2010-05-20 04:50:00   2010-05-29 05:50:00
-84.591   24.824   268   136   -84.6   24.8   22102.3229167   22153.2930556   2010-07-07 07:45:00   2010-08-27 07:02:00
-85.28   25.33   254   147   -85.3   25.35   22106.3673611   22157.9541667   2010-07-11 08:49:00   2010-08-31 22:54:00
-88.112   27.756   198   195   -88.1   27.75   22112.075   22157.9347222   2010-07-17 01:48:00   2010-08-31 22:26:00
-88.444   28.668   191   213   -88.45   28.65   22053.7847222   22086.8680556   2010-05-19 18:50:00   2010-06-21 20:50:00
-86.031   26.684   239   174   -86.05   26.7   22073.9465278   22157.9548611   2010-06-08 22:43:00   2010-08-31 22:55:00
-87.208   28.025   216   200   -87.2   28.0   22073.125   22157.9527778   2010-06-08 03:00:00   2010-08-31 22:52:00
-84.02   26.161   280   163   -84.0   26.15   22074.7930556   22157.95   2010-06-09 19:02:00   2010-08-31 22:48:00
-88.019   28.631   200   213   -88.0   28.65   22072.8090278   22088.5055556   2010-06-07 19:25:00   2010-06-23 12:08:00
-87.98   28.637   200   213   -88.0   28.65   22072.8652778   22097.6222222   2010-06-07 20:46:00   2010-07-02 14:56:00
-84.168   26.229   277   165   -84.15   26.25   22074.7944444   22144.7680556   2010-06-09 19:04:00   2010-08-18 18:26:00
-86.865   27.751   223   195   -86.85   27.75   22073.3145833   22157.9555556   2010-06-08 07:33:00   2010-08-31 22:56:00
-84.029   26.167   279   163   -84.05   26.15   22074.7798611   22157.9520833   2010-06-09 18:43:00   2010-08-31 22:51:00
-85.291   24.473   254   129   -85.3   24.45   22102.5881944   22115.4756944   2010-07-07 14:07:00   2010-07-20 11:25:00
-83.8   24.881   284   138   -83.8   24.9   22101.6743056   22110.9520833   2010-07-06 16:11:00   2010-07-15 22:51:00
-89.586   27.784   168   196   -89.6   27.8   22099.7361111   22132.7708333   2010-07-04 17:40:00   2010-08-06 18:30:00
-89.851   28.35   163   207   -89.85   28.35   22098.3402778   22126.6875   2010-07-03 08:10:00   2010-07-31 16:30:00
-86.039   26.702   239   174   -86.05   26.7   22073.94375   22157.9534722   2010-06-08 22:39:00   2010-08-31 22:53:00
-86.513   24.173   230   123   -86.5   24.15   22103.5729167   22149.3256944   2010-07-08 13:45:00   2010-08-23 07:49:00
-89.982   28.144   160   203   -90.0   28.15   22099.2638889   22133.6875   2010-07-04 06:20:00   2010-08-07 16:30:00
-83.008   25.041   300   141   -83.0   25.05   22101.15625   22108.7541667   2010-07-06 03:45:00   2010-07-13 18:06:00
-89.602   27.852   168   197   -89.6   27.85   22098.4791667   22151.7708333   2010-07-03 11:30:00   2010-08-25 18:30:00
-89.517   27.601   170   192   -89.5   27.6   22098.5486111   22145.0416667   2010-07-03 13:10:00   2010-08-19 01:00:00
-87.226   27.997   215   200   -87.25   28.0   22073.3125   22157.93125   2010-06-08 07:30:00   2010-08-31 22:21:00
-84.594   24.702   268   134   -84.6   24.7   22102.0652778   22134.6145833   2010-07-07 01:34:00   2010-08-08 14:45:00
-89.719   28.051   166   201   -89.7   28.05   22098.4166667   22130.375   2010-07-03 10:00:00   2010-08-04 09:00:00
-87.593   27.198   208   184   -87.6   27.2   22060.4027778   22115.6666667   2010-05-26 09:40:00   2010-07-20 16:00:00
-88.404   28.812   192   216   -88.4   28.8   22055.2986111   22077.9861111   2010-05-21 07:10:00   2010-06-12 23:40:00
-88.293   28.71   194   214   -88.3   28.7   22054.3958333   22092.7847222   2010-05-20 09:30:00   2010-06-27 18:50:00
-88.6   29.201   188   224   -88.6   29.2   22055.1666667   22107.3888889   2010-05-21 04:00:00   2010-07-12 09:20:00
-87.993   26.506   200   170   -88.0   26.5   22059.6458333   22103.3194444   2010-05-25 15:30:00   2010-07-08 07:40:00
-87.413   26.665   212   173   -87.4   26.65   22060.2638889   22097.0416667   2010-05-26 06:20:00   2010-07-02 01:00:00
-87.008   26.003   220   160   -87.0   26.0   22059.8888889   22116.6458333   2010-05-25 21:20:00   2010-07-21 15:30:00
-87.199   26.401   216   168   -87.2   26.4   22060.1805556   22091.0138889   2010-05-26 04:20:00   2010-06-26 00:20:00
-88.186   28.817   196   216   -88.2   28.8   22056.3611111   22075.8541667   2010-05-22 08:40:00   2010-06-10 20:30:00
-87.709   27.578   206   192   -87.7   27.6   22060.4861111   22094.0208333   2010-05-26 11:40:00   2010-06-29 00:30:00

Program Listing - drifters.py

#!/usr/bin/python2.5
# encoding: utf-8
#
#  Reads drifting-buoys-*.nc files, extracts the initial lon/lat value, and converts that to the
#  corresponding i,j values on the standard GOMEX grid.

"""
drifters.py
Created on 2011-10-27

Copyright (c) 2011 Steven K. Baum

Version 1:  
"""

import os, sys, fnmatch
import scipy,numpy,netCDF4
import time,datetime
#import matplotlib
#import octant.tools.rot2d
from netCDF4 import num2date, date2num
#from mpl_toolkits.basemap import Basemap
#from numpy import *
#from time import strftime
import numpy as np

timebase = "days since 1950-01-01 00:00:00 UTC"

nlon = 401
nlat = 261

driftdir = "/raid2/GOMEX/lagrange/drifters"

#  Extract all the drifting*.nc filenames from the directory.
dfiles = os.listdir(driftdir)
driftfiles = fnmatch.filter(dfiles,'drifting*')
#driftfiles[0] = driftfiles[0]

#  Create lons/lats for GOMEX standard grid.
lon_gomex = np.arange(-98,-77.95,.05)
lat_gomex = np.arange(18,31.05,.05)

#  In the standard GOMEX grid:
#     longitude[0] = -98.0       longitude[400] = -78.0
#     latitude[0] = 18.0         latitude[260] = 31.0
#     resolution = 0.05 deg. both ways

dfile = open("drifters_seed.txt","w")
ofile = open("test_drifters.txt","w")
tp = "   "

for file in driftfiles:

        ncfile = driftdir + "/" + file
        inf = netCDF4.Dataset(file,'r+')
        lon0 = inf.variables['LONGITUDE'][0]
        lat0 = inf.variables['LATITUDE'][0]
        juld = inf.variables['JULD'][:]
        juld0 = juld[0]
        last = juld.size - 1
        juldZ = juld[last]
        datefirst = num2date(juld0,timebase)
        datelast = num2date(juldZ,timebase)

#  Find the closest GOMEX grid i to the first longitude.
        loncnt = 0
        for lon in lon_gomex:
                if lon > lon0:
                        lonsti = loncnt
                        d1 = lon - lon0
                        d2 = lon0 - lon_gomex[lonsti-1]
                        if d1 > d2:
                                lonsti = loncnt - 1
                        break
                else:
                        loncnt = loncnt + 1

#  Find the close GOMEX grid j to the first latitude.
        latcnt = 0
        for lat in lat_gomex:
                if lat > lat0:
                        latstj = latcnt
                        d1 = lat - lat0
                        d2 = lat0 - lat_gomex[latstj-1]
                        if d1 > d2:
                                latstj = latcnt - 1
                        break
                else:
                        latcnt = latcnt + 1

        inf.close()
        outstr = str(lon0) + tp + str(lat0) + tp + str(lonsti) + tp + str(latstj) + tp + str(lon_gomex[lonsti]) + tp + str(lat_gomex[latstj]) + tp + str(juld0) + tp + str(juldZ) + tp + str(datefirst) + tp + str(datelast) + "\n"
        ofile.write(outstr)

        print >>dfile, repr(lonsti).rjust(6),repr(latstj).rjust(5),repr(1).rjust(5),repr(0).rjust(5),repr(3).rjust(5),repr(-1).rjust(5)

ofile.close()
dfile.close()

Listing N - seedType definition in init_seed.f95

  case (2) ! === seedlist method ===   
     print *,'------------------------------------------------------'
     if (varSeedFile == 1) then
        fileStamp='/seed00000000.asc'
        write (fileStamp(6:13),'(i8.8)') intstart
        fullSeedFile=trim(seedDir) // trim(fileStamp)
        print *,'===  Particles are seeded from a dynamic listfile  ==='
     else
        fullSeedFile=trim(seedFile)
        print *,'===  Particles are seeded from a given listfile    ==='
     end if
     inquire(FILE = fullSeedFile, exist=fileexists)
     chFile: if (fileexists) then
        ijkMax=0
        open(unit=34,file=fullSeedFile, ACCESS = 'SEQUENTIAL', &
             FORM = 'FORMATTED', ACTION = 'READ')
        findRecl: do
           read (unit=34, fmt=ijkform,iostat=filestat)
           if (filestat < 0) exit findRecl
           ijkMax=ijkMax+1
        end do findRecl
        allocate (ijkst(ijkMax,6))

        rewind(34)
        read_ijkst: do ijk=1,ijkMax
           read (unit=34, fmt=ijkform) ijkst(ijk,1), ijkst(ijk,2), &
                ijkst(ijk,3),  ijkst(ijk,4),  ijkst(ijk,5), ijkst(ijk,6)


           ijkst(ijk,4)=idir
           ijkst(ijk,5)=isec
           ijkst(ijk,6)=-1

        end do read_ijkst

        print *,'File name    : '//trim(fullSeedFile)
        print *,'Seed size    : ', ijkMax
     else
        print *,'======================================================'
        print *,'*** ERROR!                                         ***'
        print *,'*** Seed files does not exisit                     ***'
        print *,'File name    : '//trim(fullSeedFile)
        print *,'*** Run terminated.                                ***'
        stop
     end if chFile

An example of the input run file located at:

/raid2/GOMEX/lagrange/tracmass/tracmass/tracsvn/projects/ncarCSM/ncarCSM_run.in

is:

&INITRUNVERSION
  ! === Used to check if this file has the correct syntax
  runVerNum   = 1/
  !
  !
&INITRUNDESC
  !
  caseName     = 'CSMBGCGlobal',
  caseDesc     = 'Testing CSM fields for fun and entertainment.'/
  !
  !
&INITRUNGRID
  !
  !==subGrid:       0 = Use full grid.
  !==               1 = Define subGrid in this file.
  !==               2 = Define subGrid with sep. file and subGridID.
  !==               3 = Define subGrid with sep. file and MPI.
  subGrid      = 0,
  !
  ! === Used if SubGrid = 1
  subGridImin  = 1,
  subGridImax  = 116,
  subGridJmin  = 1,
  subGridJmax  = 100,
  ! === Used if SubGrid = 2 or 3
  SubGridFile  = '/Users/bror/svn/orm/grd/templ.asc',
  ! === Used if SubGrid = 2
  subGridID    = 1/
  !
  !
&INITRUNTIME
  !
  ! === Startval for initial dataset          ===
  intmin       = 50,
  ! === Trajectory release period (timesteps) ===
  intspin      = 1,
  ! === Number of timesteps for the run.      ===
  intrun       = 2000,
  ! === Pos if forward  neg if backward       ===
  intstep      = 1/
  !
  !
&INITRUNDATE
  !
  ! === Start time for this run
  startHour         = 0,
  startDay          = 1,
  startMon          = 3,
  startYear         = 2004/
  !
  !
&INITRUNWRITE
  !
  !==ncoor:         0 = output in model coordinates
  !==               1 = output in long/lat coordinates
  ncoor   = 0,      
  !
  !==kriva:         0 = no writing
  !==               1 = write at time intervals of gcm datasets
  !==               2 = write each grid-crossing and timne change
  !==               3 = write at each iteration (all the time)
  !==               4 = write only start and end positions
  !==               5 = write at chosen intervals
  kriva    = 1,
  !
  ! === Directory where input fields are stored
  inDataDir    = '/projData/CSMBGC/Ar_run/',
  !
  ! === Directory where output files are saved
  !outDataDir   = '/Users/doos/data/orc12/',
  outDataDir   = '/Users/bror/ormOut/',
  !
  ! === name of current trajectory run
  outDataFile     = 'csmTest'/
  ! === Name of rerun file
  ! namep='op.fw.n0'
  !
  !
$INITRUNSEED
  !
  !==nff:           1 = Follow trajectories forward
  !==               2 = Follow trajectories backward
  !==               3 = Follow trajectories both ways.
  nff      = 3,
  !==isec:          1 = Seed particles meridional(y-z)
  !==               2 = Seed particles zonal(x-z)
  !==               3 = Seed particles horiz(x-y)
  !==               4 = Seed particles in the middle of T-box
  isec     = 3,
  !==idir:          1 = follow positive direction (eastward/northward)  
  !==              -1 = follow negative direction (westward/southward)
  !==               0 = both directions
  idir     = 0,
  !
  ! number of trajectories can be set by
  ! nqua:           1 = constant number of particles in all boxes
  !                     (partQuant in # particles / gridcell)
  !                 2 = Each particle reflects water transport at seeding. 
  !                     (partQuant in m3s-1. per particle)
  !                 3 = Each particle reflects water volume at seeding.
  !                     (partQuant in m3 per particle)
  nqua     = 1, 
  partQuant= 10,    !(particles/gridcell or m3s-1/particle or m3/particle)
  !
  ! === initial directions all in MODEL COORDINATES ===
  ! Method for seeding particles.
  ! seedType:       1 = Seed an area defined by ist, jst, and kst.
  !                 2 = Use a list to define which cells to seed.
  !                 3 = Use a 2-D mask file.
  SeedType = 1,
  !
  ! ===   ===   ===   ===   ===
  ! If seedType = 1, define area where particles are seeded.
  ! -1 indicates max value in grid.
  !
  ist1     = 50,                 
  ist2     = 70,                 
  jst1     = 50,
  jst2     = 70,
  kst1     =  1,                 
  kst2     =  -1/
  !
  !
$INITRUNTEMPSALT
  !
  ! === Control trajectories by salt and temp ===
  ! ===  (active only with option tempsalt)   ===
  !
  !==Starting a trajectory
  tmin0      =  -50.,
  tmax0      =  400.,
  smin0      = -500.,
  smax0      =  400.,
  rmin0      = -100.,
  rmax0      =  500.,
  ! 
  !==Ending a trajectory
  tmine      =  -50.
  tmaxe      =  400.,
  smine      = -150.,
  smaxe      =  500.,
  rmine      = -100.,
  rmaxe      =  500./
  !
  !
$INITRUNEND     
  ienw(1) =   10,       
  iene(1) =   90,       
  jens(1) =   10,       
  jenn(1) =   90,       
  timax =  36500./  ! maximum time length of a trajectory in days

Plotting Trajectories

An NCL trajectory plotting example can be found at:

http://www.ncl.ucar.edu/Applications/traj.shtml

described as:

A simple trajectory plot. Each trajectory is a different color for clarity, and every fourth time step is marked with a circle. The start of each trajectory is marked with a green circle.

The source code for the example is at:

http://www.ncl.ucar.edu/Applications/Scripts/traj_1.ncl

and is:

;************************************************* 
; traj_1.ncl
;*************************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
;**************************************************
begin
;*************************************
; read in data
;************************************
; note this trajectory file contains data in the form of
; 9 variables x 131 timesteps x 100 trajectories
   ntime = 131

   big = fbinrecread("traj.bin",0,(/9,ntime,100/),"float")
   time= fbinrecread("traj.bin",1,(/ntime/),"float")
;********************************************
   wks = gsn_open_wks("ps","traj")               ; open workstation
   res            = True                          ; map resources
   res@gsnDraw    = False                         ; don't draw
   res@gsnFrame   = False                         ; don't advance frame
   res@vpWidthF   = 0.80                          ; make map bigger
   res@vpHeightF  = 0.80
   res@mpMaxLatF  = -20                           ; select subregion
   res@mpMinLatF  = -60 
   res@mpMinLonF  = -75
   res@mpMaxLonF  = -25

   res@tiMainString    = "Example of a trajectory plot"  ; title
   res@gsnCenterString = "markers every 4th time step"   ; center string
  
   map = gsn_csm_map_ce(wks,res)                         ; create map
;*********************************************
; trajectory parameters
;*********************************************
   xpt = new(ntime,float)            ; allocate memory
   ypt = new(ntime,float)
   traj = (/1,10,53,67,80/)          ; choose which trajectories to plot
;*********************************************
; some plot parameters
;*********************************************
   pres                  = True               ; polyline resources
   pres@gsLineThicknessF = 2.0                ; line thickness
   colors= (/"red","blue","dark green","grey","magenta"/) ; line color

   mres  = True                               ; marker resources
   first = True                               ; start of traj resources
;********************************
   c=0                                        ; counter variable
   do i = 0,dimsizes(traj)-1                  ; loop through chosen traj
      ypt = big(2,:,traj(i))                  ; extract lat from whole array
      xpt = big(1,:,traj(i))                  ; extract lon from whole array

      pres@gsLineColor = colors(c)            ; line color for this traj
      c=c+1                                   ; advance counter
      gsn_polyline(wks,map,xpt,ypt,pres)      ; draw the traj

; add markers to the trajectories

      mres@gsMarkerIndex  = 16           ; marker style (circle)
      mres@gsMarkerSizeF  = 4.0          ; marker size
      mres@gsMarkerColor  = "black"      ; maker color
      gsn_polymarker(wks,map,xpt(::4),ypt(::4),mres) ; draw every 4th marker

; create a unique marker to indicate the start of the trajectory

      first@gsMarkerSizeF  = 8.0         ; marker size
      first@gsMarkerColor  = "green"     ; marker color

      gsn_polymarker(wks,map,xpt(0),ypt(0),first) ; draw start of traj
      delete(first@gsMarkerColor)
      delete(first@gsMarkerSizeF)
   end do 
   draw(map)                                           
   frame(wks)                                         
  
end