GAMAP bugs and fixes

From Geos-chem
Jump to: navigation, search

GAMAP is now obsolete. We recommend using GCPy for analyzing output from recent GEOS-Chem versions. But we will preserve the GAMAP wiki documentation for reference.



On this page we list various bugs that have been encountered in GAMAP and how to fix them. For information about general GAMAP usage, please see our GAMAP tips and tricks wiki page.

Outstanding issues

Bugs in gc_combine_nd49.pro

Cannot combine more than 32767 data blocks

Mike Barkley wrote:

I've found a bug in version v2.12 of gamap with the program: gc_combine_nd49.pro.
The final 4th output dimension is calculated using the nblocks multiplied by the nfiles, and this value is then assigned to the dim array (which is an INT array; max value 32767).
If however, you are combing several years of data it is possible that nblocks * nfiles can exceed 32767 which cannot be assigned into the dim array. To get around this, one can simply convert the dim into a LONG array as soon as it is created (see extract of code below).
   ;; --------- From 1st file only ----------------
        IF First THEN BEGIN

           ;; Check availability of tracer/category records, and get nb
           ;; of daily timesteps
           nblocks = n_elements( datainfo )

    ....
  
           ;; Get MODELINFO and GRIDINFO structures, dimension and offset
           GetModelAndGridInfo, DataInfo[0], ModelInfo, GridInfo

           ; Previous code
           ; dim    =  dataInfo[0].dim
           ;  *** Simple correction ***
           dim    = LONG( dataInfo[0].dim )
           offset = dataInfo[0].first 

Philippe Le Sager wrote:

Thanks Mike, nice catch. It solves the problem locally: you can correctly plot and manipulate data now. However there may be some issue when you save them, since the call to ctm_make_datainfo.pro will try to assign LONG to an INT array (l. 441):
   datainfo.dim = ( dim > 1 )

So, we need another fix in create3dhstru.pro. We need to replace l. 94:

          Dim        : [ 0,0,0,0 ], $  ; dimension of data array (I,J,L,time)

with

          Dim        : [ 0L,0L,0L,0L ], $  ; dimension of data array (I,J,L,time)

Ultimately, a nice fix is to force all integers to 32-bit with:

      compile_opt defint32

But you have to put it in every routine....there must be a way to set that option globally, but I have not found it yet.

--Bob Y. 09:16, 2 December 2008 (EST)

Incorrect results when using /DMAX or /DAVG options

A minor bug could give incorrect results when either /DMAX or /DAVG was used, and timeseries from more than one month were processed at once. The output OUTTIME vector (if requested) was not correct, and if a smaller TIME span than the original series was requested, you could end up with a wrong series. The fix will be available in the next GAMAP release (v2.13). Meanwhile you can get it here.

--phs 14:55, 15 August 2008 (EDT)

I cannot get a correct GIF or JPEG image from the screen

In the current version of GAMAP (v2.12), we use a version of TVREAD to get screen dump that has a bug for 16-bit displays. If you are using screen2gif.pro or tvread.pro directly, and your images do not look like the screen output, you are probably using a 16-bit display. You can:

  • switch to true color 32-bit display if available (go to display in the control panel),
  • or update tvread.pro

--phs 16:39, 2 September 2008 (EDT)

Using PS-PTOP in GAMAP with GEOS-Chem v8+ outputs

GEOS-Chem does not write the category PS-PTOP into the tracerinfo.dat anymore (see PS-PTOP vs. PEDGE). However some routines still try to get the pressure surface with that category name. You can simply replace PS-TOP by PEDGE-$ in your program for them to work. To be able to use older GEOS-Chem output and new ones, you may use a coding similar to the one I did for ctm_column_du.pro, which you can find here.

Note that to read pressure surface from time series diagnostic (ND48-51), you need to keep PS-PTOP in your IDL routines and replace PEDGE-$ by PS-PTOP in the tracerinfo.dat files created by GEOS-Chem. This is a bug, and we expect to get rid of all references to PS-PTOP in a future release of GAMAP.

--phs 16:25, 2 September 2008 (EDT)

I cannot read a binary file larger than about 2.1 GB

Paul Hezel wrote:

I am trying to read large ctm.bpch (1 year simulation) files via GAMAP and associated routines. I've found that it fails because the filepos datatype is only a LONG, (eg. in internals/create3dhstru.pro, or in internals/ctm_read_data.pro). The following posts seem to address this issue:
If I changed filepos to either LONG64 or even ulong, in create3dhstru.pro, would that capture most of the times it runs into this? or would i have to trace and amend all the places where that variable (i assume mostly filepos) were cast from or back to a LONG datatype? Are there other issues I may be missing?

Bob Yantosca replied:

Changing the file pointer to a LONG64 would probably be the solution. But we also have to be careful that once we pull on the thread, we don't unravel the whole sweater. :-) There could conceivably be places in several internal routines that would need modification...but won't know until we look.
Also -- one of the reasons we haven't looked at this is because you might need the 64-bit version of IDL to get the LONG64 data type. Not all of our GAMAP users may have 64-bit IDL.

--Bob Y. 16:27, 26 November 2008 (EST)

Backwards compatibility problems in myct.pro with IDL 5

Nicolas Bousserez wrote:

I got into trouble when using an old version of IDL (v5.6) together with the GAMAP tools. The errors seem come from the symbol ~ in the new myct.pro routine:
   Reverse_Colors  = ~Reverse_Colors
                       ^
   % Illegal character in program text.
I also experienced some colors problem when using IDL v5.6 instead of more recent version like 6 or 7. For example I got black and white colors only for the IDL 5.6 version instead of colors for the idl 6 or 7 version when plotting.
The problem is that for the moment we only have a limited number of IDL 7 licenses so that it would be very practical for me to be able to use the older versions of IDL together with the new GAMAP tools.

Philippe Le Sager replied:

The ~ (logical not) is not part of older IDL, but you can get around it with the NOT (bitwise negation). You just need to add 2. For example,
   top = not below + 2  ; same as: ~below (in IDL6.x)
You can quickly test it in IDL 7. As for your missing colors, the first advice is to be sure you have set:
   device, decomposed=0
After that, if you still have problem double check that your idl_startup.pro file is read.

--Bob Y. 16:40, 26 November 2008 (EST)

GAMAP plot titles showing "incorrect" altitudes for plots

Ray Nassar wrote:

The altitudes of levels that GAMAP uses to label plots does not agree with the altitudes shown in the GEOS-Chem online manual, for GEOS-4 and possibly other met fields. This is a result of GAMAP using a different surface pressure to translate the sigma level into an altitude than used in the manual. One can change GAMAP to get the same altitudes with the following edits to the code in the gamap_util directory:

In CTM_GRID, add (somewhere pretty close to the beginning):

  ;====================================================================
  ; Define surface pressure
  ;====================================================================
       PSURF = 1000
       print, "PSURF =", psurf

To get a second decimal place for the altitude, edit CTM_LABEL to show:

     AltStr = StrTrim( String( Alt[0], Format='(f6.2)' ), 2 ) + ' km'

This will give the same altitudes as in the manual (to 2 decimal places) but it disables calls with alternate values of "psurf". There is likely a better solution for future versions of GAMAP, but this should be sufficient for most purposes. Furthermore, 1013 hPa would probably be the best value to use for a default "psurf" in the long term since it is the standard sea level pressure (1 atm).

Resolved issues

IDLDE does not work properly

If you are running IDL and Firefox in a Linux environment then please take note. You may have to disable Firefox in order to get the IDL Development Environment (IDLDE) to work properly with IDL 7.

The issue is that the 32-bit version of Firefox (esp. Firefox 3) may cause a conflict with a Java library that the IDLDE needs to use. Disabling the browser seems to fix the problem.

--Bob Y. 12:52, 20 October 2008 (EDT)

Corrupted binary file error

Dylan Millet wrote:

Hi guys -- I ran into a problem with bpch_link.pro:
   IDL> bpch_link,'2006*bpch','new.bpch' 
   Now Reading 20060101_2x25.ctm.bpch
   % READU: Corrupted f77 unformatted file detected. Unit: 108, File: 20060101_2x25.ctm.bpch
   % Execution halted at: BPCH_LINK         142 /users/dbm/IDL/gamap_v2.10/file_io/bpch_link.pro
   %                      BPCH_LINK         142 /users/dbm/IDL/gamap_v2.10/file_io/bpch_link.pro
   %                      $MAIN$         IDL>

Bob Yantosca replied:

You may be using an older version of bpch_link.pro. We recently fixed a problem for the little/big endian but that may not be in the std gamap as of yet. The fix is simple. You have to use the SWAP_ENDIAN keyword to all instances of OPEN_FILE:
   ; External functions
   FORWARD_FUNCTION CTM_Grid, CTM_Type, MFindFile, Little_Endian

   ; Open the output file (write as big-endian)
   Open_File, OutFile, Ilun_OUT, /F77, /GET_LUN, /Write, Swap_Endian=Little_Endian()

Specifying SWAP_ENDIAN=LITTLE_ENDIAN() in each call to OPEN_FILE will cause IDL to swap the endian ordering if you are running IDL on a little-endian machine.

NOTE: In GAMAP v2-12, all routines that read binary files now use the SWAP_ENDIAN function, so this type of error should not occur again.

--Bob Y. 10:35, 18 July 2008 (EDT)

Too many logical unit numbers used by GAMAP

Chris Holmes wrote:

A lot of gamap routines seem to use up LUNs for file io. This means that after iterating several read cycles, I can no longer open more files, even though I'm closing them. Here's what's going on:
A lot of gamap2/gamap_util/ routines are opening files with the following command
   OPEN_FILE, File, LUN, /GET_LUN
OPEN_FILE checks whether the LUN variable is undefined (lines 151 and 241) and if so it allocates one (line 241).
But then OPEN_FILE also passes the /GET_LUN keyword to OPENW or OPENR via the _EXTRA keyword (lines 244-249), which allocates another LUN without freeing the first one.
I'm not sure what the best solution is. Probably remove the /GET_LUN keyword in all instances of OPEN_FILE. You could also add a GET_LUN keyword to OPEN_FILE that does nothing except prevent it from passing it on to OPENW and OPENR.
Having resolved that, I still run out of LUNs and I've traced the problem to CTM_OPEN_FILE. Actually there's already a comment about this problem on lines 725-731. So yes I think it is a problem and should be fixed.

Philippe Le Sager wrote:

I looked at the issues you pointed out. OPEN_FILE was indeed automatically using get_lun, while the user could pass /get_lun with _extra. To avoid the double LUN assignment, I think the dummy /get_lun keyword is the best solution (users can put whatever they want in the _extra), although we could be more sophisticated and check on the field names in the extra keyword structure.
CTM_OPEN_FILE was more difficult. I did not find it well written, and it had quite a few coding errors and misleading comments. I decided to simplify it, and get rid of all the "testlun" business. It uses all LUN instead of every other one now. I am testing my version now and it seems to work fine. Feel free to look at it if you want.

This bug has now been resolved in GAMAP v2-12.

--Bob Y. 13:33, 8 September 2008 (EDT)

User Defined Color Table

When upgrading the color tables in GAMAP v2.12, one functionality was lost: you could not define your own color table. An updated myct.pro routine was added to GAMAP v2-13) which will let you set your own RGB vectors.

If you are using an older version of GAMAP, then you can get the updated myct.pro file HERE.

--Bob Y. 16:51, 3 March 2010 (EST)

Bug in BPCH2HDF

Jesse Kenyon reported that there was a bug in the BPCH2HDF routine, which converts binary punch format to HDF format. There was a typo in the code. The lines:

  ; Read all data blocks
  if ( N_Elements( Category ) gt 0 )                           $
     then CTM_Get_Data, DataInfo, DiagN, File=InFile, _EXTRA=e $
     else CTM_Get_Data, DataInfo,        File=InFile, _EXTRA=e

should now be replaced with:

  ; Read all data blocks
  if ( N_Elements( DiagN ) gt 0 )                              $
     then CTM_Get_Data, DataInfo, DiagN, File=InFile, _EXTRA=e $
     else CTM_Get_Data, DataInfo,        File=InFile, _EXTRA=e

This fix was included in GAMAP release v2-13.

--Bob Y. 16:51, 3 March 2010 (EST)

CTM_OVERLAY cannot make log-scale color plots

Chris Holmes reported that the LOG keyword was not being properly passed to the TVMAP routine call within ctm_overlay.pro. Therefore we have to add it manually, as shown below:

  ; Call TVMAP to plot the pixel or contour map, but don't advance to the
  ; next frame.  NOTE: C_COLORS is the vector of color levels for the
  ; filled contour plot.  If not passed, then C_COLORS will return the 
  ; default values generated by internally by TVMAP.  This will be needed 
  ; for overlaying station data points on a filled contour map.
  if ( OverPlot ne 1L ) then begin
     TvMap, Data, Xmid, Ymid,                                    $
        /NoAdvance,        WindowPos=Wp,      MinData=MinData,   $
        MaxData=MaxData,   FContour=FContour, C_Levels=C_Levels, $
        C_Colors=C_Colors, Log=Log,           _EXTRA=e 
  endif

This fix was included in GAMAP release v2-13.

--Bob Y. 16:51, 3 March 2010 (EST)

Bug in isopleth_map.pro

Matt Alvarado wrote:

I found an error on line 314 of isopleth_map.pro in GAMAP v2-12: It should say !D.Table_size, not !D.Tble_size.

This fix was added to GAMAP version v2-13.

--Bob Y. 16:51, 3 March 2010 (EST)